Logo Search packages:      
Sourcecode: mingw-w64 version File versions  Download package

pow.c

/*                                        pow.c
 *
 *    Power function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, z, pow();
 *
 * z = pow( x, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Computes x raised to the yth power.  Analytically,
 *
 *      x**y  =  exp( y log(x) ).
 *
 * Following Cody and Waite, this program uses a lookup table
 * of 2**-i/16 and pseudo extended precision arithmetic to
 * obtain an extra three bits of accuracy in both the logarithm
 * and the exponential.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE     -26,26       30000      4.2e-16      7.7e-17
 *    IEEE     0,8700       30000      1.5e-14      2.1e-15
 * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * pow overflow     x**y > MAXNUM      INFINITY
 * pow underflow   x**y < 1/MAXNUM       0.0
 * pow domain      x<0 and y noninteger  0.0
 *
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/

/*
Modified for mingw
2002-09-27 Danny Smith <dannysmith@users.sourceforge.net>
*/

#ifdef __MINGW32__
#include "cephes_mconf.h"
#else
#include "mconf.h"
static char fname[] = {"pow"};
#endif

#ifndef _SET_ERRNO
#define _SET_ERRNO(x)
#endif

#define SQRTH 0.70710678118654752440

#ifdef UNK
static uD P[4] = {
  { { 4.97778295871696322025E-1 } },
  { { 3.73336776063286838734E0 } },
  { { 7.69994162726912503298E0 } },
  { { 4.66651806774358464979E0 } }
};
static uD Q[4] = {
  { { 9.33340916416696166113E0 } },
  { { 2.79999886606328401649E1 } },
  { { 3.35994905342304405431E1 } },
  { { 1.39995542032307539578E1 } }
};
/* 2^(-i/16), IEEE precision */
static uD A[17] = {
  { { 1.00000000000000000000E0 } },
  { { 9.57603280698573700036E-1 } },
  { { 9.17004043204671215328E-1 } },
  { { 8.78126080186649726755E-1 } },
  { { 8.40896415253714502036E-1 } },
  { { 8.05245165974627141736E-1 } },
  { { 7.71105412703970372057E-1 } },
  { { 7.38413072969749673113E-1 } },
  { { 7.07106781186547572737E-1 } },
  { { 6.77127773468446325644E-1 } },
  { { 6.48419777325504820276E-1 } },
  { { 6.20928906036742001007E-1 } },
  { { 5.94603557501360513449E-1 } },
  { { 5.69394317378345782288E-1 } },
  { { 5.45253866332628844837E-1 } },
  { { 5.22136891213706877402E-1 } },
  { { 5.00000000000000000000E-1 } }
};
static uD B[9] = {
  { { 0.00000000000000000000E0 } },
  { { 1.64155361212281360176E-17 } },
  { { 4.09950501029074826006E-17 } },
  { { 3.97491740484881042808E-17 } },
  { { -4.83364665672645672553E-17 } },
  { { 1.26912513974441574796E-17 } },
  { { 1.99100761573282305549E-17 } },
  { { -1.52339103990623557348E-17 } },
  { { 0.00000000000000000000E0 } }
};
static uD R[7] = {
  { { 1.49664108433729301083E-5 } },
  { { 1.54010762792771901396E-4 } },
  { { 1.33335476964097721140E-3 } },
  { { 9.61812908476554225149E-3 } },
  { { 5.55041086645832347466E-2 } },
  { { 2.40226506959099779976E-1 } },
  { { 6.93147180559945308821E-1 } }
};

#define douba(k) A[k].d
#define doubb(k) B[k].d
#define MEXP 16383.0
#ifdef DENORMAL
#define MNEXP -17183.0
#else
#define MNEXP -16383.0
#endif
#endif

