/****************************  vectormath_exp.h   ******************************
* Author:        Agner Fog
* Date created:  2014-04-18
* Last modified: 2014-10-22
* Version:       1.16
* Project:       vector classes
* Description:
* Header file containing inline vector functions of logarithms, exponential 
* and power functions:
* exp         exponential function
* exmp1       exponential function minus 1
* log         natural logarithm
* log1p       natural logarithm of 1+x
* cbrt        cube root
* pow         raise vector elements to power
* pow_ratio   raise vector elements to rational power
*
* Theory, methods and inspiration based partially on these sources:
* > Moshier, Stephen Lloyd Baluk: Methods and programs for mathematical functions.
*   Ellis Horwood, 1989.
* > VDT library developed on CERN by Danilo Piparo, Thomas Hauth and
*   Vincenzo Innocente, 2012, https://svnweb.cern.ch/trac/vdt
* > Cephes math library by Stephen L. Moshier 1992,
*   http://www.netlib.org/cephes/
*
* For detailed instructions, see vectormath_common.h and VectorClass.pdf
*
* (c) Copyright 2014 GNU General Public License http://www.gnu.org/licenses
******************************************************************************/

#ifndef VECTORMATH_EXP_H
#define VECTORMATH_EXP_H  1 

#include "vectormath_common.h"  


/******************************************************************************
*                 Exponential functions
******************************************************************************/

// Helper functions, used internally:

// This function calculates pow(2,n) where n must be an integer. Does not check for overflow or underflow
static inline Vec2d vm_pow2n (Vec2d const & n) {
    const double pow2_52 = 4503599627370496.0;   // 2^52
    const double bias = 1023.0;                  // bias in exponent
    Vec2d a = n + (bias + pow2_52);              // put n + bias in least significant bits
    Vec2q b = reinterpret_i(a);                  // bit-cast to integer
    Vec2q c = b << 52;                           // shift left 52 places to get into exponent field
    Vec2d d = reinterpret_d(c);                  // bit-cast back to double
    return d;
}

