diff --git a/src/basic_op.cpp b/src/basic_op.cpp index 08e18a863..588208653 100644 --- a/src/basic_op.cpp +++ b/src/basic_op.cpp @@ -32,6 +32,9 @@ #include "typetraits.hpp" #include "sigfpehandler.hpp" + +#include "gdl_util.hpp" //for gdl::powI + using namespace std; #if defined(USE_EIGEN) @@ -3131,34 +3134,6 @@ Data_* Data_::ModInvS(BaseGDL* r) { return this; } -// Pow -// C++ defines pow only for floats and doubles -//template T pow( const T r, const TT l) - -template T pow(const T r, const T l) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) - typedef T TT; - - if (l == 0) return 1; - if (l < 0) return 0; - - const int nBits = sizeof (TT) * 8; - - T arr = r; - T res = 1; - TT mask = 1; - for (SizeT i = 0; i < nBits; ++i) { - if (l & mask) res *= arr; - mask <<= 1; - if (l < mask) return res; - arr *= arr; - } - - return res; -} - -// power of value: left=left ^ right -// integral types - template Data_* Data_::Pow(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) Data_* right = static_cast (r); @@ -3216,17 +3191,9 @@ Data_* Data_::Pow(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,_ return this; } -// PowInt and PowIntNew can only be called for GDL_FLOAT and GDL_DOUBLE - +// anygdl (except complex) power to a GDL_LONG value left=left ^ right template Data_* Data_::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) - assert(0); - return this; -} -// floats power of value with GDL_LONG: left=left ^ right - -template<> -Data_* Data_::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) DLongGDL* right = static_cast (r); ULong rEl = right->N_Elements(); @@ -3237,11 +3204,11 @@ Data_* Data_::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION_ DLong r0 = (*right)[0]; if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], r0); + for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = gdl::powI((*this)[i], r0); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], r0); + for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = gdl::powI((*this)[i], r0); } return this; } @@ -3250,97 +3217,44 @@ Data_* Data_::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION_ Ty s0 = (*this)[ 0]; if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = gdl::powI(s0, (*right)[ i]); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = gdl::powI(s0, (*right)[ i]); } return res; } if (nEl <= rEl) { if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = gdl::powI((*this)[i], (*right)[i]); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = gdl::powI((*this)[i], (*right)[i]); } return this; } else { Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = gdl::powI((*this)[i], (*right)[i]); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = gdl::powI((*this)[i], (*right)[i]); } return res; } } - template<> -Data_* Data_::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) - DLongGDL* right = static_cast (r); - - ULong rEl = right->N_Elements(); - ULong nEl = N_Elements(); - assert(rEl); - assert(nEl); - if (r->StrictScalar()) { - DLong r0 = (*right)[0]; - - if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], r0); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], r0); - } - return this; - } - if (StrictScalar()) { - Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); - Ty s0 = (*this)[ 0]; - - if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); - } - return res; - } - if (nEl <= rEl) { - - if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], (*right)[i]); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = pow((*this)[i], (*right)[i]); - } - return this; - } else { - Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); - - if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); - } - return res; - } +Data_* Data_::PowInt(BaseGDL* r) { + assert(0); + throw GDLException("Internal error: Data_::PowInt called.", true, false); + return NULL; } -// floats inverse power of value: left=right ^ left - template<> Data_* Data_::PowInv(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) Data_* right = static_cast (r); diff --git a/src/basic_op_new.cpp b/src/basic_op_new.cpp index fc729ab3f..a1b6d5acf 100644 --- a/src/basic_op_new.cpp +++ b/src/basic_op_new.cpp @@ -33,6 +33,8 @@ #include "sigfpehandler.hpp" +#include "gdl_util.hpp" //for gdl::powI + using namespace std; // binary operators @@ -2301,36 +2303,6 @@ Data_* Data_::ModInvSNew(BaseGDL* r) { } - -// Pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow -// C++ defines pow only for floats and doubles -// in basic_op.cpp: -// template T pow( const T r, const T l) -// { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) -// typedef T TT; -// -// if( l == 0) return 1; -// if( l < 0) return 0; -// -// const int nBits = sizeof(TT) * 8; -// -// T arr = r; -// T res = 1; -// TT mask = 1; -// for( SizeT i=0; i Data_* Data_::PowNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) Data_* right = static_cast (r); @@ -2352,7 +2324,7 @@ Data_* Data_::PowNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,_ return res; } // inverse power of value: left=right ^ left - +//PowNew: operands A and B are both arrays. n_el(A) > n_el(B), so this is B and right is A (inverted mode) template Data_* Data_::PowInvNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) Data_* right = static_cast (r); @@ -2374,18 +2346,9 @@ Data_* Data_::PowInvNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE_ return res; } -// PowInt and PowIntNew can only be called for GDL_FLOAT and GDL_DOUBLE - +// anygdl (except complex) power to a GDL_LONG value left=left ^ right template -Data_* Data_::PowIntNew(BaseGDL* r) { - assert(0); - throw GDLException("Internal error: Data_::PowIntNew called.", true, false); - return NULL; -} -// floats power of value with GDL_LONG: left=left ^ right - -template<> -Data_* Data_::PowIntNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) +Data_* Data_::PowIntNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) DLongGDL* right = static_cast (r); ULong rEl = right->N_Elements(); @@ -2393,14 +2356,14 @@ Data_* Data_::PowIntNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTI assert(rEl); assert(nEl); if (r->StrictScalar()) { - Data_* res = new Data_(Dim(), BaseGDL::NOZERO); + Data_* res = new Data_(this->Dim(), BaseGDL::NOZERO); DLong r0 = (*right)[0]; if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], r0); + for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = gdl::powI((*this)[i], r0); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], r0); + for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = gdl::powI((*this)[i], r0); } return res; } @@ -2408,94 +2371,43 @@ Data_* Data_::PowIntNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTI Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); Ty s0 = (*this)[ 0]; if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = gdl::powI(s0, (*right)[ i]); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = gdl::powI(s0, (*right)[ i]); } return res; } if (nEl <= rEl) { - Data_* res = new Data_(Dim(), BaseGDL::NOZERO); + Data_* res = new Data_(this->Dim(), BaseGDL::NOZERO); if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = gdl::powI((*this)[i], (*right)[i]); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = gdl::powI((*this)[i], (*right)[i]); } return res; } else { Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = gdl::powI((*this)[i], (*right)[i]); } else { TRACEOMP(__FILE__, __LINE__) #pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); + for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = gdl::powI((*this)[i], (*right)[i]); } return res; } } - template<> -Data_* Data_::PowIntNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) - DLongGDL* right = static_cast (r); - - ULong rEl = right->N_Elements(); - ULong nEl = N_Elements(); - assert(rEl); - assert(nEl); - if (r->StrictScalar()) { - Data_* res = new Data_(Dim(), BaseGDL::NOZERO); - DLong r0 = (*right)[0]; - if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], r0); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], r0); - } - return res; - } - if (StrictScalar()) { - Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); - Ty s0 = (*this)[ 0]; - if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[ i] = pow(s0, (*right)[ i]); - } - return res; - } - if (nEl <= rEl) { - Data_* res = new Data_(Dim(), BaseGDL::NOZERO); - if ((GDL_NTHREADS=parallelize( nEl))==1) { - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); - } - return res; - } else { - Data_* res = new Data_(right->Dim(), BaseGDL::NOZERO); - if ((GDL_NTHREADS=parallelize( rEl))==1) { - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); - } else { - TRACEOMP(__FILE__, __LINE__) -#pragma omp parallel for num_threads(GDL_NTHREADS) - for (OMPInt i = 0; i < rEl; ++i) (*res)[i] = pow((*this)[i], (*right)[i]); - } - return res; - } +Data_* Data_::PowIntNew(BaseGDL* r) { + assert(0); + throw GDLException("Internal error: Data_::PowIntNew called.", true, false); + return NULL; } -// floats power of value: left=left ^ right - template<> Data_* Data_::PowNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) Data_* right = static_cast (r); @@ -2977,6 +2889,7 @@ Data_* Data_::PowInvNew(BaseGDL* r) { } // scalar versions +//PowInvSNew: A^B: operand B is singleton and A array template Data_* Data_::PowSNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) @@ -3001,6 +2914,7 @@ Data_* Data_::PowSNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__, return res; } // inverse power of value: left=right ^ left +//PowInvSNew: A^B: operand A is singleton and B array template Data_* Data_::PowInvSNew(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__) diff --git a/src/gdl_util.hpp b/src/gdl_util.hpp new file mode 100644 index 000000000..2f96ef164 --- /dev/null +++ b/src/gdl_util.hpp @@ -0,0 +1,93 @@ +#ifndef GDL_UTIL_HPP_ +#define GDL_UTIL_HPP_ +#include "basegdl.hpp" + +//templates must be in the same namespace as where thet are used. +namespace gdl { +//in practice due to specialization this is for integer types only. + template T1 powI(const T1 x, const DLong y) { +//All x integer types + if (y == 0) return 1; + if (y < 0) { + return x==1?1:0; //do not forget 1 + } + + const int nBits = sizeof (DLong) * 8; + + T1 arr = x; + T1 res = 1; + DLong mask = 1; + for (SizeT i = 0; i < nBits; ++i) { + if (y & mask) res *= arr; + mask <<= 1; + if (y < mask) return res; + arr *= arr; + } + + return res; +} + + //make specializations inline as they are seen by more than one .obj +template <> +inline DFloat powI(const DFloat x, const DLong yy) { +//All x integer types + if (yy == 0) return 1; + const int nBits = sizeof (DLong) * 8; + + + + DLong y = yy; + DFloat arr = x; + if (y < 0) { + arr=1/arr; + y=-y; + } + DFloat res = 1; + DLong mask = 1; + for (SizeT i = 0; i < nBits; ++i) { + if (y & mask) res *= arr; + mask <<= 1; + if (y < mask) return res; + arr *= arr; + } + + return res; +} +template <> +inline DDouble powI(const DDouble x, const DLong yy) { +//All x integer types + if (yy == 0) return 1; + const int nBits = sizeof (DLong) * 8; + + + + DLong y = yy; + DDouble arr = x; + if (y < 0) { + arr=1/arr; + y=-y; + } + DDouble res = 1; + DLong mask = 1; + for (SizeT i = 0; i < nBits; ++i) { + if (y & mask) res *= arr; + mask <<= 1; + if (y < mask) return res; + arr *= arr; + } + + return res; +} +template <> +inline DComplex powI(const DComplex x, const DLong y) { + assert(0); + throw GDLException("powI() not defined for DComplex"); +} +template <> +inline DComplexDbl powI(const DComplexDbl x, const DLong y) { + assert(0); + throw GDLException("powI() not defined for DComplexDbl"); +} +} + +#endif diff --git a/src/gdlwxstream.cpp b/src/gdlwxstream.cpp index 2149e3b3a..e9ea2c494 100644 --- a/src/gdlwxstream.cpp +++ b/src/gdlwxstream.cpp @@ -569,15 +569,19 @@ bool GDLWXStream::GetGin(PLGraphicsIn *gin, int mode) { DOWN, //3 UP //4 }; + Update(); wxMouseState mouse=wxGetMouseState(); wxPoint mousePoint = wxGetMousePosition (); unsigned int state=0, ostate=0; //start with null state unsigned int button=0; //gets the button state 1,2,3 int x,y; - x = mouse.GetX(); - y = mouse.GetY(); - gin->pX = x; - gin->pY = y; + x = 0; + y = 0; + if (container->GetScreenRect().Contains(wxGetMousePosition())) { //if cursor is in the window... + mouse=wxGetMouseState(); + x = mouse.GetX(); + y = mouse.GetY(); + } //state is like the button state value, combination of all buttonstate masks of X11. It is merely to know what buttons are down if (mouse.LeftIsDown()) {button = 1; ostate |= 1;} @@ -615,7 +619,7 @@ bool GDLWXStream::GetGin(PLGraphicsIn *gin, int mode) { if (mode==UP && ButtonRelease) goto end; if (!ButtonRelease && (mode==WAIT || mode==DOWN) ) goto end; } - if ( (pow((float)x-gin->pX,2)+pow((float)y-gin->pY,2) > 0) && mode==CHANGE) { + if ( ( (x-gin->pX)*(x-gin->pX)+(y-gin->pY)*(y-gin->pY) > 0 ) && mode==CHANGE) { gin->pX = x; gin->pY = y; goto end; diff --git a/src/math_fun_jmg.cpp b/src/math_fun_jmg.cpp index 3d9f6cfd5..18f83b5c8 100644 --- a/src/math_fun_jmg.cpp +++ b/src/math_fun_jmg.cpp @@ -28,6 +28,7 @@ #include "math_utl.hpp" #include "math_fun_jmg.hpp" #include "graphicsdevice.hpp" +#include "gdl_util.hpp" //#define GDL_DEBUG #undef GDL_DEBUG @@ -2042,52 +2043,6 @@ namespace lib { } -/*-------------------------------------------------------------------------*/ -/** - @brief Same as pow(x,y) but for integer values of y only (faster). - @param x A double number. - @param p An integer power. - @return x to the power p. - - This is much faster than the math function due to the integer. Some - compilers make this optimization already, some do not. - - p can be positive, negative or null. - */ -/*--------------------------------------------------------------------------*/ -double ipow(double x, int p) -{ - double r, recip ; - - /* Get rid of trivial cases */ - switch (p) { - case 0: - return 1.00 ; - - case 1: - return x ; - - case 2: - return x*x ; - - case 3: - return x*x*x ; - - case -1: - return 1.00 / x ; - - case -2: - return (1.00 / x) * (1.00 / x) ; - } - if (p>0) { - r = x ; - while (--p) r *= x ; - } else { - r = recip = 1.00 / x ; - while (++p) r *= recip ; - } - return r; -} /*-------------------------------------------------------------------------*/ /** @@ -2111,7 +2066,7 @@ double poly2d_compute( z = 0.00 ; for (i=0 ; inc ; i++) { - z += p->c[i] * ipow(x, p->px[i]) * ipow(y, p->py[i]) ; + z += p->c[i] * gdl::powI(x, p->px[i]) * gdl::powI(y, p->py[i]) ; } return z ; } @@ -2173,10 +2128,10 @@ DDouble * generate_interpolation_kernel(int kernel_type, DDouble cubicParameter) for (i=1 ; iKeywordSet(UPIx)) wait=UP; if(actStream->GetGin(&gin, wait)==false) return; // outside window report -1 -1 at least for DEVICE values - if (gin.pX < 0 || gin.pX > actStream->xPageSize() || gin.pY < 0 || gin.pY > actStream->yPageSize()) - { - gin.pX = -1; - gin.pY = -1; - } + bool out=(gin.pX < 0 || gin.pX > actStream->xPageSize() || gin.pY < 0 || gin.pY > actStream->yPageSize()); if (e->KeywordSet(DEVICEIx)) { + if (out) { gin.pX = -1; gin.pY = -1;} DLongGDL* xLong; DLongGDL* yLong; xLong = new DLongGDL(gin.pX); @@ -219,6 +216,7 @@ void cursor(EnvT* e){ } else { + if (out) { gin.dX = 0; gin.dY = 0;} DDoubleGDL* x; DDoubleGDL* y; if (e->KeywordSet(NORMALIx)) @@ -228,6 +226,7 @@ void cursor(EnvT* e){ } else { // default (/data) + if (out) { gin.dX = 0; gin.dY = 0;} DDouble tempx,tempy; #ifdef USE_LIBPROJ bool mapSet = false; diff --git a/src/prognodeexpr.cpp b/src/prognodeexpr.cpp index 3025807c0..168e7932c 100644 --- a/src/prognodeexpr.cpp +++ b/src/prognodeexpr.cpp @@ -1268,114 +1268,115 @@ BaseGDL* MOD_OPNode::Eval() BaseGDL* POWNode::Eval() { BaseGDL* res; - Guard e1( op1->Eval()); - Guard e2( op2->Eval()); + Guard e1(op1->Eval()); + Guard e2(op2->Eval()); // special handling for aTy == complex && bTy != complex - DType aTy=e1->Type(); - DType bTy=e2->Type(); + DType aTy = e1->Type(); + DType bTy = e2->Type(); if( aTy == GDL_STRING) { - e1.reset( e1->Convert2( GDL_FLOAT, BaseGDL::COPY)); - aTy = GDL_FLOAT; - } + e1.reset(e1->Convert2(GDL_FLOAT, BaseGDL::COPY)); + aTy = GDL_FLOAT; + } if( bTy == GDL_STRING) { - e2.reset( e2->Convert2( GDL_FLOAT, BaseGDL::COPY)); - bTy = GDL_FLOAT; - } + e2.reset(e2->Convert2(GDL_FLOAT, BaseGDL::COPY)); + bTy = GDL_FLOAT; + } + //powers of complex must use std::pow always, no "integer power" refinement possible. if( ComplexType( aTy)) { if( IntType( bTy)) { - e2.reset( e2.release()->Convert2( GDL_LONG)); - res = e1->Pow( e2.get()); - if( res == e1.get()) - e1.release(); - return res; - } + e2.reset(e2.release()->Convert2(GDL_LONG)); + res = e1->Pow(e2.get()); + if (res == e1.get()) + e1.release(); + return res; + } if( aTy == GDL_COMPLEX) { if( bTy == GDL_DOUBLE) { - e1.reset( e1.release()->Convert2( GDL_COMPLEXDBL)); - aTy = GDL_COMPLEXDBL; + e1.reset(e1.release()->Convert2(GDL_COMPLEXDBL)); + aTy = GDL_COMPLEXDBL; } else if( bTy == GDL_FLOAT) { - res = e1->Pow( e2.get()); - if( res == e1.get()) - e1.release(); - return res; - } - } + res = e1->Pow(e2.get()); + if (res == e1.get()) + e1.release(); + return res; + } + } if( aTy == GDL_COMPLEXDBL) { if( bTy == GDL_FLOAT) { - e2.reset( e2.release()->Convert2( GDL_DOUBLE)); - bTy = GDL_DOUBLE; - } + e2.reset(e2.release()->Convert2(GDL_DOUBLE)); + bTy = GDL_DOUBLE; + } if( bTy == GDL_DOUBLE) { - res = e1->Pow( e2.get()); - if( res == e1.get()) - e1.release(); - return res; - } - } + res = e1->Pow(e2.get()); + if (res == e1.get()) + e1.release(); + return res; + } } + } - if( IntType( bTy) && FloatType( aTy)) - { - e2.reset( e2.release()->Convert2( GDL_LONG)); + // GD: simplify: force use of PowInt for all integer powers, not only integer powers of float types as previously + if (IntType(bTy)) + { + e2.reset(e2.release()->Convert2(GDL_LONG)); - res = e1->PowInt( e2.get()); - if( res == e1.get()) - e1.release(); + res = e1->PowInt(e2.get()); + if (res == e1.get()) + e1.release(); - return res; - } - - DType convertBackT; + return res; + } + + DType convertBackT; // convert back - if( IntType( bTy) && (DTypeOrder[ bTy] > DTypeOrder[ aTy])) -// if( IntType( bTy) && (DTypeOrder[ bTy] > DTypeOrder[ aTy])) + if (IntType(bTy) && (DTypeOrder[ bTy] > DTypeOrder[ aTy])) + // if( IntType( bTy) && (DTypeOrder[ bTy] > DTypeOrder[ aTy])) convertBackT = aTy; else convertBackT = GDL_UNDEF; - + // first operand determines type JMG // AdjustTypes(e2,e1); // order crucial here (for converting back) - AdjustTypes(e1,e2); // order crucial here (for converting back) + AdjustTypes(e1, e2); // order crucial here (for converting back) if( e1->StrictScalar()) { - res= e2->PowInvS(e1.get()); // scalar+scalar or array+scalar - e2.release(); + res = e2->PowInvS(e1.get()); // scalar+scalar or array+scalar + e2.release(); } else if( e2->StrictScalar()) { - res= e1->PowS(e2.get()); // array+scalar - e1.release(); + res = e1->PowS(e2.get()); // array+scalar + e1.release(); } else if( e1->N_Elements() <= e2->N_Elements()) { - res= e1->Pow(e2.get()); // smaller_array + larger_array or same size - e1.release(); + res = e1->Pow(e2.get()); // smaller_array + larger_array or same size + e1.release(); } else { - res= e2->PowInv(e1.get()); // smaller + larger - e2.release(); - } + res = e2->PowInv(e1.get()); // smaller + larger + e2.release(); + } if( convertBackT != GDL_UNDEF) { - res = res->Convert2( convertBackT, BaseGDL::CONVERT); - } - endPOW: + res = res->Convert2(convertBackT, BaseGDL::CONVERT); + } return res; } // BaseGDL* DECNode::Eval() @@ -2883,231 +2884,229 @@ BaseGDL* POWNCNode::Eval() BaseGDL *e1, *e2; if( op1NC) { - e1 = op1->EvalNC(); + e1 = op1->EvalNC(); } else { - e1 = op1->Eval(); - g1.reset( e1); - } + e1 = op1->Eval(); + g1.reset(e1); + } if( op2NC) { - e2 = op2->EvalNC(); + e2 = op2->EvalNC(); } else { - e2 = op2->Eval(); - g2.reset( e2); - } + e2 = op2->Eval(); + g2.reset(e2); + } - DType aTy=e1->Type(); - DType bTy=e2->Type(); + DType aTy = e1->Type(); + DType bTy = e2->Type(); if( aTy == GDL_STRING) { - e1 = e1->Convert2( GDL_FLOAT, BaseGDL::COPY); - g1.reset( e1); // deletes old e1 - aTy = GDL_FLOAT; - } + e1 = e1->Convert2(GDL_FLOAT, BaseGDL::COPY); + g1.reset(e1); // deletes old e1 + aTy = GDL_FLOAT; + } if( bTy == GDL_STRING) { - e2 = e2->Convert2( GDL_FLOAT, BaseGDL::COPY); - g2.reset( e2); // deletes old e2 - bTy = GDL_FLOAT; - } - + e2 = e2->Convert2(GDL_FLOAT, BaseGDL::COPY); + g2.reset(e2); // deletes old e2 + bTy = GDL_FLOAT; + } + //powers of complex must use std::pow always, no "integer power" refinement possible. if( ComplexType(aTy)) { if( IntType( bTy)) { if( bTy != GDL_LONG) { - e2 = e2->Convert2( GDL_LONG, BaseGDL::COPY); - g2.reset( e2); - } + e2 = e2->Convert2(GDL_LONG, BaseGDL::COPY); + g2.reset(e2); + } if( g1.get() == NULL) { - return e1->PowNew( e2); + return e1->PowNew(e2); } else { - res = g1->Pow( e2); - if( res == g1.get()) - g1.release(); - return res; - } - } + res = g1->Pow(e2); + if (res == g1.get()) + g1.release(); + return res; + } + } if( aTy == GDL_COMPLEX) { if( bTy == GDL_DOUBLE) { - e1 = e1->Convert2( GDL_COMPLEXDBL, BaseGDL::COPY); - g1.reset( e1); - aTy = GDL_COMPLEXDBL; + e1 = e1->Convert2(GDL_COMPLEXDBL, BaseGDL::COPY); + g1.reset(e1); + aTy = GDL_COMPLEXDBL; } else if( bTy == GDL_FLOAT) { if( g1.get() == NULL) { - return e1->PowNew( e2); + return e1->PowNew(e2); } else { - res = g1->Pow( e2); - if( res == g1.get()) - g1.release(); - return res; - } - } - } + res = g1->Pow(e2); + if (res == g1.get()) + g1.release(); + return res; + } + } + } if( aTy == GDL_COMPLEXDBL) { if( bTy == GDL_FLOAT) { - e2 = e2->Convert2( GDL_DOUBLE, BaseGDL::COPY); - g2.reset( e2); - bTy = GDL_DOUBLE; - } + e2 = e2->Convert2(GDL_DOUBLE, BaseGDL::COPY); + g2.reset(e2); + bTy = GDL_DOUBLE; + } if( bTy == GDL_DOUBLE) { if( g1.get() == NULL) { - return e1->PowNew( e2); + return e1->PowNew(e2); } else { - res = g1->Pow( e2); - if( res == g1.get()) - g1.release(); - return res; - } - } - } + res = g1->Pow(e2); + if (res == g1.get()) + g1.release(); + return res; + } + } } + } - if( IntType( bTy) && FloatType( aTy)) - { - if( bTy != GDL_LONG) - { - e2 = e2->Convert2( GDL_LONG, BaseGDL::COPY); - g2.reset( e2); - } - if( g1.get() == NULL) - res = e1->PowIntNew( e2); - else - { - res = g1->PowInt( e2); - if( res == g1.get()) - g1.release(); - } - return res; + // GD: simplify: force use of PowInt for all integer powers, not only integer powers of float types as previously + if (IntType(bTy)) { + if (bTy != GDL_LONG) { + e2 = e2->Convert2(GDL_LONG, BaseGDL::COPY); + g2.reset(e2); } + if (g1.get() == NULL) + res = e1->PowIntNew(e2); + else { + res = g1->PowInt(e2); + if (res == g1.get()) + g1.release(); + } + return res; + } + + DType convertBackT; - DType convertBackT; - bool aTyGEbTy = DTypeOrder[aTy] >= DTypeOrder[bTy]; // convert back - if( IntType( bTy) && !aTyGEbTy) + if (IntType(bTy) && !aTyGEbTy) convertBackT = aTy; else convertBackT = GDL_UNDEF; if( aTy != bTy) { - if( aTyGEbTy) // crucial: '>' -> '>=' - { + if (aTyGEbTy) // crucial: '>' -> '>=' + { if( DTypeOrder[aTy] > 100) { - throw GDLException( "Expressions of this type cannot be converted."); - } + throw GDLException("Expressions of this type cannot be converted."); + } - // convert e2 to e1 - e2 = e2->Convert2( aTy, BaseGDL::COPY); - g2.reset( e2); // delete former e2 + // convert e2 to e1 + e2 = e2->Convert2(aTy, BaseGDL::COPY); + g2.reset(e2); // delete former e2 } else // bTy > aTy (order) - { + { if( DTypeOrder[bTy] > 100) { - throw GDLException( "Expressions of this type cannot be converted."); - } + throw GDLException("Expressions of this type cannot be converted."); + } - // convert e1 to e2 - e1 = e1->Convert2( bTy, BaseGDL::COPY); - g1.reset( e1); // delete former e1 - } + // convert e1 to e2 + e1 = e1->Convert2(bTy, BaseGDL::COPY); + g1.reset(e1); // delete former e1 } + } // AdjustTypes(e2,e1); // order crucial here (for converting back) if( e1->StrictScalar()) { - if( g2.get() == NULL) - res = e2->PowInvSNew( e1); + if (g2.get() == NULL) + res = e2->PowInvSNew(e1); else { - g2.release(); - res = e2->PowInvS(e1); // scalar+scalar or array+scalar - } + g2.release(); + res = e2->PowInvS(e1); // scalar+scalar or array+scalar + } } else if( e2->StrictScalar()) { if( g1.get() == NULL) { - res = e1->PowSNew( e2); -// res = e1->PowS(e2); // array+scalar + res = e1->PowSNew(e2); + // res = e1->PowS(e2); // array+scalar } else { - g1.release(); - res= e1->PowS(e2); // array+scalar - } + g1.release(); + res = e1->PowS(e2); // array+scalar + } } else if( e1->N_Elements() == e2->N_Elements()) { if( g1.get() != NULL) { - g1.release(); - res = e1->Pow(e2); + g1.release(); + res = e1->Pow(e2); } else if( g2.get() != NULL) { - g2.release(); - res = e2->PowInv(e1); - res->SetDim( e1->Dim()); + g2.release(); + res = e2->PowInv(e1); + res->SetDim(e1->Dim()); } else { - res = e1->PowNew( e2); -// res = e1->Pow(e2); - } + res = e1->PowNew(e2); + // res = e1->Pow(e2); + } } else if( e1->N_Elements() < e2->N_Elements()) { if( g1.get() == NULL) { - res = e1->PowNew( e2); -// res= e1->Pow(e2); // smaller_array + larger_array or same size + res = e1->PowNew(e2); + // res= e1->Pow(e2); // smaller_array + larger_array or same size } else { - g1.release(); - res= e1->Pow(e2); // smaller_array + larger_array or same size - } + g1.release(); + res = e1->Pow(e2); // smaller_array + larger_array or same size + } } else { - if( g2.get() == NULL) - res = e2->PowInvNew( e1); + if (g2.get() == NULL) + res = e2->PowInvNew(e1); else { - g2.release(); - res= e2->PowInv(e1); // smaller + larger - } - } + g2.release(); + res = e2->PowInv(e1); // smaller + larger + } + } if( convertBackT != GDL_UNDEF) { - res = res->Convert2( convertBackT, BaseGDL::CONVERT); - } + res = res->Convert2(convertBackT, BaseGDL::CONVERT); + } return res; } diff --git a/src/specializations.hpp b/src/specializations.hpp index 165b794b3..686865ce1 100644 --- a/src/specializations.hpp +++ b/src/specializations.hpp @@ -317,18 +317,6 @@ Data_* Data_::Pow( BaseGDL* r); template<> Data_* Data_::PowInv( BaseGDL* r); template<> -Data_* Data_::PowInt( BaseGDL* r); -template<> -Data_* Data_::PowIntNew( BaseGDL* r); -template<> -Data_* Data_::PowInt( BaseGDL* r); -template<> -Data_* Data_::PowIntNew( BaseGDL* r); -template<> -Data_* Data_::PowInt( BaseGDL* r); -template<> -Data_* Data_::PowIntNew( BaseGDL* r); -template<> Data_* Data_::Pow( BaseGDL* r); template<> Data_* Data_::PowInv( BaseGDL* r);