#ifdef IBMPC
static const uD P[4] = {
  { { 0x5cf0,0x7f5b,0xdb99,0x3fdf } },
  { { 0xdf15,0xea9e,0xddef,0x400d } },
  { { 0xeb6f,0x7f78,0xccbd,0x401e } },
  { { 0x9b74,0xb65c,0xaa83,0x4012 } }
};
static const uD Q[4] = {
  { { 0x914e,0x9b20,0xaab4,0x4022 } },
  { { 0xc9f5,0x41c1,0xffff,0x403b } },
  { { 0x6402,0x1b17,0xccbc,0x4040 } },
  { { 0xe92e,0x918a,0xffc5,0x402b } }
};
static const uD A[17] = {
  { { 0x0000,0x0000,0x0000,0x3ff0 } },
  { { 0x90da,0xa2a4,0xa4af,0x3fee } },
  { { 0xa487,0xdcfb,0x5818,0x3fed } },
  { { 0x529c,0xdd85,0x199b,0x3fec } },
  { { 0xd3ad,0x995a,0xe89f,0x3fea } },
  { { 0xf090,0x82a3,0xc491,0x3fe9 } },
  { { 0xa0db,0x422a,0xace5,0x3fe8 } },
  { { 0x0187,0x73eb,0xa114,0x3fe7 } },
  { { 0x3bcd,0x667f,0xa09e,0x3fe6 } },
  { { 0x5429,0xdd48,0xab07,0x3fe5 } },
  { { 0x2a27,0xd536,0xbfda,0x3fe4 } },
  { { 0x3422,0x4c12,0xdea6,0x3fe3 } },
  { { 0xb715,0x0a31,0x06fe,0x3fe3 } },
  { { 0x6238,0x6e75,0x387a,0x3fe2 } },
  { { 0x517b,0x3c7d,0x72b8,0x3fe1 } },
  { { 0x890f,0x6cf9,0xb558,0x3fe0 } },
  { { 0x0000,0x0000,0x0000,0x3fe0 } }
};
static const uD B[9] = {
  { { 0x0000,0x0000,0x0000,0x0000 } },
  { { 0x3707,0xd75b,0xed02,0x3c72 } },
  { { 0xcc81,0x345d,0xa1cd,0x3c87 } },
  { { 0x4b27,0x5686,0xe9f1,0x3c86 } },
  { { 0x6456,0x13b2,0xdd34,0xbc8b } },
  { { 0x42e2,0xafec,0x4397,0x3c6d } },
  { { 0x82e4,0xd231,0xf46a,0x3c76 } },
  { { 0x8a76,0xb9d7,0x9041,0xbc71 } },
  { { 0x0000,0x0000,0x0000,0x0000 } }
};
static const uD R[7] = {
  { { 0x937f,0xd7f2,0x6307,0x3eef } },
  { { 0x9259,0x60fc,0x2fbe,0x3f24 } },
  { { 0xef1d,0xc84a,0xd87e,0x3f55 } },
  { { 0x33b7,0x6ef1,0xb2ab,0x3f83 } },
  { { 0x1a92,0xd704,0x6b08,0x3fac } },
  { { 0xc56d,0xff82,0xbfbd,0x3fce } },
  { { 0x39ef,0xfefa,0x2e42,0x3fe6 } }
};

#define douba(k) (A[(k)].d)
#define doubb(k) (B[(k)].d)
#define MEXP 16383.0
#ifdef DENORMAL
#define MNEXP -17183.0
#else
#define MNEXP -16383.0
#endif
#endif