static inline Vec4f vm_pow2n (Vec4f const & n) {
    const float pow2_23 =  8388608.0;            // 2^23
    const float bias = 127.0;                    // bias in exponent
    Vec4f a = n + (bias + pow2_23);              // put n + bias in least significant bits
    Vec4i b = reinterpret_i(a);                  // bit-cast to integer
    Vec4i c = b << 23;                           // shift left 23 places to get into exponent field
    Vec4f d = reinterpret_f(c);                  // bit-cast back to float
    return d;
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec4d vm_pow2n (Vec4d const & n) {
    const double pow2_52 = 4503599627370496.0;   // 2^52
    const double bias = 1023.0;                  // bias in exponent
    Vec4d a = n + (bias + pow2_52);              // put n + bias in least significant bits
    Vec4q b = reinterpret_i(a);                  // bit-cast to integer
    Vec4q c = b << 52;                           // shift left 52 places to get value into exponent field
    Vec4d d = reinterpret_d(c);                  // bit-cast back to double
    return d;
}

static inline Vec8f vm_pow2n (Vec8f const & n) {
    const float pow2_23 =  8388608.0;            // 2^23
    const float bias = 127.0;                    // bias in exponent
    Vec8f a = n + (bias + pow2_23);              // put n + bias in least significant bits
    Vec8i b = reinterpret_i(a);                  // bit-cast to integer
    Vec8i c = b << 23;                           // shift left 23 places to get into exponent field
    Vec8f d = reinterpret_f(c);                  // bit-cast back to float
    return d;
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec8d vm_pow2n (Vec8d const & n) {
    const double pow2_52 = 4503599627370496.0;   // 2^52
    const double bias = 1023.0;                  // bias in exponent
    Vec8d a = n + (bias + pow2_52);              // put n + bias in least significant bits
    Vec8q b = Vec8q(reinterpret_i(a));           // bit-cast to integer
    Vec8q c = b << 52;                           // shift left 52 places to get value into exponent field
    Vec8d d = Vec8d(reinterpret_d(c));           // bit-cast back to double
    return d;
}

static inline Vec16f vm_pow2n (Vec16f const & n) {
    const float pow2_23 =  8388608.0;            // 2^23
    const float bias = 127.0;                    // bias in exponent
    Vec16f a = n + (bias + pow2_23);             // put n + bias in least significant bits
    Vec16i b = Vec16i(reinterpret_i(a));         // bit-cast to integer
    Vec16i c = b << 23;                          // shift left 23 places to get into exponent field
    Vec16f d = Vec16f(reinterpret_f(c));         // bit-cast back to float
    return d;
}

#endif // MAX_VECTOR_SIZE >= 512


// Template for exp function, double precision
// The limit of abs(x) is defined by max_x below
// This function does not produce denormals
// Template parameters:
// VTYPE: double vector type
// BTYPE: boolean vector type
// M1: 0 for exp, 1 for expm1
// BA: 0 for exp, 1 for 0.5*exp, 2 for pow(2,x), 10 for pow(10,x)

#if 1  // choose method

// Taylor expansion
template<class VTYPE, class BTYPE, int M1, int BA> 
static inline VTYPE exp_d(VTYPE const & initial_x) {    

    // Taylor coefficients, 1/n!
    // Not using minimax approximation because we prioritize precision close to x = 0
    const double p2  = 1./2.;
    const double p3  = 1./6.;
    const double p4  = 1./24.;
    const double p5  = 1./120.; 
    const double p6  = 1./720.; 
    const double p7  = 1./5040.; 
    const double p8  = 1./40320.; 
    const double p9  = 1./362880.; 
    const double p10 = 1./3628800.; 
    const double p11 = 1./39916800.; 
    const double p12 = 1./479001600.; 
    const double p13 = 1./6227020800.; 

    // maximum abs(x), value depends on BA, defined below
    // The lower limit of x is slightly more restrictive than the upper limit.
    // We are specifying the lower limit, except for BA = 1 because it is not used for negative x
    double max_x;

    // data vectors
    VTYPE x, r, z, n2;
    BTYPE inrange;                               // boolean vector

    if (BA <= 1) { // exp(x)
        max_x = BA == 0 ? 708.39 : 709.7; // lower limit for 0.5*exp(x) is -707.6, but we are using 0.5*exp(x) only for positive x in hyperbolic functions
        const double ln2d_hi = 0.693145751953125;
        const double ln2d_lo = 1.42860682030941723212E-6;
        x  = initial_x;
        r  = round(initial_x*VM_LOG2E);
        // subtraction in two steps for higher precision
        x = nmul_add(r, ln2d_hi, x);             //  x -= r * ln2d_hi;
        x = nmul_add(r, ln2d_lo, x);             //  x -= r * ln2d_lo;
    }
    else if (BA == 2) { // pow(2,x)
        max_x = 1022.0;
        r  = round(initial_x);
        x  = initial_x - r;
        x *= VM_LN2;
    }
    else if (BA == 10) { // pow(10,x)
        max_x = 307.65;
        const double log10_2_hi = 0.30102999554947019;     // log10(2) in two parts
        const double log10_2_lo = 1.1451100899212592E-10;
        x  = initial_x;
        r  = round(initial_x*(VM_LOG2E*VM_LN10));
        // subtraction in two steps for higher precision
        x  = nmul_add(r, log10_2_hi, x);         //  x -= r * log10_2_hi;
        x  = nmul_add(r, log10_2_lo, x);         //  x -= r * log10_2_lo;
        x *= VM_LN10;
    }
    else  {  // undefined value of BA
        return 0.;
    }

    z = polynomial_13m(x, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13);

    if (BA == 1) r--;  // 0.5 * exp(x)

    // multiply by power of 2 
    n2 = vm_pow2n(r);

    if (M1 == 0) {
        // exp
        z = (z + 1.0) * n2;
    }
    else {
        // expm1
        z = mul_add(z, n2, n2 - 1.0);            // z = z * n2 + (n2 - 1.0);
    }

    // check for overflow
    inrange  = abs(initial_x) < max_x;
    // check for INF and NAN
    inrange &= is_finite(initial_x);

    if (horizontal_and(inrange)) {
        // fast normal path
        return z;
    }
    else {
        // overflow, underflow and NAN
        r = select(sign_bit(initial_x), 0.-M1, infinite_vec<VTYPE>()); // value in case of +/- overflow or INF
        z = select(inrange, z, r);                                     // +/- underflow
        z = select(is_nan(initial_x), initial_x, z);                   // NAN goes through
        return z;
    }
}

#else

// Pade expansion uses less code and fewer registers, but is slower
template<class VTYPE, class BTYPE, int M1, int BA> 
static inline VTYPE exp_d(VTYPE const & initial_x) {

    // define constants
    const double ln2p1   = 0.693145751953125;
    const double ln2p2   = 1.42860682030941723212E-6;
    const double log2e   = VM_LOG2E;
    const double max_exp = 708.39;
    // coefficients of pade polynomials
    const double P0exp = 9.99999999999999999910E-1;
    const double P1exp = 3.02994407707441961300E-2;
    const double P2exp = 1.26177193074810590878E-4;
    const double Q0exp = 2.00000000000000000009E0;
    const double Q1exp = 2.27265548208155028766E-1;
    const double Q2exp = 2.52448340349684104192E-3;
    const double Q3exp = 3.00198505138664455042E-6;

    VTYPE x, r, xx, px, qx, y, n2;               // data vectors
    BTYPE inrange;                               // boolean vector

    x = initial_x;
    r = round(initial_x*log2e);

    // subtraction in one step would gives loss of precision
    x -= r * ln2p1;
    x -= r * ln2p2;

    xx = x * x;

    // px = x * P(x^2).
    px = polynomial_2(xx, P0exp, P1exp, P2exp) * x;

    // Evaluate Q(x^2).
    qx = polynomial_3(xx, Q0exp, Q1exp, Q2exp, Q3exp);

    // e^x = 1 + 2*P(x^2)/( Q(x^2) - P(x^2) )
    y = (2.0 * px) / (qx - px);

    // Get 2^n in double.
    // n  = round_to_int64_limited(r);
    // n2 = exp2(n);
    n2 = vm_pow2n(r);  // this is faster

    if (M1 == 0) {
        // exp
        y = (y + 1.0) * n2;
    }
    else {
        // expm1
        y = y * n2 + (n2 - 1.0);
    }

    // overflow
    inrange  = abs(initial_x) < max_exp;
    // check for INF and NAN
    inrange &= is_finite(initial_x);

    if (horizontal_and(inrange)) {
        // fast normal path
        return y;
    }
    else {
        // overflow, underflow and NAN
        r = select(sign_bit(initial_x), 0.-M1, infinite_vec<VTYPE>()); // value in case of overflow or INF
        y = select(inrange, y, r);                                     // +/- overflow
        y = select(is_nan(initial_x), initial_x, y);                   // NAN goes through
        return y;
    }
}
#endif

// instances of exp_d template
static inline Vec2d exp(Vec2d const & x) {
    return exp_d<Vec2d, Vec2db, 0, 0>(x);
}

static inline Vec2d expm1(Vec2d const & x) {
    return exp_d<Vec2d, Vec2db, 1, 0>(x);
}

static inline Vec2d exp2(Vec2d const & x) {
    return exp_d<Vec2d, Vec2db, 0, 2>(x);
}

static inline Vec2d exp10(Vec2d const & x) {
    return exp_d<Vec2d, Vec2db, 0, 10>(x);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec4d exp(Vec4d const & x) {
    return exp_d<Vec4d, Vec4db, 0, 0>(x);
}

static inline Vec4d expm1(Vec4d const & x) {
    return exp_d<Vec4d, Vec4db, 1, 0>(x);
}

static inline Vec4d exp2(Vec4d const & x) {
    return exp_d<Vec4d, Vec4db, 0, 2>(x);
}

static inline Vec4d exp10(Vec4d const & x) {
    return exp_d<Vec4d, Vec4db, 0, 10>(x);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec8d exp(Vec8d const & x) {
    return exp_d<Vec8d, Vec8db, 0, 0>(x);
}

static inline Vec8d expm1(Vec8d const & x) {
    return exp_d<Vec8d, Vec8db, 1, 0>(x);
}

static inline Vec8d exp2(Vec8d const & x) {
    return exp_d<Vec8d, Vec8db, 0, 2>(x);
}

static inline Vec8d exp10(Vec8d const & x) {
    return exp_d<Vec8d, Vec8db, 0, 10>(x);
}

#endif // MAX_VECTOR_SIZE >= 512

// Template for exp function, single precision
// The limit of abs(x) is defined by max_x below
// This function does not produce denormals
// Template parameters:
// VTYPE: float vector type
// BTYPE: boolean vector type
// M1: 0 for exp, 1 for expm1
// BA: 0 for exp, 1 for 0.5*exp, 2 for pow(2,x), 10 for pow(10,x)

template<class VTYPE, class BTYPE, int M1, int BA> 
static inline VTYPE exp_f(VTYPE const & initial_x) {

    // Taylor coefficients
    const float P0expf   =  1.f/2.f;
    const float P1expf   =  1.f/6.f;
    const float P2expf   =  1.f/24.f;
    const float P3expf   =  1.f/120.f; 
    const float P4expf   =  1.f/720.f; 
    const float P5expf   =  1.f/5040.f; 

    VTYPE x, r, x2, z, n2;                       // data vectors        
    BTYPE inrange;                               // boolean vector

    // maximum abs(x), value depends on BA, defined below
    // The lower limit of x is slightly more restrictive than the upper limit.
    // We are specifying the lower limit, except for BA = 1 because it is not used for negative x
    float max_x;

    if (BA <= 1) { // exp(x)
        const float ln2f_hi  =  0.693359375f;
        const float ln2f_lo  = -2.12194440e-4f;
        max_x = (BA == 0) ? 87.3f : 89.0f;

        x = initial_x;
        r = round(initial_x*float(VM_LOG2E));
        x = nmul_add(r, VTYPE(ln2f_hi), x);      //  x -= r * ln2f_hi;
        x = nmul_add(r, VTYPE(ln2f_lo), x);      //  x -= r * ln2f_lo;
    }
    else if (BA == 2) {                          // pow(2,x)
        max_x = 126.f;
        r = round(initial_x);
        x = initial_x - r;
        x = x * (float)VM_LN2;
    }
    else if (BA == 10) {                         // pow(10,x)
        max_x = 37.9f;
        const float log10_2_hi = 0.301025391f;   // log10(2) in two parts
        const float log10_2_lo = 4.60503907E-6f;
        x = initial_x;
        r = round(initial_x*float(VM_LOG2E*VM_LN10));
        x = nmul_add(r, VTYPE(log10_2_hi), x);   //  x -= r * log10_2_hi;
        x = nmul_add(r, VTYPE(log10_2_lo), x);   //  x -= r * log10_2_lo;
        x = x * (float)VM_LN10;
    }
    else  {  // undefined value of BA
        return 0.;
    }

    x2 = x * x;
    z = polynomial_5(x,P0expf,P1expf,P2expf,P3expf,P4expf,P5expf);    
    z = mul_add(z, x2, x);                       // z *= x2;  z += x;

    if (BA == 1) r--;                            // 0.5 * exp(x)

    // multiply by power of 2 
    n2 = vm_pow2n(r);

    if (M1 == 0) {
        // exp
        z = (z + 1.0f) * n2;
    }
    else {
        // expm1
        z = mul_add(z, n2, n2 - 1.0f);           //  z = z * n2 + (n2 - 1.0f);
    }

    // check for overflow
    inrange  = abs(initial_x) < max_x;
    // check for INF and NAN
    inrange &= is_finite(initial_x);

    if (horizontal_and(inrange)) {
        // fast normal path
        return z;
    }
    else {
        // overflow, underflow and NAN
        r = select(sign_bit(initial_x), 0.f-M1, infinite_vec<VTYPE>()); // value in case of +/- overflow or INF
        z = select(inrange, z, r);                                      // +/- underflow
        z = select(is_nan(initial_x), initial_x, z);                    // NAN goes through
        return z;
    }
}

// instances of exp_f template
static inline Vec4f exp(Vec4f const & x) {
    return exp_f<Vec4f, Vec4fb, 0, 0>(x);
}

static inline Vec4f expm1(Vec4f const & x) {
    return exp_f<Vec4f, Vec4fb, 1, 0>(x);
}

static inline Vec4f exp2(Vec4f const & x) {
    return exp_f<Vec4f, Vec4fb, 0, 2>(x);
}

static inline Vec4f exp10(Vec4f const & x) {
    return exp_f<Vec4f, Vec4fb, 0, 10>(x);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec8f exp(Vec8f const & x) {
    return exp_f<Vec8f, Vec8fb, 0, 0>(x);
}

static inline Vec8f expm1(Vec8f const & x) {
    return exp_f<Vec8f, Vec8fb, 1, 0>(x);
}

static inline Vec8f exp2(Vec8f const & x) {
    return exp_f<Vec8f, Vec8fb, 0, 2>(x);
}

static inline Vec8f exp10(Vec8f const & x) {
    return exp_f<Vec8f, Vec8fb, 0, 10>(x);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec16f exp(Vec16f const & x) {
    return exp_f<Vec16f, Vec16fb, 0, 0>(x);
}

static inline Vec16f expm1(Vec16f const & x) {
    return exp_f<Vec16f, Vec16fb, 1, 0>(x);
}

static inline Vec16f exp2(Vec16f const & x) {
    return exp_f<Vec16f, Vec16fb, 0, 2>(x);
}

static inline Vec16f exp10(Vec16f const & x) {
    return exp_f<Vec16f, Vec16fb, 0, 10>(x);
}

#endif // MAX_VECTOR_SIZE >= 512


/******************************************************************************
*                 Logarithm functions
******************************************************************************/

// Helper functions: fraction_2(x) = fraction(x)*0.5

// Modified fraction function:
// Extract the fraction part of a floating point number, and divide by 2
// The fraction function is defined in vectorf128.h etc.
// fraction_2(x) = fraction(x)*0.5
// This version gives half the fraction without extra delay
// Does not work for x = 0
static inline Vec4f fraction_2(Vec4f const & a) {
    Vec4ui t1 = _mm_castps_si128(a);   // reinterpret as 32-bit integer
    Vec4ui t2 = Vec4ui((t1 & 0x007FFFFF) | 0x3F000000); // set exponent to 0 + bias
    return _mm_castsi128_ps(t2);
}

static inline Vec2d fraction_2(Vec2d const & a) {
    Vec2uq t1 = _mm_castpd_si128(a);   // reinterpret as 64-bit integer
    Vec2uq t2 = Vec2uq((t1 & 0x000FFFFFFFFFFFFFll) | 0x3FE0000000000000ll); // set exponent to 0 + bias
    return _mm_castsi128_pd(t2);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec8f fraction_2(Vec8f const & a) {
#if defined (VECTORI256_H) && VECTORI256_H > 2  // 256 bit integer vectors are available, AVX2
    Vec8ui t1 = _mm256_castps_si256(a);   // reinterpret as 32-bit integer
    Vec8ui t2 = (t1 & 0x007FFFFF) | 0x3F000000; // set exponent to 0 + bias
    return _mm256_castsi256_ps(t2);
#else
    return Vec8f(fraction_2(a.get_low()), fraction_2(a.get_high()));
#endif
}

static inline Vec4d fraction_2(Vec4d const & a) {
#if VECTORI256_H > 1  // AVX2
    Vec4uq t1 = _mm256_castpd_si256(a);   // reinterpret as 64-bit integer
    Vec4uq t2 = Vec4uq((t1 & 0x000FFFFFFFFFFFFFll) | 0x3FE0000000000000ll); // set exponent to 0 + bias
    return _mm256_castsi256_pd(t2);
#else
    return Vec4d(fraction_2(a.get_low()), fraction_2(a.get_high()));
#endif
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec16f fraction_2(Vec16f const & a) {
#if INSTRSET >= 9                    // 512 bit integer vectors are available, AVX512
    return _mm512_getmant_ps(a, _MM_MANT_NORM_p5_1, _MM_MANT_SIGN_zero);
    //return Vec16f(_mm512_getmant_ps(a, _MM_MANT_NORM_1_2, _MM_MANT_SIGN_zero)) * 0.5f;
#else
    return Vec16f(fraction_2(a.get_low()), fraction_2(a.get_high()));
#endif
}

static inline Vec8d fraction_2(Vec8d const & a) {
#if INSTRSET >= 9                    // 512 bit integer vectors are available, AVX512
    return _mm512_getmant_pd(a, _MM_MANT_NORM_p5_1, _MM_MANT_SIGN_zero);
    //return Vec8d(_mm512_getmant_pd(a, _MM_MANT_NORM_1_2, _MM_MANT_SIGN_zero)) * 0.5;
#else
    return Vec8d(fraction_2(a.get_low()), fraction_2(a.get_high()));
#endif
}

#endif // MAX_VECTOR_SIZE >= 512


// Helper functions: exponent_f(x) = exponent(x) as floating point number

union vm_ufi {
    float f;
    uint32_t i;
};

union vm_udi {
    double d;
    uint64_t i;
};

// extract exponent of a positive number x as a floating point number
static inline Vec4f exponent_f(Vec4f const & x) {
    const float pow2_23 =  8388608.0f;           // 2^23
    const float bias = 127.f;                    // bias in exponent
    const vm_ufi upow2_23 = {pow2_23};
    Vec4ui a = reinterpret_i(x);                 // bit-cast x to integer
    Vec4ui b = a >> 23;                          // shift down exponent to low bits
    Vec4ui c = b | Vec4ui(upow2_23.i);           // insert new exponent
    Vec4f  d = reinterpret_f(c);                 // bit-cast back to double
    Vec4f  e = d - (pow2_23 + bias);             // subtract magic number and bias
    return e;
}

static inline Vec2d exponent_f(Vec2d const & x) {
    const double pow2_52 = 4503599627370496.0;   // 2^52
    const double bias = 1023.0;                  // bias in exponent
    const vm_udi upow2_52 = {pow2_52};

    Vec2uq a = reinterpret_i(x);                 // bit-cast x to integer
    Vec2uq b = a >> 52;                          // shift down exponent to low bits
    Vec2uq c = b | Vec2uq(upow2_52.i);           // insert new exponent
    Vec2d  d = reinterpret_d(c);                 // bit-cast back to double
    Vec2d  e = d - (pow2_52 + bias);             // subtract magic number and bias
    return e;
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec8f exponent_f(Vec8f const & x) {
    const float pow2_23 =  8388608.0f;           // 2^23
    const float bias = 127.f;                    // bias in exponent
    const vm_ufi upow2_23 = {pow2_23};
    Vec8ui a = reinterpret_i(x);                 // bit-cast x to integer
    Vec8ui b = a >> 23;                          // shift down exponent to low bits
    Vec8ui c = b | Vec8ui(upow2_23.i);           // insert new exponent
    Vec8f  d = reinterpret_f(c);                 // bit-cast back to double
    Vec8f  e = d - (pow2_23 + bias);             // subtract magic number and bias
    return e;
} 

// extract exponent of a positive number x as a floating point number
static inline Vec4d exponent_f(Vec4d const & x) {
    const double pow2_52 = 4503599627370496.0;   // 2^52
    const double bias = 1023.0;                  // bias in exponent
    const vm_udi upow2_52 = {pow2_52};

    Vec4uq a = reinterpret_i(x);                 // bit-cast x to integer
    Vec4uq b = a >> 52;                          // shift down exponent to low bits
    Vec4uq c = b | Vec4uq(upow2_52.i);           // insert new exponent
    Vec4d  d = reinterpret_d(c);                 // bit-cast back to double
    Vec4d  e = d - (pow2_52 + bias);             // subtract magic number and bias
    return e;
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec16f exponent_f(Vec16f const & x) {
#if INSTRSET >= 9                                // AVX512
    return _mm512_getexp_ps(x);
#else
    return Vec16f(exponent_f(x.get_low()), exponent_f(x.get_high()));
#endif
} 

// extract exponent of a positive number x as a floating point number
static inline Vec8d exponent_f(Vec8d const & x) {
#if INSTRSET >= 9                                // AVX512
    return _mm512_getexp_pd(x);
#else
    return Vec8d(exponent_f(x.get_low()), exponent_f(x.get_high()));
#endif
}

#endif // MAX_VECTOR_SIZE >= 512


// log function, double precision
// template parameters:
// VTYPE: f.p. vector type
// BTYPE: boolean vector type
// M1: 0 for log, 1 for log1p
template<class VTYPE, class BTYPE, int M1> 
static inline VTYPE log_d(VTYPE const & initial_x) {

    // define constants
    const double ln2_hi =  0.693359375;
    const double ln2_lo = -2.121944400546905827679E-4;
    const double P0log  =  7.70838733755885391666E0;
    const double P1log  =  1.79368678507819816313E1;
    const double P2log  =  1.44989225341610930846E1;
    const double P3log  =  4.70579119878881725854E0;
    const double P4log  =  4.97494994976747001425E-1;
    const double P5log  =  1.01875663804580931796E-4;
    const double Q0log  =  2.31251620126765340583E1;
    const double Q1log  =  7.11544750618563894466E1;
    const double Q2log  =  8.29875266912776603211E1;
    const double Q3log  =  4.52279145837532221105E1;
    const double Q4log  =  1.12873587189167450590E1;

    VTYPE x1, x, x2, px, qx, res, fe;            // data vectors
    BTYPE blend, overflow, underflow;            // boolean vectors

    if (M1 == 0) {
        x1 = initial_x;                          // log(x)
    }
    else {
        x1 = initial_x + 1.0;                    // log(x+1)
    }
    // separate mantissa from exponent 
    // VTYPE x  = fraction(x1) * 0.5;
    x  = fraction_2(x1);
    fe = exponent_f(x1);

    blend = x > VM_SQRT2*0.5;
    x  = if_add(!blend, x, x);                   // conditional add
    fe = if_add(blend, fe, 1.);                  // conditional add

    if (M1 == 0) {
        // log(x). Expand around 1.0
        x -= 1.0;
    }
    else {
        // log(x+1). Avoid loss of precision when adding 1 and later subtracting 1 if exponent = 0
        x = select(fe==0., initial_x, x - 1.0);
    }

    // rational form 
    px  = polynomial_5 (x, P0log, P1log, P2log, P3log, P4log, P5log);
    x2  = x * x;
    px *= x * x2;
    qx  = polynomial_5n(x, Q0log, Q1log, Q2log, Q3log, Q4log);
    res = px / qx ;

    // add exponent
    res  = mul_add(fe, ln2_lo, res);             // res += fe * ln2_lo;
    res += nmul_add(x2, 0.5, x);                 // res += x  - 0.5 * x2;
    res  = mul_add(fe, ln2_hi, res);             // res += fe * ln2_hi;

    overflow  = !is_finite(x1);
    underflow = x1 < VM_SMALLEST_NORMAL;         // denormals not supported by this functions

    if (!horizontal_or(overflow | underflow)) {
        // normal path
        return res;
    }
    else {
        // overflow and underflow
        res = select(underflow, nan_vec<VTYPE>(NAN_LOG), res);                   // x1  < 0 gives NAN
        res = select(x1 == 0. || is_subnormal(x1), -infinite_vec<VTYPE>(), res); // x1 == 0 gives -INF
        res = select(overflow,  x1, res);                                        // INF or NAN goes through
        res = select(is_inf(x1)&sign_bit(x1), nan_vec<VTYPE>(NAN_LOG), res);     // -INF gives NAN
        return res;
    }
}


static inline Vec2d log(Vec2d const & x) {
    return log_d<Vec2d, Vec2db, 0>(x);
}

static inline Vec2d log1p(Vec2d const & x) {
    return log_d<Vec2d, Vec2db, 1>(x);
}

static inline Vec2d log2(Vec2d const & x) {
    return VM_LOG2E * log_d<Vec2d, Vec2db, 0>(x);
}

static inline Vec2d log10(Vec2d const & x) {
    return VM_LOG10E * log_d<Vec2d, Vec2db, 0>(x);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec4d log(Vec4d const & x) {
    return log_d<Vec4d, Vec4db, 0>(x);
}

static inline Vec4d log1p(Vec4d const & x) {
    return log_d<Vec4d, Vec4db, 1>(x);
}

static inline Vec4d log2(Vec4d const & x) {
    return VM_LOG2E * log_d<Vec4d, Vec4db, 0>(x);
}

static inline Vec4d log10(Vec4d const & x) {
    return VM_LOG10E * log_d<Vec4d, Vec4db, 0>(x);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec8d log(Vec8d const & x) {
    return log_d<Vec8d, Vec8db, 0>(x);
}

static inline Vec8d log1p(Vec8d const & x) {
    return log_d<Vec8d, Vec8db, 1>(x);
}

static inline Vec8d log2(Vec8d const & x) {
    return VM_LOG2E * log_d<Vec8d, Vec8db, 0>(x);
}

static inline Vec8d log10(Vec8d const & x) {
    return VM_LOG10E * log_d<Vec8d, Vec8db, 0>(x);
}

#endif // MAX_VECTOR_SIZE >= 512



// log function, single precision
// template parameters:
// VTYPE: f.p. vector type
// ITYPE: integer vector type with same element size
// BTYPE: boolean vector type
// BTYPEI: boolean vector type for ITYPE
// M1: 0 for log, 1 for log1p
template<class VTYPE, class ITYPE, class BTYPE, class BTYPEI, int M1> 
static inline VTYPE log_f(VTYPE const & initial_x) {

    // define constants
    const float ln2f_hi =  0.693359375f;
    const float ln2f_lo = -2.12194440E-4f;
    const float P0logf  =  3.3333331174E-1f;
    const float P1logf  = -2.4999993993E-1f;
    const float P2logf  =  2.0000714765E-1f;
    const float P3logf  = -1.6668057665E-1f;
    const float P4logf  =  1.4249322787E-1f;
    const float P5logf  = -1.2420140846E-1f;
    const float P6logf  =  1.1676998740E-1f;
    const float P7logf  = -1.1514610310E-1f;
    const float P8logf  =  7.0376836292E-2f;

    VTYPE x1, x, res, x2, fe;                    // data vectors
    ITYPE e;                                     // integer vector
    BTYPE blend, overflow, underflow;            // boolean vectors

    if (M1 == 0) {
        x1 = initial_x;                          // log(x)
    }
    else {
        x1 = initial_x + 1.0f;                   // log(x+1)
    }

    // separate mantissa from exponent 
    x = fraction_2(x1);
    e = exponent(x1);

    blend = x > float(VM_SQRT2*0.5);
    x  = if_add(!blend, x, x);                   // conditional add
    e  = if_add(BTYPEI(blend),  e, ITYPE(1));    // conditional add
    fe = to_float(e);

    if (M1 == 0) {
        // log(x). Expand around 1.0
        x -= 1.0f;
    }
    else {
        // log(x+1). Avoid loss of precision when adding 1 and later subtracting 1 if exponent = 0
        x = select(BTYPE(e==0), initial_x, x - 1.0f);
    }

    // Taylor expansion
    res = polynomial_8(x, P0logf, P1logf, P2logf, P3logf, P4logf, P5logf, P6logf, P7logf, P8logf);
    x2  = x*x;
    res *= x2*x;

    // add exponent
    res  = mul_add(fe, ln2f_lo, res);            // res += ln2f_lo  * fe;
    res += nmul_add(x2, 0.5f, x);                // res += x - 0.5f * x2;
    res  = mul_add(fe, ln2f_hi, res);            // res += ln2f_hi  * fe;

    overflow  = !is_finite(x1);
    underflow = x1 < VM_SMALLEST_NORMALF;        // denormals not supported by this functions

    if (!horizontal_or(overflow | underflow)) {
        // normal path
        return res;
    }
    else {
        // overflow and underflow
        res = select(underflow, nan_vec<VTYPE>(NAN_LOG), res);                    // x1 < 0 gives NAN
        res = select(x1 == 0.f || is_subnormal(x1), -infinite_vec<VTYPE>(), res); // x1 == 0 or denormal gives -INF
        res = select(overflow,  x1, res);                                         // INF or NAN goes through
        res = select(is_inf(x1)&sign_bit(x1), nan_vec<VTYPE>(NAN_LOG), res);      // -INF gives NAN
        return res;
    }
}

static inline Vec4f log(Vec4f const & x) {
    return log_f<Vec4f, Vec4i, Vec4fb, Vec4ib, 0>(x);
}

static inline Vec4f log1p(Vec4f const & x) {
    return log_f<Vec4f, Vec4i, Vec4fb, Vec4ib, 1>(x);
}

static inline Vec4f log2(Vec4f const & x) {
    return float(VM_LOG2E) * log_f<Vec4f, Vec4i, Vec4fb, Vec4ib, 0>(x);
}

static inline Vec4f log10(Vec4f const & x) {
    return float(VM_LOG10E) * log_f<Vec4f, Vec4i, Vec4fb, Vec4ib, 0>(x);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec8f log(Vec8f const & x) {
    return log_f<Vec8f, Vec8i, Vec8fb, Vec8ib, 0>(x);
}

static inline Vec8f log1p(Vec8f const & x) {
    return log_f<Vec8f, Vec8i, Vec8fb, Vec8ib, 1>(x);
}

static inline Vec8f log2(Vec8f const & x) {
    return float(VM_LOG2E) * log_f<Vec8f, Vec8i, Vec8fb, Vec8ib, 0>(x);
}

static inline Vec8f log10(Vec8f const & x) {
    return float(VM_LOG10E) * log_f<Vec8f, Vec8i, Vec8fb, Vec8ib, 0>(x);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec16f log(Vec16f const & x) {
    return log_f<Vec16f, Vec16i, Vec16fb, Vec16ib, 0>(x);
}

static inline Vec16f log1p(Vec16f const & x) {
    return log_f<Vec16f, Vec16i, Vec16fb, Vec16ib, 1>(x);
}

static inline Vec16f log2(Vec16f const & x) {
    return float(VM_LOG2E) * log_f<Vec16f, Vec16i, Vec16fb, Vec16ib, 0>(x);
}

static inline Vec16f log10(Vec16f const & x) {
    return float(VM_LOG10E) * log_f<Vec16f, Vec16i, Vec16fb, Vec16ib, 0>(x);
}

#endif // MAX_VECTOR_SIZE >= 512


/******************************************************************************
*           Cube root and reciprocal cube root
******************************************************************************/

// cube root template, double precision
// template parameters:
// VTYPE:  f.p. vector type
// ITYPE:  uint32_t integer vector type with same total number of bits
// ITYPE2: uint64_t integer vector type with same total number of bits
// BTYPE:  boolean vector type
// CR:     -1 for reciprocal cube root, 1 for cube root, 2 for cube root squared
template<class VTYPE, class ITYPE, class ITYPE2, class BTYPE, int CR> 
static inline VTYPE cbrt_d(VTYPE const & x) {
    const int iter = 7;     // iteration count of x^(-1/3) loop
    int i;
    VTYPE  xa, xa3, a, a2;
    ITYPE  m1, m2;
    BTYPE  underflow;
    ITYPE2 q1(0x5540000000000000ULL);            // exponent bias
    ITYPE2 q2(0x0005555500000000ULL);            // exponent multiplier for 1/3
    ITYPE2 q3(0x0010000000000000ULL);            // denormal limit
    const double one_third  = 1./3.;
    const double four_third = 4./3.;

    xa  = abs(x);
    xa3 = one_third*xa;

    // multiply exponent by -1/3
    m1 = reinterpret_i(xa);
    m2 = ITYPE(q1) - (m1 >> 20) * ITYPE(q2);
    a  = reinterpret_d(m2);
    underflow = BTYPE(ITYPE2(m1) < q3);          // true if denormal or zero

    // Newton Raphson iteration
    for (i = 0; i < iter-1; i++) {
        a2 = a * a;
        a = nmul_add(xa3, a2*a2, four_third*a);  // a = four_third*a - xa3*a2*a2;
    }
    // last iteration with better precision
    a2 = a * a;    
    a = mul_add(one_third, nmul_add(xa, a2*a2, a), a); // a = a + one_third*(a - xa*a2*a2);

    if (CR == -1) {  // reciprocal cube root
        // (note: gives wrong sign when input is INF)
        // generate INF if underflow
        a = select(underflow, infinite_vec<VTYPE>(), a);
        // get sign
        a = sign_combine(a, x);
    }
    else if (CR == 1) {     // cube root
        a = a * a * x;
        // generate 0 if underflow
        a = select(underflow, 0., a);
    }
    else if (CR == 2) {     // cube root squared
        // (note: gives wrong sign when input is INF)
        a = a * xa;
        // generate 0 if underflow
        a = select(underflow, 0., a);
    }
    return a;
}

// template instances for cbrt and reciprocal_cbrt

// cube root
static inline Vec2d cbrt(Vec2d const & x) {
    return cbrt_d<Vec2d, Vec4ui, Vec2uq, Vec2db, 1> (x);
}

// reciprocal cube root
static inline Vec2d reciprocal_cbrt(Vec2d const & x) {
    return cbrt_d<Vec2d, Vec4ui, Vec2uq, Vec2db, -1> (x);
}

// square cube root
static inline Vec2d square_cbrt(Vec2d const & x) {
    return cbrt_d<Vec2d, Vec4ui, Vec2uq, Vec2db, 2> (x);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec4d cbrt(Vec4d const & x) {
    return cbrt_d<Vec4d, Vec8ui, Vec4uq, Vec4db, 1> (x);
}

static inline Vec4d reciprocal_cbrt(Vec4d const & x) {
    return cbrt_d<Vec4d, Vec8ui, Vec4uq, Vec4db, -1> (x);
}

static inline Vec4d square_cbrt(Vec4d const & x) {
    return cbrt_d<Vec4d, Vec8ui, Vec4uq, Vec4db, 2> (x);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec8d cbrt(Vec8d const & x) {
    return cbrt_d<Vec8d, Vec16ui, Vec8uq, Vec8db, 1> (x);
}

static inline Vec8d reciprocal_cbrt(Vec8d const & x) {
    return cbrt_d<Vec8d, Vec16ui, Vec8uq, Vec8db, -1> (x);
}

static inline Vec8d square_cbrt(Vec8d const & x) {
    return cbrt_d<Vec8d, Vec16ui, Vec8uq, Vec8db, 2> (x);
}

#endif // MAX_VECTOR_SIZE >= 512


// cube root template, single precision
// template parameters:
// VTYPE:  f.p. vector type
// ITYPE:  uint32_t integer vector type
// BTYPE:  boolean vector type
// CR:     -1 for reciprocal cube root, 1 for cube root, 2 for cube root squared
template<class VTYPE, class ITYPE, class BTYPE, int CR> 
static inline VTYPE cbrt_f(VTYPE const & x) {

    const int iter = 6;                          // iteration count of x^(-1/3) loop
    int i;
    VTYPE  xa, xa3, a, a2;
    ITYPE  m1, m2;
    BTYPE  underflow;
    ITYPE  q1(0x54800000U);                      // exponent bias
    ITYPE  q2(0x002AAAAAU);                      // exponent multiplier for 1/3
    ITYPE  q3(0x00800000U);                      // denormal limit
    const  float one_third  = float(1./3.);
    const  float four_third = float(4./3.);

    xa  = abs(x);
    xa3 = one_third*xa;

    // multiply exponent by -1/3
    m1 = reinterpret_i(xa);
    m2 = q1 - (m1 >> 23) * q2;
    a  = reinterpret_f(m2);

    underflow = BTYPE(m1 < q3);                  // true if denormal or zero

    // Newton Raphson iteration
    for (i = 0; i < iter-1; i++) {
        a2 = a*a;        
        a = nmul_add(xa3, a2*a2, four_third*a);  // a = four_third*a - xa3*a2*a2;
    }
    // last iteration with better precision
    a2 = a*a;    
    a = mul_add(one_third, nmul_add(xa, a2*a2, a), a); //a = a + one_third*(a - xa*a2*a2);

    if (CR == -1) {                              // reciprocal cube root
        // generate INF if underflow
        a = select(underflow, infinite_vec<VTYPE>(), a);
        // get sign
        a = sign_combine(a, x);
    }
    else if (CR == 1) {                          // cube root
        a = a * a * x;
        // generate 0 if underflow
        a = select(underflow, 0., a);
    }
    else if (CR == 2) {                          // cube root squared
        a = a * xa;
        // generate 0 if underflow
        a = select(underflow, 0., a);
    }
    return a;
}

// template instances for cbrt and reciprocal_cbrt

// cube root
static inline Vec4f cbrt(Vec4f const & x) {
    return cbrt_f<Vec4f, Vec4ui, Vec4fb, 1> (x);
}

// reciprocal cube root
static inline Vec4f reciprocal_cbrt(Vec4f const & x) {
    return cbrt_f<Vec4f, Vec4ui, Vec4fb, -1> (x);
}

// square cube root
static inline Vec4f square_cbrt(Vec4f const & x) {
    return cbrt_f<Vec4f, Vec4ui, Vec4fb, 2> (x);
}

#if MAX_VECTOR_SIZE >= 256

static inline Vec8f cbrt(Vec8f const & x) {
    return cbrt_f<Vec8f, Vec8ui, Vec8fb, 1> (x);
}

static inline Vec8f reciprocal_cbrt(Vec8f const & x) {
    return cbrt_f<Vec8f, Vec8ui, Vec8fb, -1> (x);
}

static inline Vec8f square_cbrt(Vec8f const & x) {
    return cbrt_f<Vec8f, Vec8ui, Vec8fb, 2> (x);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

static inline Vec16f cbrt(Vec16f const & x) {
    return cbrt_f<Vec16f, Vec16ui, Vec16fb, 1> (x);
}

static inline Vec16f reciprocal_cbrt(Vec16f const & x) {
    return cbrt_f<Vec16f, Vec16ui, Vec16fb, -1> (x);
}

static inline Vec16f square_cbrt(Vec16f const & x) {
    return cbrt_f<Vec16f, Vec16ui, Vec16fb, 2> (x);
}

#endif // MAX_VECTOR_SIZE >= 512


// ****************************************************************************
//                pow template, double precision
// ****************************************************************************
// Calculate x to the power of y.

// Precision is important here because rounding errors get multiplied by y.
// The logarithm is calculated with extra precision, and the exponent is 
// calculated separately.
// The logarithm is calculated by Pad\E9 approximation with 6'th degree 
// polynomials. A 7'th degree would be preferred for best precision by high y.
// The alternative method: log(x) = z + z^3*R(z)/S(z), where z = 2(x-1)/(x+1)
// did not give better precision.

// Template parameters:
// VTYPE:  data vector type
// ITYPE:  signed integer vector type
// BTYPE:  boolean vector type
template <class VTYPE, class ITYPE, class BTYPE>
static inline VTYPE pow_template_d(VTYPE const & x0, VTYPE const & y) {

    // define constants
    const double ln2d_hi = 0.693145751953125;           // log(2) in extra precision, high bits
    const double ln2d_lo = 1.42860682030941723212E-6;   // low bits of log(2)
    const double log2e   = VM_LOG2E;                    // 1/log(2)
    const double pow2_52 = 4503599627370496.0;          // 2^52

    // coefficients for Pad\E9 polynomials
    const double P0logl =  2.0039553499201281259648E1;
    const double P1logl =  5.7112963590585538103336E1;
    const double P2logl =  6.0949667980987787057556E1;
    const double P3logl =  2.9911919328553073277375E1;
    const double P4logl =  6.5787325942061044846969E0;
    const double P5logl =  4.9854102823193375972212E-1;
    const double P6logl =  4.5270000862445199635215E-5;
    const double Q0logl =  6.0118660497603843919306E1;
    const double Q1logl =  2.1642788614495947685003E2;
    const double Q2logl =  3.0909872225312059774938E2;
    const double Q3logl =  2.2176239823732856465394E2;
    const double Q4logl =  8.3047565967967209469434E1;
    const double Q5logl =  1.5062909083469192043167E1;

    // Taylor coefficients for exp function, 1/n!
    const double p2  = 1./2.;
    const double p3  = 1./6.;
    const double p4  = 1./24.;
    const double p5  = 1./120.; 
    const double p6  = 1./720.; 
    const double p7  = 1./5040.; 
    const double p8  = 1./40320.; 
    const double p9  = 1./362880.; 
    const double p10 = 1./3628800.; 
    const double p11 = 1./39916800.; 
    const double p12 = 1./479001600.; 
    const double p13 = 1./6227020800.; 

    // data vectors
    VTYPE x, x1, x2;
    VTYPE px, qx, ef, yr, v, z, z1;
    VTYPE lg, lg1, lg2;
    VTYPE lgerr, x2err;
    VTYPE e1, e2, e3, ee;
    // integer vectors
    ITYPE ei, ej, yodd;
    // boolean vectors
    BTYPE blend, xzero, xnegative;
    BTYPE overflow, underflow, xfinite, yfinite, efinite;

    // remove sign
    x1 = abs(x0);

    // Separate mantissa from exponent 
    // This gives the mantissa * 0.5
    x  = fraction_2(x1);

    // reduce range of x = +/- sqrt(2)/2
    blend = x > VM_SQRT2*0.5;
    x  = if_add(!blend, x, x);                   // conditional add

    // Pade approximation
    // Higher precision than in log function. Still higher precision wanted
    x -= 1.0;
    x2 = x*x;
    px = polynomial_6  (x, P0logl, P1logl, P2logl, P3logl, P4logl, P5logl, P6logl);
    px *= x * x2;
    qx = polynomial_6n (x, Q0logl, Q1logl, Q2logl, Q3logl, Q4logl, Q5logl);
    lg1 = px / qx;
 
    // extract exponent
    ef = exponent_f(x1);
    ef = if_add(blend, ef, 1.);                  // conditional add

    // multiply exponent by y
    // nearest integer e1 goes into exponent of result, remainder yr is added to log
    e1 = round(ef * y);
    yr = mul_sub_x(ef, y, e1);                   // calculate remainder yr. precision very important here

    // add initial terms to Pade expansion
    lg = nmul_add(0.5, x2, x) + lg1;             // lg = (x - 0.5 * x2) + lg1;
    // calculate rounding errors in lg
    // rounding error in multiplication 0.5*x*x
    x2err = mul_sub_x(0.5*x, x, 0.5*x2);
    // rounding error in additions and subtractions
    lgerr = mul_add(0.5, x2, lg - x) - lg1;      // lgerr = ((lg - x) + 0.5 * x2) - lg1;

    // extract something for the exponent
    e2 = round(lg * y * VM_LOG2E);
    // subtract this from lg, with extra precision
    v = mul_sub_x(lg, y, e2 * ln2d_hi);
    v = nmul_add(e2, ln2d_lo, v);                // v -= e2 * ln2d_lo;

    // add remainder from ef * y
    v = mul_add(yr, VM_LN2, v);                  // v += yr * VM_LN2;

    // correct for previous rounding errors
    v = nmul_add(lgerr + x2err, y, v);           // v -= (lgerr + x2err) * y;

    // exp function

    // extract something for the exponent if possible
    x = v;
    e3 = round(x*log2e);
    // high precision multiplication not needed here because abs(e3) <= 1
    x = nmul_add(e3, VM_LN2, x);                 // x -= e3 * VM_LN2;

    z = polynomial_13m(x, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13);
    z = z + 1.0;

    // contributions to exponent
    ee = e1 + e2 + e3;
    ei = round_to_int64_limited(ee);
    // biased exponent of result:
    ej = ei + (ITYPE(reinterpret_i(z)) >> 52);
    // check exponent for overflow and underflow
    overflow  = BTYPE(ej >= 0x07FF) | (ee >  3000.);
    underflow = BTYPE(ej <= 0x0000) | (ee < -3000.);

    // add exponent by integer addition
    z = reinterpret_d(ITYPE(reinterpret_i(z)) + (ei << 52));

    // check for special cases
    xfinite   = is_finite(x0);
    yfinite   = is_finite(y);
    efinite   = is_finite(ee);
    xzero     = is_zero_or_subnormal(x0);
    xnegative = x0  < 0.;

    // check for overflow and underflow
    if (horizontal_or(overflow | underflow)) {
        // handle errors
        z = select(underflow, VTYPE(0.), z);
        z = select(overflow, infinite_vec<VTYPE>(), z);
    }

    // check for x == 0
    z = select(xzero, select(y < 0., infinite_vec<VTYPE>(), select(y == 0., VTYPE(1.), VTYPE(0.))), z);

    // check for x < 0. y must be integer
    if (horizontal_or(xnegative)) {
        // test if y odd
        yodd = ITYPE(reinterpret_i(abs(y) + pow2_52)) << 63;     // convert y to integer and shift bit 0 to position of sign bit
        z1 = z | (x0 & VTYPE(reinterpret_d(yodd)));              // apply sign if y odd
        z1 = select(y == round(y), z1, nan_vec<VTYPE>(NAN_POW)); // NAN if y not integer
        z = select(xnegative, z1, z);
    }

    // check for range errors
    if (horizontal_and(xfinite & yfinite & efinite)) {
        // fast return if no special cases
        return z;
    }
    // handle special error cases
    z = select(yfinite & efinite, z, select(x1 == 1., VTYPE(1.), select((x1 > 1.) ^ sign_bit(y), infinite_vec<VTYPE>(), 0.)));
    yodd = ITYPE(reinterpret_i(abs(y) + pow2_52)) << 63; // same as above
    z = select(xfinite, z, select(y == 0., VTYPE(1.), select(y < 0., VTYPE(0.), infinite_vec<VTYPE>() | ( VTYPE(reinterpret_d(yodd)) & x0))));
    z = select(is_nan(x0), select(is_nan(y), x0 | y, x0), select(is_nan(y), y, z));
    return z;
}; 


//This template is in vectorf128.h to prevent implicit conversion of float y to int when float version is not defined:
//template <typename TT> static Vec2d pow(Vec2d const & a, TT n);

// instantiations of pow_template_d:
template <>
inline Vec2d pow<Vec2d const &>(Vec2d const & x, Vec2d const & y) {
    return pow_template_d<Vec2d, Vec2q, Vec2db>(x, y);
}

template <>
inline Vec2d pow<double>(Vec2d const & x, double y) {
    return pow_template_d<Vec2d, Vec2q, Vec2db>(x, y);
}
template <>
inline Vec2d pow<float>(Vec2d const & x, float y) {
    return pow_template_d<Vec2d, Vec2q, Vec2db>(x, (double)y);
}

#if MAX_VECTOR_SIZE >= 256

template <>
inline Vec4d pow<Vec4d const &>(Vec4d const & x, Vec4d const & y) {
    return pow_template_d<Vec4d, Vec4q, Vec4db>(x, y);
}

template <>
inline Vec4d pow<double>(Vec4d const & x, double y) {
    return pow_template_d<Vec4d, Vec4q, Vec4db>(x, y);
}

template <>
inline Vec4d pow<float>(Vec4d const & x, float y) {
    return pow_template_d<Vec4d, Vec4q, Vec4db>(x, (double)y);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

template <>
inline Vec8d pow<Vec8d const &>(Vec8d const & x, Vec8d const & y) {
    return pow_template_d<Vec8d, Vec8q, Vec8db>(x, y);
}

template <>
inline Vec8d pow<double>(Vec8d const & x, double y) {
    return pow_template_d<Vec8d, Vec8q, Vec8db>(x, y);
}

template <>
inline Vec8d pow<float>(Vec8d const & x, float y) {
    return pow_template_d<Vec8d, Vec8q, Vec8db>(x, (double)y);
}

#endif // MAX_VECTOR_SIZE >= 512



// ****************************************************************************
//                pow template, single precision
// ****************************************************************************

// Template parameters:
// VTYPE:  data vector type
// ITYPE:  signed integer vector type
// BTYPE:  boolean vector type
// Calculate x to the power of y
template <class VTYPE, class ITYPE, class BTYPE>
static inline VTYPE pow_template_f(VTYPE const & x0, VTYPE const & y) {

    // define constants
    const float ln2f_hi  =  0.693359375f;
    const float ln2f_lo  = -2.12194440e-4f;
    //const float max_expf =  87.3f;
    const float log2e    =  float(VM_LOG2E);     // 1/log(2)
    const float pow2_23  =  8388608.0f;          // 2^23

    const float P0logf  =  3.3333331174E-1f;
    const float P1logf  = -2.4999993993E-1f;
    const float P2logf  =  2.0000714765E-1f;
    const float P3logf  = -1.6668057665E-1f;
    const float P4logf  =  1.4249322787E-1f;
    const float P5logf  = -1.2420140846E-1f;
    const float P6logf  =  1.1676998740E-1f;
    const float P7logf  = -1.1514610310E-1f;
    const float P8logf  =  7.0376836292E-2f;

    // Taylor coefficients for exp function, 1/n!
    const float p2expf   =  1.f/2.f;
    const float p3expf   =  1.f/6.f;
    const float p4expf   =  1.f/24.f;
    const float p5expf   =  1.f/120.f; 
    const float p6expf   =  1.f/720.f; 
    const float p7expf   =  1.f/5040.f; 

    // data vectors
    VTYPE x, x1, x2;
    VTYPE ef, yr, v, z, z1;
    VTYPE lg, lg1;
    VTYPE lgerr, x2err;
    VTYPE e1, e2, e3, ee;
    // integer vectors
    ITYPE ei, ej, yodd;
    // boolean vectors
    BTYPE blend, xzero, xnegative;
    BTYPE overflow, underflow, xfinite, yfinite, efinite;

    // remove sign
    x1 = abs(x0);

    // Separate mantissa from exponent 
    // This gives the mantissa * 0.5
    x  = fraction_2(x1);

    // reduce range of x = +/- sqrt(2)/2
    blend = x > float(VM_SQRT2 * 0.5);
    x  = if_add(!blend, x, x);                   // conditional add

    // Taylor expansion, high precision
    x   -= 1.0f;
    x2   = x * x;
    lg1  = polynomial_8(x, P0logf, P1logf, P2logf, P3logf, P4logf, P5logf, P6logf, P7logf, P8logf);
    lg1 *= x2 * x; 
 
    // extract exponent
    ef = exponent_f(x1);
    ef = if_add(blend, ef, 1.0f);                // conditional add

    // multiply exponent by y
    // nearest integer e1 goes into exponent of result, remainder yr is added to log
    e1 = round(ef * y);
    yr = mul_sub_x(ef, y, e1);                   // calculate remainder yr. precision very important here

    // add initial terms to expansion
    lg = nmul_add(0.5f, x2, x) + lg1;            // lg = (x - 0.5f * x2) + lg1;

    // calculate rounding errors in lg
    // rounding error in multiplication 0.5*x*x
    x2err = mul_sub_x(0.5f*x, x, 0.5f * x2);
    // rounding error in additions and subtractions
    lgerr = mul_add(0.5f, x2, lg - x) - lg1;     // lgerr = ((lg - x) + 0.5f * x2) - lg1;

    // extract something for the exponent
    e2 = round(lg * y * float(VM_LOG2E));
    // subtract this from lg, with extra precision
    v = mul_sub_x(lg, y, e2 * ln2f_hi);
    v = nmul_add(e2, ln2f_lo, v);                // v -= e2 * ln2f_lo;

    // correct for previous rounding errors
    v -= mul_sub(lgerr + x2err, y, yr * float(VM_LN2)); // v -= (lgerr + x2err) * y - yr * float(VM_LN2) ;

    // exp function

    // extract something for the exponent if possible
    x = v;
    e3 = round(x*log2e);
    // high precision multiplication not needed here because abs(e3) <= 1
    x = nmul_add(e3, float(VM_LN2), x);          // x -= e3 * float(VM_LN2);

    // Taylor polynomial
    x2  = x  * x;
    z = polynomial_5(x, p2expf, p3expf, p4expf, p5expf, p6expf, p7expf)*x2 + x + 1.0f;

    // contributions to exponent
    ee = e1 + e2 + e3;
    ei = round_to_int(ee);
    // biased exponent of result:
    ej = ei + (ITYPE(reinterpret_i(z)) >> 23);
    // check exponent for overflow and underflow
    overflow  = BTYPE(ej >= 0x0FF) | (ee >  300.f);
    underflow = BTYPE(ej <= 0x000) | (ee < -300.f);

    // add exponent by integer addition
    z = reinterpret_f(ITYPE(reinterpret_i(z)) + (ei << 23)); // the extra 0x10000 is shifted out here

    // check for special cases
    xfinite   = is_finite(x0);
    yfinite   = is_finite(y);
    efinite   = is_finite(ee);
    xzero     = is_zero_or_subnormal(x0);
    xnegative = x0  < 0.f;

    // check for overflow and underflow
    if (horizontal_or(overflow | underflow)) {
        // handle errors
        z = select(underflow, VTYPE(0.f), z);
        z = select(overflow, infinite_vec<VTYPE>(), z);
    }

    // check for x == 0
    z = select(xzero, select(y < 0.f, infinite_vec<VTYPE>(), select(y == 0.f, VTYPE(1.), VTYPE(0.f))), z);

    // check for x < 0. y must be integer
    if (horizontal_or(xnegative)) {
        // test if y odd
        yodd = ITYPE(reinterpret_i(abs(y) + pow2_23)) << 31;     // convert y to integer and shift bit 0 to position of sign bit
        z1 = z | (x0 & VTYPE(reinterpret_f(yodd)));              // apply sign if y odd
        z1 = select(y == round(y), z1, nan_vec<VTYPE>(NAN_POW)); // NAN if y not integer
        z = select(xnegative, z1, z);
    }

    // check for range errors
    if (horizontal_and(xfinite & yfinite & efinite)) {
        // fast return if no special cases
        return z;
    }
    // handle special error cases
    z = select(yfinite & efinite, z, select(x1 == 1.f, VTYPE(1.f), select((x1 > 1.f) ^ sign_bit(y), infinite_vec<VTYPE>(), 0.f)));
    yodd = ITYPE(reinterpret_i(abs(y) + pow2_23)) << 31; // same as above
    z = select(xfinite, z, select(y == 0.f, VTYPE(1.f), select(y < 0.f, VTYPE(0.f), infinite_vec<VTYPE>() | (VTYPE(reinterpret_f(yodd)) & x0))));
    z = select(is_nan(x0), select(is_nan(y), x0 | y, x0), select(is_nan(y), y, z));
    return z;
}

//This template is in vectorf128.h to prevent implicit conversion of float y to int when float version is not defined:
//template <typename TT> static Vec4f pow(Vec4f const & a, TT n);

template <>
inline Vec4f pow<Vec4f const &>(Vec4f const & x, Vec4f const & y) {
    return pow_template_f<Vec4f, Vec4i, Vec4fb>(x, y);
}

template <>
inline Vec4f pow<float>(Vec4f const & x, float y) {
    return pow_template_f<Vec4f, Vec4i, Vec4fb>(x, y);
}

template <>
inline Vec4f pow<double>(Vec4f const & x, double y) {
    return pow_template_f<Vec4f, Vec4i, Vec4fb>(x, (float)y);
}

#if MAX_VECTOR_SIZE >= 256

template <>
inline Vec8f pow<Vec8f const &>(Vec8f const & x, Vec8f const & y) {
    return pow_template_f<Vec8f, Vec8i,  Vec8fb>(x, y);
}

template <>
inline Vec8f pow<float>(Vec8f const & x, float y) {
    return pow_template_f<Vec8f, Vec8i,  Vec8fb>(x, y);
}
template <>
inline Vec8f pow<double>(Vec8f const & x, double y) {
    return pow_template_f<Vec8f, Vec8i,  Vec8fb>(x, (float)y);
}

#endif // MAX_VECTOR_SIZE >= 256

#if MAX_VECTOR_SIZE >= 512

inline Vec16f pow(Vec16f const & x, Vec16f const & y) {
    return pow_template_f<Vec16f, Vec16i,  Vec16fb>(x, y);
}

inline Vec16f pow(Vec16f const & x, float y) {
    return pow_template_f<Vec16f, Vec16i,  Vec16fb>(x, y);
}

inline Vec16f pow(Vec16f const & x, double y) {
    return pow_template_f<Vec16f, Vec16i,  Vec16fb>(x, (float)y);
}

#endif // MAX_VECTOR_SIZE >= 512


// *************************************************************
//             power function with rational exponent
// *************************************************************
// Power function with rational exponent: x^(a/b)
// Template must be defined as class to allow partial template specialization
template <int a, int b>
class Power_rational {
public:
    // overloaded member function for each vector type
    static Vec4f pow(Vec4f const & x) {
        Vec4f y = x;
        // negative x allowed when b odd or a even
        // (if a is even then either b is odd or a/b can be reduced, 
        // but we can check a even anyway at no cost to be sure)
        if (a == 0) return 1.f;
        if ((b | ~a) & 1) y = abs(y);
        y = ::pow(y, float(double(a)/double(b)));
        if (a & b & 1) y = sign_combine(y, x);          // apply sign if a and b both odd
        if ((a ^ b) >= 0) y = select(x == 0.f, 0.f, y); // zero allowed for positive a and b
        return y;
    }
    static Vec2d pow(Vec2d const & x) {
        Vec2d y = x;
        if (a == 0) return 1.;
        if ((b | ~a) & 1) y = abs(y);
        y = ::pow(y, double((long double)a/(long double)b));
        if (a & b & 1) y = sign_combine(y, x);
        if ((a ^ b) >= 0) y = select(x == 0., 0., y);
        return y;
    }
#if MAX_VECTOR_SIZE >= 256
    static Vec8f pow(Vec8f const & x) {
        Vec8f y = x;
        if (a == 0) return 1.f;
        if ((b | ~a) & 1) y = abs(y);
        y = ::pow(y, float(double(a)/double(b)));
        if (a & b & 1) y = sign_combine(y, x);
        if ((a ^ b) >= 0) y = select(x == 0.f, 0.f, y);
        return y;
    }
    static Vec4d pow(Vec4d const & x) {
        Vec4d y = x;
        if (a == 0) return 1.;
        if ((b | ~a) & 1) y = abs(y);
        y = ::pow(y, double((long double)a/(long double)b));
        if (a & b & 1) y = sign_combine(y, x);
        if ((a ^ b) >= 0) y = select(x == 0., 0., y);
        return y;
    }
#endif // MAX_VECTOR_SIZE >= 256
#if MAX_VECTOR_SIZE >= 512
    static Vec16f pow(Vec16f const & x) {
        Vec16f y = x;
        if (a == 0) return 1.f;
        if ((b | ~a) & 1) y = abs(y);
        y = ::pow(y, float(double(a)/double(b)));
        if (a & b & 1) y = sign_combine(y, x);
        if ((a ^ b) >= 0) y = select(x == 0.f, 0.f, y);
        return y;
    }
    static Vec8d pow(Vec8d const & x) {
        Vec8d y = x;
        if (a == 0) return 1.;
        if ((b | ~a) & 1) y = abs(y);
        y = ::pow(y, double((long double)a/(long double)b));
        if (a & b & 1) y = sign_combine(y, x);
        if ((a ^ b) >= 0) y = select(x == 0., 0., y);
        return y;
    }
#endif // MAX_VECTOR_SIZE >= 512
};

// partial specialization for b = 0
template<int a>
class Power_rational<a,0> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {return nan_vec<VTYPE>(NAN_LOG);}
};

// partial specialization for b = 1
template<int a>
class Power_rational<a,1> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {return pow_n<a>(x);}
};

// partial specialization for b = 2
template<int a>
class Power_rational<a,2> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {
        VTYPE y = pow_n<(a > 0 ? a/2 : (a-1)/2)>(x);
        if (a & 1) y *= sqrt(x);
        return y;
    }
};

// full specialization for a = 1, b = 2
template<>
class Power_rational<1,2> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {        
        return sqrt(x);
    }
};

// full specialization for a = -1, b = 2
template<>
class Power_rational<-1,2> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {        
        // (this is faster than iteration method on modern CPUs)
        return VTYPE(1.f) / sqrt(x);
    }
};

// partial specialization for b = 3
template<int a>
class Power_rational<a,3> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {
        VTYPE t;
        switch (a % 3) {
        case -2:
            t = reciprocal_cbrt(x);
            t *= t;
            if (a == -2) return t;
            t = t / pow_n<(-a-2)/3>(x);
            break;
        case -1:
            t = reciprocal_cbrt(x);
            if (a == -1) return t;
            t = t / pow_n<(-a-1)/3>(x);
            break;
        case  0:
            t = pow_n<a/3>(x);
            break;
        case  1:
            t = cbrt(x);
            if (a == 1) return t;
            t = t * pow_n<a/3>(x);
            break;
        case  2:
            t = square_cbrt(x);
            if (a == 2) return t;
            t = t * pow_n<a/3>(x);
            break;
        }
        return t;
    }
};

// partial specialization for b = 4
template<int a>
class Power_rational<a,4> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {
        VTYPE t, s1, s2;
        s1 = sqrt(x);
        if (a & 1) s2 = sqrt(s1);
        switch (a % 4) {
        case -3:
            t = s2 / pow_n<1+(-a)/4>(x);
            break;
        case -2:
            t = s1 / pow_n<1+(-a)/4>(x);
            break;
        case -1:
            if (a != -1) s2 *= pow_n<(-a)/4>(x);
            t = VTYPE(1.f) / s2;
            break;
        case  0: default:
            t = pow_n<a/4>(x);
            break;
        case  1:
            t = s2;
            if (a != 1) t *= pow_n<a/4>(x);
            break;
        case  2:
            t = s1;
            if (a != 2) t *= pow_n<a/4>(x);
            break;
        case  3:
            t = s1 * s2;
            if (a != 3) t *= pow_n<a/4>(x);
            break;
        }
        return t;
    }
};

// partial specialization for b = 6
template<int a>
class Power_rational<a,6> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {
        VTYPE t, s1, s2, s3;
        switch (a % 6) {
        case -5:
            t = reciprocal_cbrt(x);
            t = t * t * sqrt(t);
            if (a != -5) t /= pow_n<(-a)/6>(x);
            break;
        case -4:
            t = reciprocal_cbrt(x);
            t *= t;
            if (a != -4) t /= pow_n<(-a)/6>(x);
            break;
        case -3:
            t = pow_n<a/6>(x);
            t /= sqrt(x);
            break;
        case -2:
            t = reciprocal_cbrt(x);
            if (a != -2) t /= pow_n<(-a)/6>(x);
            break;
        case -1:
            t = sqrt(reciprocal_cbrt(x));
            if (a != -1) t /= pow_n<(-a)/6>(x);
            break;
        case  0: default:
            t = pow_n<a/6>(x);
            break;
        case  1:
            t = sqrt(cbrt(x));
            if (a != 1) t *= pow_n<a/6>(x);
            break;
        case  2:
            t = cbrt(x);
            if (a != 2) t *= pow_n<a/6>(x);
            break;
        case  3:
            t = sqrt(x);
            if (a != 3) t *= pow_n<a/6>(x);
            break;
        case  4:
            t = square_cbrt(x);
            if (a != 4) t *= pow_n<a/6>(x);
            break;
        case  5:
            t = cbrt(x);
            t = t * t * sqrt(t);
            if (a != 5) t *= pow_n<a/6>(x);
            break;
        }
        return t;
    }
};

// partial specialization for b = 8
template<int a>
class Power_rational<a,8> {
public:
    template<class VTYPE>
    static VTYPE pow(VTYPE const & x) {
        VTYPE t, s1, s2, s3;
        s1 = sqrt(x);               // x^(1/2)
        if (a & 3) s2 = sqrt(s1);   // x^(1/4)
        if (a & 1) s3 = sqrt(s2);   // x^(1/8)
        switch (a % 8) {
        case -7:
            t = s3 / pow_n<1+(-a)/8>(x);
            break;
        case -6:
            t = s2 / pow_n<1+(-a)/8>(x);
            break;
        case -5:
            t = s3 * (s2 / pow_n<1+(-a)/8>(x));
            break;
        case -4:
            t = s1 / pow_n<1+(-a)/8>(x);
            break;
        case -3:
            t = s3 * (s1 / pow_n<1+(-a)/8>(x));
            break;
        case -2:
            if (a != -2) s2 *= pow_n<(-a)/8>(x);
            t = VTYPE(1.f) / s2;
            break;
        case -1:
            if (a != -1) s3 *= pow_n<(-a)/8>(x);
            t = VTYPE(1.f) / s3;
            break;
        case  0: default:
            t = pow_n<a/8>(x);
            break;
        case  1:
            t = s3;
            if (a != 1) t *= pow_n<a/8>(x);
            break;
        case  2:
            t = s2;
            if (a != 2) t *= pow_n<a/8>(x);
            break;
        case  3:
            t = s2 * s3;
            if (a != 3) t *= pow_n<a/8>(x);
            break;
        case  4:
            t = s1;
            if (a != 4) t *= pow_n<a/8>(x);
            break;
        case  5:
            t = s1 * s3;
            if (a != 5) t *= pow_n<a/8>(x);
            break;
        case  6:
            t = s1 * s2;
            if (a != 6) t *= pow_n<a/8>(x);
            break;
        case  7:
            t = s2 * s3;
            if (a != 7) s1 *= pow_n<a/8>(x);
            t *= s1;
            break;

        }
        return t;
    }
};

// macro to call template class member function pow
#define pow_ratio(x, a, b) (Power_rational<(b)<0 ? -(a):(a), (b)<0 ? -(b):(b)> ().pow(x))


/******************************************************************************
*                 Detect NAN codes
*
* These functions return the code hidden in a NAN. The sign bit is ignored
******************************************************************************/

Vec4i nan_code(Vec4f const & x) {
    Vec4i  a = reinterpret_i(x);
    Vec4ib b = (a & 0x7F800000) == 0x7F800000;   // check if NAN/INF
    return a & 0x007FFFFF & Vec4i(b);            // isolate NAN code bits
}

// This function returns the code hidden in a NAN. The sign bit is ignored
Vec2q nan_code(Vec2d const & x) {
    Vec2q  a = reinterpret_i(x);
    Vec2q const m = 0x7FF0000000000000;
    Vec2q const n = 0x000FFFFFFFFFFFFF;
    Vec2qb b = (a & m) == m;                     // check if NAN/INF
    return a & n & Vec2q(b);                     // isolate NAN code bits
}

#if MAX_VECTOR_SIZE >= 256

// This function returns the code hidden in a NAN. The sign bit is ignored
Vec8i nan_code(Vec8f const & x) {
    Vec8i  a = reinterpret_i(x);
    Vec8ib b = (a & 0x7F800000) == 0x7F800000;   // check if NAN/INF
    return a & 0x007FFFFF & Vec8i(b);            // isolate NAN code bits
}

// This function returns the code hidden in a NAN. The sign bit is ignored
Vec4q nan_code(Vec4d const & x) {
    Vec4q  a = reinterpret_i(x);
    Vec4q const m = 0x7FF0000000000000;
    Vec4q const n = 0x000FFFFFFFFFFFFF;
    Vec4qb b = (a & m) == m;                     // check if NAN/INF
    return a & n & Vec4q(b);                     // isolate NAN code bits
}

#endif // MAX_VECTOR_SIZE >= 256 
#if MAX_VECTOR_SIZE >= 512

// This function returns the code hidden in a NAN. The sign bit is ignored
Vec16i nan_code(Vec16f const & x) {
    Vec16i  a = Vec16i(reinterpret_i(x));
    Vec16ib b = (a & 0x7F800000) == 0x7F800000;  // check if NAN/INF
    return a & 0x007FFFFF & Vec16i(b);           // isolate NAN code bits
}

// This function returns the code hidden in a NAN. The sign bit is ignored
Vec8q nan_code(Vec8d const & x) {
    Vec8q  a = Vec8q(reinterpret_i(x));
    Vec8q const m = 0x7FF0000000000000;
    Vec8q const n = 0x000FFFFFFFFFFFFF;
    Vec8qb b = (a & m) == m;                     // check if NAN/INF
    return a & n & Vec8q(b);                     // isolate NAN code bits
}

#endif // MAX_VECTOR_SIZE >= 512

#endif  // VECTORMATH_EXP_H