Function nc_pow

numexpr3/complex_functions.hpp:833–877

RAM: algorithm branches, and is really not designed for a modern CPU... TODO: we should be able to branch the exp>-100 && exp<100 test before the for loop.

// TODO: we should be able to branch the exp>-100 && exp<100 test before the 
// for loop.
static void
nc_pow( npy_intp n, npy_complex64 *a, npy_intp sb1, npy_complex64 *b, npy_intp sb2, npy_complex64 *r)
{
    npy_intp exp, mask;
    npy_complex64 p, aa;
    sb1 /= sizeof(npy_complex64);
    sb2 /= sizeof(npy_complex64);
                 
    for( int I = 0; I < n; I++ ) {
        if (b[I*sb2].real == 0. && b[I*sb2].imag == 0.) {
            r[I].real = 1.;
            r[I].imag = 0.;
            continue;
        }
        if (a[I*sb1].real == 0. && a[I*sb1].imag == 0.) {
            r[I].real = 0.;
            r[I].imag = 0.;
            continue;
        }
        if (b[I*sb2].imag == 0 && (exp=(npy_intp)b[I*sb2].real) == b[I*sb2].real) {
            if (exp > -100 && exp < 100) {
                mask = 1;
                exp = (exp < 0) ? -exp : exp;
                      
                aa = Z_1;
                p.real = a[I*sb1].real; p.imag = a[I*sb1].imag;
                while (1) {
                    aa = (exp & mask) ? _return_mul( aa, p ) : aa;
                         
                    mask <<= 1;
                    if (n < mask || mask <= 0) break;
                    p = _return_mul( p, p );
                }
                r[I].real = aa.real; 
                r[I].imag = aa.imag;
                r[I] = (b[I*sb2].real < 0) ? _return_div( Z_1, r[I] ) : r[I];
                continue;
            }
        }
        
        _inline_log( a[I*sb1], r[I] );
        _inline_mul(r[I], b[I*sb2], r[I] );
        _inline_exp(r[I], r[I] );
    }
}
        
static void
nc_pow( npy_intp n, npy_complex128 *a, npy_intp sb1, npy_complex128 *b, npy_intp sb2, npy_complex128 *r)