#ifdef MIEEE
static uD P[4] = {
  { { 0x3fdf,0xdb99,0x7f5b,0x5cf0 } },
  { { 0x400d,0xddef,0xea9e,0xdf15 } },
  { { 0x401e,0xccbd,0x7f78,0xeb6f } },
  { { 0x4012,0xaa83,0xb65c,0x9b74 } }
};
static uD Q[4] = {
  { { 0x4022,0xaab4,0x9b20,0x914e } },
  { { 0x403b,0xffff,0x41c1,0xc9f5 } },
  { { 0x4040,0xccbc,0x1b17,0x6402 } },
  { { 0x402b,0xffc5,0x918a,0xe92e } }
};
static unsigned short A[17] = {
  { { 0x3ff0,0x0000,0x0000,0x0000 } },
  { { 0x3fee,0xa4af,0xa2a4,0x90da } },
  { { 0x3fed,0x5818,0xdcfb,0xa487 } },
  { { 0x3fec,0x199b,0xdd85,0x529c } },
  { { 0x3fea,0xe89f,0x995a,0xd3ad } },
  { { 0x3fe9,0xc491,0x82a3,0xf090 } },
  { { 0x3fe8,0xace5,0x422a,0xa0db } },
  { { 0x3fe7,0xa114,0x73eb,0x0187 } },
  { { 0x3fe6,0xa09e,0x667f,0x3bcd } },
  { { 0x3fe5,0xab07,0xdd48,0x5429 } },
  { { 0x3fe4,0xbfda,0xd536,0x2a27 } },
  { { 0x3fe3,0xdea6,0x4c12,0x3422 } },
  { { 0x3fe3,0x06fe,0x0a31,0xb715 } },
  { { 0x3fe2,0x387a,0x6e75,0x6238 } },
  { { 0x3fe1,0x72b8,0x3c7d,0x517b } },
  { { 0x3fe0,0xb558,0x6cf9,0x890f } },
  { { 0x3fe0,0x0000,0x0000,0x0000 } }
};
static uD B[9] = {
  { { 0x0000,0x0000,0x0000,0x0000 } },
  { { 0x3c72,0xed02,0xd75b,0x3707 } },
  { { 0x3c87,0xa1cd,0x345d,0xcc81 } },
  { { 0x3c86,0xe9f1,0x5686,0x4b27 } },
  { { 0xbc8b,0xdd34,0x13b2,0x6456 } },
  { { 0x3c6d,0x4397,0xafec,0x42e2 } },
  { { 0x3c76,0xf46a,0xd231,0x82e4 } },
  { { 0xbc71,0x9041,0xb9d7,0x8a76 } },
  { { 0x0000,0x0000,0x0000,0x0000 } }
};
static uD R[7] = {
  { { 0x3eef,0x6307,0xd7f2,0x937f } },
  { { 0x3f24,0x2fbe,0x60fc,0x9259 } },
  { { 0x3f55,0xd87e,0xc84a,0xef1d } },
  { { 0x3f83,0xb2ab,0x6ef1,0x33b7 } },
  { { 0x3fac,0x6b08,0xd704,0x1a92 } },
  { { 0x3fce,0xbfbd,0xff82,0xc56d } },
  { { 0x3fe6,0x2e42,0xfefa,0x39ef } }
};

#define douba(k) (A[(k)].d)
#define doubb(k) (B[(k)].d)
#define MEXP 16383.0
#ifdef DENORMAL
#define MNEXP -17183.0
#else
#define MNEXP -16383.0
#endif
#endif

/* log2(e) - 1 */
#define LOG2EA 0.44269504088896340736

#define F W
#define Fa Wa
#define Fb Wb
#define G W
#define Ga Wa
#define Gb u
#define H W
#define Ha Wb
#define Hb Wb

#ifdef __MINGW32__
static __inline__ double reduc( double );
extern double __powi ( double, int );
extern double pow ( double x,  double y);

#else /* __MINGW32__ */

extern double floor ( double );
extern double fabs ( double );
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, uD *, int );
extern double __powi ( double, int );
extern int signbit ( double );
extern int isnan ( double );
extern int isfinite ( double );
static double reduc ( double );
extern double MAXNUM;
#ifdef INFINITIES
extern double INFINITY;
#endif
#ifdef NANS
extern double NAN;
#endif
#ifdef MINUSZERO
extern double NEGZERO;
#endif

#endif /* __MINGW32__ */

