windows-server-2003/base/crts/fpw32/tran/pow.c

/***
*pow.c - raise to a power
*
*       Copyright (c) 1991-2001, Microsoft Corporation. All rights reserved.
*
*Purpose:
*
*Revision History:
*        8-15-91  GDP   written
*       12-20-91  GDP   support IEEE exceptions & denormals
*        1-11-92  GDP   special handling of small powers
*                       special handling of u1, u2 when cancellation occurs
*        3-22-92  GDP   changed handling of int exponents, pow(0, neg)
*                       added check to avoid internal overflow due to large y
*        6-23-92  GDP   adjusted special return values according to NCEG spec
*       02-06-95  JWM   Mac merge
*       02-07-95  JWM   powhlp() usage restored to Intel version.
*       10-07-97  RDL   Added IA64.
*
*******************************************************************************/

#include <math.h>
#include <trans.h>
#include <float.h>

#if defined(_M_IA64)
#pragma function(pow)
#endif

static double _reduce(double);

static double const a1[18] = {
    0.00000000000000000000e+000,   /* dummy element */
    1.00000000000000000000e+000,
    9.57603280698573646910e-001,
    9.17004043204671231754e-001,
    8.78126080186649741555e-001,
    8.40896415253714543073e-001,
    8.05245165974627154042e-001,
    7.71105412703970411793e-001,
    7.38413072969749655712e-001,
    7.07106781186547524436e-001,
    6.77127773468446364133e-001,
    6.48419777325504832961e-001,
    6.20928906036742024317e-001,
    5.94603557501360533344e-001,
    5.69394317378345826849e-001,
    5.45253866332628829604e-001,
    5.22136891213706920173e-001,
    5.00000000000000000000e-001
};

static double const a2[9] = {
    0.00000000000000000000e+000,   /* dummy element */
   -5.31259064517897172664e-017,
    1.47993596544271355242e-017,
    1.23056946577104753260e-017,
   -1.74014448683923461658e-017,
    3.84891771232354074073e-017,
    2.33103467084383453312e-017,
    4.45607092891542322377e-017,
    4.27717757045531499216e-017
};

static double const log2inv = 1.44269504088896340739e+0; //  1/log(2)
static double const K = 0.44269504088896340736e+0;

static double const p1 = 0.83333333333333211405e-1;
static double const p2 = 0.12500000000503799174e-1;
static double const p3 = 0.22321421285924258967e-2;
static double const p4 = 0.43445775672163119635e-3;

#define P(v) (((p4 * v + p3) * v + p2) * v + p1)

static double const q1 = 0.69314718055994529629e+0;
static double const q2 = 0.24022650695909537056e+0;
static double const q3 = 0.55504108664085595326e-1;
static double const q4 = 0.96181290595172416964e-2;
static double const q5 = 0.13333541313585784703e-2;
static double const q6 = 0.15400290440989764601e-3;
static double const q7 = 0.14928852680595608186e-4;

#define Q(w)   ((((((q7 * w + q6) * w + q5) * w + q4) * w + \
                              q3) * w + q2) * w + q1)


/*
 * Thresholds for over/underflow that results in an adjusted value
 * too big/small to be represented as a double. An infinity or 0
 * is delivered to the trap handler instead
 */

static _dbl const _ovfx ={SET_DBL(0x40e40000,0)}; // 16*log2(XMAX*2^IEEE_ADJ)
static _dbl const _uflx ={SET_DBL(0xc0e3fc00,0)}; // 16*log2(XMIN*2^(-IEEE_ADJ))

#define OVFX _ovfx.dbl
#define UFLX _uflx.dbl

#define INT_POW_LIMIT   128.0

static double ymax = 1e20;

static double _reduce(double x)
{
    return 0.0625 * _frnd( 16.0 * x);
}

