Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /*                                                      powi.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 = powi( 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.  *    DEC       .04,26     -26,26    100000       2.7e-16     4.3e-17
  32.  *    IEEE      .04,26     -26,26     50000       2.0e-15     3.8e-16
  33.  *    IEEE        1,2    -1022,1023   50000       8.6e-14     1.6e-14
  34.  *
  35.  * Returns MAXNUM on overflow, zero on underflow.
  36.  *
  37.  */
  39. /*                                                      powi.c  */
  40.  
  41. /*
  42. Cephes Math Library Release 2.8:  June, 2000
  43. Copyright 1984, 1995, 2000 by Stephen L. Moshier
  44. */
  45.  
  46. /*
  47. Modified for float from powi.c and adapted to mingw
  48. 2002-10-01 Danny Smith <dannysmith@users.sourceforge.net>
  49. */
  50.  
  51. #ifdef __MINGW32__
  52. #include "cephes_mconf.h"
  53. #else
  54. #include "mconf.h"
  55. #ifdef ANSIPROT
  56. extern float logf ( float );
  57. extern float frexpf ( float, int * );
  58. extern int signbitf ( float );
  59. #else
  60. float logf(), frexpf();
  61. int signbitf();
  62. #endif
  63. extern float NEGZEROF, INFINITYF, MAXNUMF, MAXLOGF, MINLOGF, LOGE2F;
  64. #endif /* __MINGW32__ */
  65.  
  66. #ifndef _SET_ERRNO
  67. #define _SET_ERRNO(x)
  68. #endif
  69.  
  70. float __powif( float x, int nn )
  71. {
  72. int n, e, sign, asign, lx;
  73. float w, y, s;
  74.  
  75. /* See pow.c for these tests.  */
  76. if( x == 0.0F )
  77.         {
  78.         if( nn == 0 )
  79.                 return( 1.0F );
  80.         else if( nn < 0 )
  81.             return( INFINITYF );
  82.         else
  83.           {
  84.             if( nn & 1 )
  85.               return( x );
  86.             else
  87.               return( 0.0 );
  88.           }
  89.         }
  90.  
  91. if( nn == 0 )
  92.         return( 1.0 );
  93.  
  94. if( nn == -1 )
  95.         return( 1.0/x );
  96.  
  97. if( x < 0.0 )
  98.         {
  99.         asign = -1;
  100.         x = -x;
  101.         }
  102. else
  103.         asign = 0;
  104.  
  105.  
  106. if( nn < 0 )
  107.         {
  108.         sign = -1;
  109.         n = -nn;
  110.         }
  111. else
  112.         {
  113.         sign = 1;
  114.         n = nn;
  115.         }
  116.  
  117. /* Even power will be positive. */
  118. if( (n & 1) == 0 )
  119.         asign = 0;
  120.  
  121. /* Overflow detection */
  122.  
  123. /* Calculate approximate logarithm of answer */
  124. s = frexpf( x, &lx );
  125. e = (lx - 1)*n;
  126. if( (e == 0) || (e > 64) || (e < -64) )
  127.         {
  128.         s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
  129.         s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
  130.         }
  131. else
  132.         {
  133.         s = LOGE2F * e;
  134.         }
  135.  
  136. if( s > MAXLOGF )
  137.         {
  138.         mtherr( "__powif", OVERFLOW );
  139.         _SET_ERRNO(ERANGE);
  140.         y = INFINITYF;
  141.         goto done;
  142.         }
  143.  
  144. #if DENORMAL
  145. if( s < MINLOGF )
  146.         {
  147.         y = 0.0;
  148.         goto done;
  149.         }
  150.  
  151. /* Handle tiny denormal answer, but with less accuracy
  152.  * since roundoff error in 1.0/x will be amplified.
  153.  * The precise demarcation should be the gradual underflow threshold.
  154.  */
  155. if( (s < (-MAXLOGF+2.0)) && (sign < 0) )
  156.         {
  157.         x = 1.0/x;
  158.         sign = -sign;
  159.         }
  160. #else
  161. /* do not produce denormal answer */
  162. if( s < -MAXLOGF )
  163.         return(0.0);
  164. #endif
  165.  
  166.  
  167. /* First bit of the power */
  168. if( n & 1 )
  169.         y = x;
  170.                
  171. else
  172.         y = 1.0;
  173.  
  174. w = x;
  175. n >>= 1;
  176. while( n )
  177.         {
  178.         w = w * w;      /* arg to the 2-to-the-kth power */
  179.         if( n & 1 )     /* if that bit is set, then include in product */
  180.                 y *= w;
  181.         n >>= 1;
  182.         }
  183.  
  184. if( sign < 0 )
  185.         y = 1.0/y;
  186.  
  187. done:
  188.  
  189. if( asign )
  190.         {
  191.         /* odd power of negative number */
  192.         if( y == 0.0 )
  193.                 y = NEGZEROF;
  194.         else
  195.                 y = -y;
  196.         }
  197. return(y);
  198. }
  199.