double pow(double x, double y)
{
      double w, z, W, Wa, Wb, ya, yb, u;
/*    double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
      double aw, ay, wy;
      int e, i, nflg, iyflg, yoddint;

      if (y == 0.0)
            return (1.0);
#ifdef NANS
      if (isnan(x) || isnan(y))
      {
            _SET_ERRNO (EDOM);
            return (x + y);
      }
#endif
      if (y == 1.0)
            return (x);

#ifdef INFINITIES
      if (!isfinite(y) && (x == 1.0 || x == -1.0))
      {
            mtherr( "pow", DOMAIN );
#ifdef NANS
            return( NAN );
#else
            return( INFINITY );
#endif
      }
#endif

      if (x == 1.0)
            return (1.0);

      if (y >= MAXNUM)
      {
            _SET_ERRNO (ERANGE);
#ifdef INFINITIES
            if (x > 1.0)
                  return (INFINITY);
#else
            if (x > 1.0)
                  return (MAXNUM);
#endif
            if (x > 0.0 && x < 1.0)
                  return (0.0);
            if (x < -1.0)
            {
#ifdef INFINITIES
                  return (INFINITY);
#else
                  return (MAXNUM);
#endif
            }
            if (x > -1.0 && x < 0.0)
                  return (0.0);
      }
      if (y <= -MAXNUM)
      {
            _SET_ERRNO (ERANGE);
            if (x > 1.0)
            return (0.0);
#ifdef INFINITIES
            if (x > 0.0 && x < 1.0)
                  return (INFINITY);
#else
            if (x > 0.0 && x < 1.0)
                  return (MAXNUM);
#endif
            if (x < -1.0)
                  return (0.0);
#ifdef INFINITIES
            if (x > -1.0 && x < 0.0)
                  return (INFINITY);
#else
            if (x > -1.0 && x < 0.0)
                  return (MAXNUM);
#endif
      }
      if (x >= MAXNUM)
      {
#if INFINITIES
            if (y > 0.0)
                  return (INFINITY);
#else
            if (y > 0.0)
                  return (MAXNUM);
#endif
            return (0.0);
      }
      /* Set iyflg to 1 if y is an integer.  */
      iyflg = 0;
      w = floor(y);
      if (w == y)
            iyflg = 1;

/* Test for odd integer y.  */
      yoddint = 0;
      if (iyflg)
      {
            ya = fabs(y);
            ya = floor(0.5 * ya);
            yb = 0.5 * fabs(w);
            if (ya != yb)
                  yoddint = 1;
      }

      if (x <= -MAXNUM)
      {
            if (y > 0.0)
            {
#ifdef INFINITIES
                  if (yoddint)
                        return (-INFINITY);
                  return (INFINITY);
#else
                  if (yoddint)
                        return (-MAXNUM);
                  return (MAXNUM);
#endif
            }
            if (y < 0.0)
            {
#ifdef MINUSZERO
                  if (yoddint)
                        return (NEGZERO);
#endif
                  return (0.0);
            }
      }

      nflg = 0;   /* flag = 1 if x<0 raised to integer power */
      if (x <= 0.0)
      {
            if (x == 0.0)
            {
                  if (y < 0.0)
                  {
#ifdef MINUSZERO
                        if (signbit(x) && yoddint)
                              return (-INFINITY);
#endif
#ifdef INFINITIES
                        return (INFINITY);
#else
                        return (MAXNUM);
#endif
                  }
                  if (y > 0.0)
                  {
#ifdef MINUSZERO
                        if (signbit(x) && yoddint)
                              return (NEGZERO);
#endif
                        return (0.0);
                  }
                  return (1.0);
            }
            else
            {
                  if (iyflg == 0)
                  { /* noninteger power of negative number */
                        mtherr(fname, DOMAIN);
                        _SET_ERRNO (EDOM);
#ifdef NANS
                        return (NAN);
#else
                        return (0.0L);
#endif
                  }
                  nflg = 1;
            }
      }

      /* Integer power of an integer.  */

      if (iyflg)
      {
            i = w;
            w = floor(x);
            if( (w == x) && (fabs(y) < 32768.0) )
            {
                  w = __powi(x, (int) y);
                  return (w);
            }
      }

      if (nflg)
            x = fabs(x);

      /* For results close to 1, use a series expansion.  */
      w = x - 1.0;
      aw = fabs(w);
      ay = fabs(y);
      wy = w * y;
      ya = fabs(wy);
      if ((aw <= 1.0e-3 && ay <= 1.0)
         || (ya <= 1.0e-3 && ay >= 1.0))
      {
            z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
                  + 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
            goto done;
      }
      /* These are probably too much trouble.  */
#if 0
      w = y * log(x);
      if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
      {
            z = ((((((w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
            goto done;
      }

      if (ya <= 1.0e-3 && aw <= 1.0e-4)
      {
            z = (((((wy*1./720.
                            + (-w*1./48. + 1./120.) )*wy
                           + ((w*17./144. - 1./12.)*w + 1./24.) )*wy
                          + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
                         + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
                        + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
                        + wy + 1.0;
            goto done;
      }
#endif

      /* separate significand from exponent */
      x = frexp(x, &e);

#if 0
      /* For debugging, check for gross overflow. */
      if ((e * y)  > (MEXP + 1024))
            goto overflow;
#endif

      /* Find significand of x in antilog table A[]. */
      i = 1;
      if (x <= douba(9))
            i = 9;
      if (x <= douba(i + 4))
            i += 4;
      if (x <= douba(i + 2))
            i += 2;
      if (x >= douba(1))
            i = -1;
      i += 1;

      /* Find (x - A[i])/A[i]
       * in order to compute log(x/A[i]):
       *
       * log(x) = log( a x/a ) = log(a) + log(x/a)
       *
       * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
       */
      x -= douba(i);
      x -= doubb(i/2);
      x /= douba(i);

      /* rational approximation for log(1+v):
       *
       * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
       */
      z = x*x;
      w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
      w = w - ldexp( z, -1 );   /*  w - 0.5 * z  */

      /* Convert to base 2 logarithm:
       * multiply by log2(e)
       */
      w = w + LOG2EA * w;
      /* Note x was not yet added in
       * to above rational approximation,
       * so do it now, while multiplying
       * by log2(e).
       */
      z = w + LOG2EA * x;
      z = z + x;

      /* Compute exponent term of the base 2 logarithm. */
      w = -i;
      w = ldexp( w, -4 );     /* divide by 16 */
      w += e;
      /* Now base 2 log of x is w + z. */

      /* Multiply base 2 log by y, in extended precision. */

      /* separate y into large part ya
       * and small part yb less than 1/16
       */
      ya = reduc(y);
      yb = y - ya;

      F = z * y  +  w * yb;
      Fa = reduc(F);
      Fb = F - Fa;

      G = Fa + w * ya;
      Ga = reduc(G);
      Gb = G - Ga;

      H = Fb + Gb;
      Ha = reduc(H);
      w = ldexp(Ga + Ha, 4);

      /* Test the power of 2 for overflow */
      if (w > MEXP)
      {
#ifndef INFINITIES
            mtherr(fname, OVERFLOW);
#endif
#ifdef INFINITIES
            if (nflg && yoddint)
                  return (-INFINITY);
            return (INFINITY);
#else
            if (nflg && yoddint)
                  return (-MAXNUM);
            return (MAXNUM);
#endif
      }

      if (w < (MNEXP - 1))
      {
#ifndef DENORMAL
            mtherr(fname, UNDERFLOW);
#endif
#ifdef MINUSZERO
            if (nflg && yoddint)
                  return (NEGZERO);
#endif
            return (0.0);
      }

      e = w;
      Hb = H - Ha;

      if (Hb > 0.0)
      {
            e += 1;
            Hb -= 0.0625;
      }

      /* Now the product y * log2(x)  =  Hb + e/16.0.
       *
       * Compute base 2 exponential of Hb,
       * where -0.0625 <= Hb <= 0.
       */
      z = Hb * polevl(Hb, R, 6);  /*    z  =  2**Hb - 1    */

      /* Express e/16 as an integer plus a negative number of 16ths.
       * Find lookup table entry for the fractional power of 2.
       */
      if (e < 0)
            i = 0;
      else
            i = 1;
      i = e/16 + i;
      e = 16*i - e;
      w = douba(e);
      z = w + w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
      z = ldexp(z, i);  /* multiply by integer power of 2 */

done:
      /* Negate if odd integer power of negative number */
      if (nflg && yoddint)
      {
#ifdef MINUSZERO
            if (z == 0.0)
                  z = NEGZERO;
            else
#endif
                  z = -z;
      }
      return (z);
}


/* Find a multiple of 1/16 that is within 1/16 of x. */
static double reduc(double x)
{
      double t;

      t = ldexp(x, 4);
      t = floor(t);
      t = ldexp(t, -4);
      return (t);
}


Generated by  Doxygen 1.6.0   Back to index