/***
*double pow(double x, double y) - x raised to the power of y
*
*Purpose:
*   Calculate x^y
*   Algorithm from Cody & Waite
*
*Entry:
*
*Exit:
*
*Exceptions:
*
*   All 5 IEEE exceptions may occur
*
*******************************************************************************/
double pow(double x, double y)
{
    uintptr_t savedcw;
    int m,mprim;
    int p,pprim;
    int i,iw1;
    int iy;
    int newexp;
    double diw1;
    double sign;
    double g,z,bigz,v,rz,result;
    double u1,u2,y1,y2,w,w1,w2;
    double savedx;

    /* save user fp control word */
    savedcw = _maskfp();
    savedx = x;         // save original value of first argument

    if (_fpclass(y) & (_FPCLASS_NZ | _FPCLASS_PZ)) {
        RETURN(savedcw, 1.0);
    }

    /* Check for zero^y */
    if (_fpclass(x) & (_FPCLASS_NZ | _FPCLASS_PZ)) { /* x==0? */
        int type;

        type = _d_inttype(y);

        if (y < 0.0) {
            result = (type == _D_ODD ? _copysign(D_INF,x) : D_INF);
            return _except2(FP_Z,OP_POW,savedx,y,result,savedcw|ISW_ZERODIVIDE);
        }
        else if (y > 0.0) {
            result = (type == _D_ODD ? x : 0.0);
                RETURN(savedcw, result);
        }
    }

    /* check for infinity or NAN */
    if (IS_D_SPECIAL(x) || IS_D_SPECIAL(y)) {
      double absx = fabs(x);

      if (IS_D_SNAN(x) || IS_D_SNAN(y)) {
        return _except2(FP_I,OP_POW,savedx,y,_d_snan2(x,y),savedcw | (ISW_INVALID>>5) );
      }
      if (IS_D_QNAN(x) || IS_D_QNAN(y)){
        return _handle_qnan2(OP_POW,x,y,savedcw | (ISW_INVALID>>5) );
      }
      
      /* there is at least one infinite argument ... */
      if (_powhlp(x, y, &result)) { /* removed "<" 0. */
        return _except2(FP_I,OP_POW,savedx,y,result,savedcw | (ISW_INVALID>>5) );
      }
      RETURN(savedcw, result);
    }

    sign = 1.0;
    if (x < 0) {
        switch (_d_inttype(y)) {
        case _D_ODD: /* y is an odd integral value */
            sign = -1.0;
            /* NO BREAK */
        case _D_EVEN:
            x = -x;
            break;
        default: /* y is not an integral value */
            return _except2(FP_I,OP_POW,savedx,y,D_IND,savedcw|(ISW_INVALID>>5));
        }
    }

    //
    // This is here in order to prevent internal overflows
    // due to a large value of y
    // The following relation holds on overflow with a scaled
    // result out of range
    // (lg stands for log base 2)
    //      |y| * |lg(x)| > MAXEXP + IEEE_ADJUST        <=>
    //      |y| >  2560 / |lg(x)|
    // The values of lg(x) closer to 0 are:
    //      x                                    lg(x)
    //  3fefffffffffffff    (0,99...9)        -1.601e-16
    //  3ff0000000000000    (1.0)              0.0
    //  3ff0000000000001    (1.00...1)         3.203e-16
    //
    // So if |y| > 2560/1.6e-16 = 1.6e19 overflow occurs
    // We set ymax to 1e20 in order to have a safety margin
    //

    if (ABS(y) > ymax) {
        if (y < 0) {
            y = -y;
            //
            // this may cause an underflow
            // there is no problem with fp sw pollution because
            // a FP_U exception is going to be raised anyway.
            //
            x = 1.0 / x;
        }
        if (x > 1.0) {
            return _except2(FP_O | FP_P,OP_POW,savedx,y,sign*D_INF,savedcw|ISW_OVERFLOW);
        }
        else if (x < 1.0){
            return _except2(FP_U | FP_P,OP_POW,savedx,y,sign*0.0,savedcw|ISW_UNDERFLOW);
        }
        else {
            RETURN(savedcw, sign*1.0);
        }
    }


    /* determine m, g */
    g = _decomp(x, &m);


    /* handle small integer powers
     * for small integer powers this is faster that Cody&Waite's
     * algorithm, and yields better precision
     * Without this piece of code there was not enough precision
     * to satisfy all requirements of the 'paranoia' test.
     * We choose INT_POW_LIMIT such that (1) no overflow or underflow
     * occurs while computing bigz (g is in the range
     * [0.5, 1.0) or (1.0, 2.0] so INT_POW_LIMIT should be less than
     * approximately 10^3) and (2) no extraordinary loss of precision
     * occurs because of repeated multiplications (this practically
     * restricts the maximum INT_POW_LIMIT to 128).
     */

    if (y <= INT_POW_LIMIT &&
        _d_inttype(x) != _D_NOINT &&
        _d_inttype(y) != _D_NOINT &&
        y > 0.0 ) {

        iy = (int)y;
        mprim = m * iy;

        for (bigz=1 ; iy ; iy >>= 1, g *= g) {
            if (iy & 0x1)
                bigz *= g;
        }

        newexp = _get_exp(bigz) + mprim;
        if (newexp > MAXEXP + IEEE_ADJUST) {
            return _except2(FP_O | FP_P, OP_POW, savedx, y, sign*bigz*D_INF, savedcw);
        }
        if (newexp < MINEXP - IEEE_ADJUST) {
            return _except2(FP_U | FP_P, OP_POW, savedx, y, sign*bigz*0.0, savedcw);
        }

    }


    else {

        /* determine p using binary search */
        p = 1;
        if (g <= a1[9])
            p = 9;
        if (g <= a1[p+4])
            p += 4;
        if (g <= a1[p+2])
            p += 2;


        /* C&W's algorithm is not very accurate when m*16-p == 1,
         * because there is cancellation between u1 and u2.
         * Handle this separately.
         */
        if (ABS(m*16-p) == 1) {
            u1 = log(x) * log2inv;
            u2 = 0.0;
        }
        else {
            /* determine z */
            z = ( (g - a1[p+1]) - a2[(p+1)/2] ) / ( g + a1[p+1] );
            z += z;


            /* determine u2 */
            v = z * z;
            rz = P(v) * v * z;
            rz += K * rz;
            u2 = (rz + z * K) + z;

            u1 = (m * 16 - p) * 0.0625;
        }

        /* determine w1, w2 */
        y1 = _reduce(y);
        y2 = y - y1;
        w = u2 * y + u1 * y2;
        w1 = _reduce(w);
        w2 = w - w1;
        w = w1 + u1 * y1;
        w1 = _reduce(w);
        w2 += w - w1;
        w = _reduce(w2);
        diw1 = 16 * (w1 + w); /* iw1 might overflow here, so use diw1 */
        w2 -= w;

        if (diw1 > OVFX) {
            return _except2(FP_O | FP_P,OP_POW,savedx,y,sign*D_INF,savedcw | ISW_OVERFLOW);
        }
        if (diw1 < UFLX) {
            return _except2(FP_U | FP_P,OP_POW,savedx,y,sign*0.0,savedcw | ISW_UNDERFLOW);
        }

        iw1 = (int) diw1;        /* now it is safe to cast to int */


        /* make sure w2 <= 0 */
        if (w2 > 0) {
            iw1 += 1;
            w2 -= 0.0625;
        }

        /* determine mprim, pprim */
        i = iw1 < 0 ? 0 : 1;
        mprim = iw1 / 16 + i;
        pprim = 16 * mprim - iw1;

        /* determine 2^w2 */
        bigz = Q(w2) * w2;

        /* determine  final result */
        bigz = a1[pprim + 1] + a1[pprim + 1] * bigz;
        newexp = _get_exp(bigz) + mprim;
    }


    if (newexp > MAXEXP) {
        result = sign * _set_exp(bigz, newexp - IEEE_ADJUST);
        return _except2(FP_O | FP_P, OP_POW, savedx, y, sign*D_INF, savedcw|ISW_OVERFLOW);
    }
    if (newexp < MINEXP) {
        result = sign * _set_exp(bigz, newexp + IEEE_ADJUST);
        return _except2(FP_U | FP_P, OP_POW, savedx, y, sign*0.0, savedcw|ISW_UNDERFLOW);
    }

    result = sign * _set_exp(bigz, newexp);
    RETURN_INEXACT2(OP_POW, savedx, y, result, savedcw|ISW_INEXACT);

}