Mercurial > hg > octave-nkf

comparison src/xpow.cc @ 9123:f39b98237d5c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use more robust and less aggressive scalar .^ range optimization
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 16 Apr 2009 07:02:31 +0200
parents 10bed8fbec99
children ab80681c44d9
comparison
equal deleted inserted replaced
9122:8ca06fd9c6ef 9123:f39b98237d5c
690 } 690 }
691 691
692 return result; 692 return result;
693 } 693 }
694 694
695 static inline bool
696 same_sign (double a, double b)
697 {
698 return (a >= 0 && b >= 0) || (a <= 0 && b <= 0);
699 }
700
695 octave_value 701 octave_value
696 elem_xpow (double a, const Range& r) 702 elem_xpow (double a, const Range& r)
697 { 703 {
698 octave_value retval; 704 octave_value retval;
699 705
700 if (r.nelem () <= 0) 706 // Only optimize powers with ranges that are integer and monotonic in
701 retval = Matrix (); 707 // magnitude.
702 else if (a < 0.0 && ! r.all_elements_are_ints ()) 708 if (r.nelem () > 1 && r.all_elements_are_ints ()
703 { 709 && same_sign (r.base (), r.limit ()))
704 ComplexMatrix mat (1, r.nelem ()); 710 {
705 Complex atmp (a); 711 octave_idx_type n = r.nelem ();
706 Complex *pmat = mat.fortran_vec (); 712 Matrix result (1, n);
707 713 if (same_sign (r.base (), r.inc ()))
708 pmat[0] = std::pow (atmp, r.base ()); 714 {
709 Complex mul = std::pow (atmp, r.inc ()); 715 double base = std::pow (a, r.base ());
710 for (octave_idx_type i = 1; i < r.nelem (); i++) 716 double inc = std::pow (a, r.inc ());
711 pmat[i] = pmat[i-1] * mul; 717 result(0) = base;
712 718 for (octave_idx_type i = 1; i < n; i++)
713 retval = mat; 719 result(i) = (base *= inc);
714 } 720 }
715 else 721 else
716 { 722 {
717 Matrix mat (1, r.nelem ()); 723 // Don't use Range::limit () here.
718 double *pmat = mat.fortran_vec (); 724 double limit = std::pow (a, r.base () + (n-1) * r.inc ());
719 725 double inc = std::pow (a, -r.inc ());
720 double base = std::pow (a, r.base ()); 726 result(n-1) = limit;
721 pmat[0] = base; 727 for (octave_idx_type i = n-2; i >= 0; i--)
722 double mul = std::pow (a, r.inc ()); 728 result(i) = (limit *= inc);
723 for (octave_idx_type i = 1; i < r.nelem (); i++) 729 }
724 pmat[i] = (base *= mul); 730
725 731 retval = result;
726 retval = mat; 732 }
727 } 733 else
734 retval = elem_xpow (a, r.matrix_value ());
728 735
729 return retval; 736 return retval;
730 } 737 }
731 738
732 // -*- 3 -*- 739 // -*- 3 -*-
929 octave_value 936 octave_value
930 elem_xpow (const Complex& a, const Range& r) 937 elem_xpow (const Complex& a, const Range& r)
931 { 938 {
932 octave_value retval; 939 octave_value retval;
933 940
934 if (r.nelem () <= 0) 941 // Only optimize powers with ranges that are integer and monotonic in
935 retval = Matrix (); 942 // magnitude.
936 else 943 if (r.nelem () > 1 && r.all_elements_are_ints ()
937 { 944 && same_sign (r.base (), r.limit ()))
938 ComplexMatrix mat (1, r.nelem ()); 945 {
939 Complex *pmat = mat.fortran_vec (); 946 octave_idx_type n = r.nelem ();
940 947 ComplexMatrix result (1, n);
941 pmat[0] = std::pow (a, r.base ()); 948
942 Complex mul = std::pow (a, r.inc ()); 949 if (same_sign (r.base (), r.inc ()))
943 for (octave_idx_type i = 1; i < r.nelem (); i++) 950 {
944 pmat[i] = pmat[i-1] * mul; 951 Complex base = std::pow (a, r.base ());
945 952 Complex inc = std::pow (a, r.inc ());
946 retval = mat; 953 result(0) = base;
947 } 954 for (octave_idx_type i = 1; i < n; i++)
955 result(i) = (base *= inc);
956 }
957 else
958 {
959 // Don't use Range::limit () here.
960 Complex limit = std::pow (a, r.base () + (n-1) * r.inc ());
961 Complex inc = std::pow (a, -r.inc ());
962 result(n-1) = limit;
963 for (octave_idx_type i = n-2; i >= 0; i--)
964 result(i) = (limit *= inc);
965 }
966
967 retval = result;
968 }
969 else
970 retval = elem_xpow (a, r.matrix_value ());
971
948 972
949 return retval; 973 return retval;
950 } 974 }
951 975
952 // -*- 9 -*- 976 // -*- 9 -*-