libm/float/powif.c

   1 /*                                                      powif.c
   2  *
   3  *      Real raised to integer power
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, powif();
  10  * int n;
  11  *
  12  * y = powif( x, n );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns argument x raised to the nth power.
  19  * The routine efficiently decomposes n as a sum of powers of
  20  * two. The desired power is a product of two-to-the-kth
  21  * powers of x.  Thus to compute the 32767 power of x requires
  22  * 28 multiplications instead of 32767 multiplications.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *
  29  *                      Relative error:
  30  * arithmetic   x domain   n domain  # trials      peak         rms
  31  *    IEEE      .04,26     -26,26    100000       1.1e-6      2.0e-7
  32  *    IEEE        1,2      -128,128  100000       1.1e-5      1.0e-6
  33  *
  34  * Returns MAXNUMF on overflow, zero on underflow.
  35  *
  36  */
  37 \f
  38 /*                                                      powi.c  */
  39
  40 /*
  41 Cephes Math Library Release 2.2:  June, 1992
  42 Copyright 1984, 1987, 1989 by Stephen L. Moshier
  43 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  44 */
  45
  46 #include <math.h>
  47 extern float MAXNUMF, MAXLOGF, MINLOGF, LOGE2F;
  48
  49 float frexpf( float, int * );
  50
  51 float powif( float x, int nn )
  52 {
  53 int n, e, sign, asign, lx;
  54 float w, y, s;
  55
  56 if( x == 0.0 )
  57         {
  58         if( nn == 0 )
  59                 return( 1.0 );
  60         else if( nn < 0 )
  61                 return( MAXNUMF );
  62         else
  63                 return( 0.0 );
  64         }
  65
  66 if( nn == 0 )
  67         return( 1.0 );
  68
  69
  70 if( x < 0.0 )
  71         {
  72         asign = -1;
  73         x = -x;
  74         }
  75 else
  76         asign = 0;
  77
  78
  79 if( nn < 0 )
  80         {
  81         sign = -1;
  82         n = -nn;
  83 /*
  84         x = 1.0/x;
  85 */
  86         }
  87 else
  88         {
  89         sign = 0;
  90         n = nn;
  91         }
  92
  93 /* Overflow detection */
  94
  95 /* Calculate approximate logarithm of answer */
  96 s = frexpf( x, &lx );
  97 e = (lx - 1)*n;
  98 if( (e == 0) || (e > 64) || (e < -64) )
  99         {
 100         s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
 101         s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
 102         }
 103 else
 104         {
 105         s = LOGE2F * e;
 106         }
 107
 108 if( s > MAXLOGF )
 109         {
 110         mtherr( "powi", OVERFLOW );
 111         y = MAXNUMF;
 112         goto done;
 113         }
 114
 115 if( s < MINLOGF )
 116         return(0.0);
 117
 118 /* Handle tiny denormal answer, but with less accuracy
 119  * since roundoff error in 1.0/x will be amplified.
 120  * The precise demarcation should be the gradual underflow threshold.
 121  */
 122 if( s < (-MAXLOGF+2.0) )
 123         {
 124         x = 1.0/x;
 125         sign = 0;
 126         }
 127
 128 /* First bit of the power */
 129 if( n & 1 )
 130         y = x;
 131
 132 else
 133         {
 134         y = 1.0;
 135         asign = 0;
 136         }
 137
 138 w = x;
 139 n >>= 1;
 140 while( n )
 141         {
 142         w = w * w;      /* arg to the 2-to-the-kth power */
 143         if( n & 1 )     /* if that bit is set, then include in product */
 144                 y *= w;
 145         n >>= 1;
 146         }
 147
 148
 149 done:
 150
 151 if( asign )
 152         y = -y; /* odd power of negative number */
 153 if( sign )
 154         y = 1.0/y;
 155 return(y);
 156 }