Skip to content

Commit

Permalink
Merge pull request #1620 from GillesDuvert/integerPowerOptimizations
Browse files Browse the repository at this point in the history 
Integer power optimizations
  • Loading branch information
@GillesDuvert
GillesDuvert authored Aug 27, 2023
2 parents 052747a + 1ef252f commit 7bf88b8
Show file tree
Hide file tree
Showing 8 changed files with 333 additions and 467 deletions.
118 changes: 16 additions & 102 deletions src/basic_op.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,9 @@
#include "typetraits.hpp"

#include "sigfpehandler.hpp"

#include "gdl_util.hpp" //for gdl::powI

using namespace std;

#if defined(USE_EIGEN)
Expand Down Expand Up @@ -3131,34 +3134,6 @@ Data_<SpDObj>* Data_<SpDObj>::ModInvS(BaseGDL* r) {
return this;
}

// Pow
// C++ defines pow only for floats and doubles
//template <typename T, typename TT> T pow( const T r, const TT l)

template <typename T> 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<class Sp>
Data_<Sp>* Data_<Sp>::Pow(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
Data_* right = static_cast<Data_*> (r);
Expand Down Expand Up @@ -3216,17 +3191,9 @@ Data_<SpDFloat>* Data_<SpDFloat>::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<class Sp>
Data_<Sp>* Data_<Sp>::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
assert(0);
return this;
}
// floats power of value with GDL_LONG: left=left ^ right

template<>
Data_<SpDFloat>* Data_<SpDFloat>::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
DLongGDL* right = static_cast<DLongGDL*> (r);

ULong rEl = right->N_Elements();
Expand All @@ -3237,11 +3204,11 @@ Data_<SpDFloat>* Data_<SpDFloat>::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;
}
Expand All @@ -3250,97 +3217,44 @@ Data_<SpDFloat>* Data_<SpDFloat>::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_<SpDDouble>* Data_<SpDDouble>::PowInt(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
DLongGDL* right = static_cast<DLongGDL*> (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_<SpDString>* Data_<SpDString>::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_<SpDFloat>* Data_<SpDFloat>::PowInv(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
Data_* right = static_cast<Data_*> (r);
Expand Down
Loading

0 comments on commit 7bf88b8

Please sign in to comment.