Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      __powil.c
  2.  *
  3.  *      Real raised to integer power, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, __powil();
  10.  * int n;
  11.  *
  12.  * y = __powil( 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     .001,1000  -1022,1023  50000       4.3e-17     7.8e-18
  32.  *    IEEE        1,2     -1022,1023  20000       3.9e-17     7.6e-18
  33.  *    IEEE     .99,1.01     0,8700    10000       3.6e-16     7.2e-17
  34.  *
  35.  * Returns INFINITY on overflow, zero on underflow.
  36.  *
  37.  */
  38.  
  39. /*                                                      __powil.c       */
  40.  
  41. /*
  42. Cephes Math Library Release 2.2:  December, 1990
  43. Copyright 1984, 1990 by Stephen L. Moshier
  44. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  45. */
  46.  
  47. /*
  48. Modified for mingw
  49. 2002-07-22 Danny Smith <dannysmith@users.sourceforge.net>
  50. */
  51.  
  52. #ifdef __MINGW32__
  53. #include "cephes_mconf.h"
  54. #else
  55. #include "mconf.h"
  56. extern long double MAXNUML, MAXLOGL, MINLOGL;
  57. extern long double LOGE2L;
  58. #ifdef ANSIPROT
  59. extern long double frexpl ( long double, int * );
  60. #else
  61. long double frexpl();
  62. #endif
  63. #endif /* __MINGW32__ */
  64.  
  65. #ifndef _SET_ERRNO
  66. #define _SET_ERRNO(x)
  67. #endif
  68.  
  69. long double __powil( x, nn )
  70. long double x;
  71. int nn;
  72. {
  73. long double w, y;
  74. long double s;
  75. int n, e, sign, asign, lx;
  76.  
  77. if( x == 0.0L )
  78.         {
  79.         if( nn == 0 )
  80.                 return( 1.0L );
  81.         else if( nn < 0 )
  82.                 return( INFINITYL );
  83.         else
  84.                 return( 0.0L );
  85.         }
  86.  
  87. if( nn == 0 )
  88.         return( 1.0L );
  89.  
  90.  
  91. if( x < 0.0L )
  92.         {
  93.         asign = -1;
  94.         x = -x;
  95.         }
  96. else
  97.         asign = 0;
  98.  
  99.  
  100. if( nn < 0 )
  101.         {
  102.         sign = -1;
  103.         n = -nn;
  104.         }
  105. else
  106.         {
  107.         sign = 1;
  108.         n = nn;
  109.         }
  110.  
  111. /* Overflow detection */
  112.  
  113. /* Calculate approximate logarithm of answer */
  114. s = x;
  115. s = frexpl( s, &lx );
  116. e = (lx - 1)*n;
  117. if( (e == 0) || (e > 64) || (e < -64) )
  118.         {
  119.         s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);
  120.         s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
  121.         }
  122. else
  123.         {
  124.         s = LOGE2L * e;
  125.         }
  126.  
  127. if( s > MAXLOGL )
  128.         {
  129.         mtherr( "__powil", OVERFLOW );
  130.         _SET_ERRNO(ERANGE);
  131.         y = INFINITYL;
  132.         goto done;
  133.         }
  134.  
  135. if( s < MINLOGL )
  136.         {
  137.         mtherr( "__powil", UNDERFLOW );
  138.         _SET_ERRNO(ERANGE);
  139.         return(0.0L);
  140.         }
  141. /* Handle tiny denormal answer, but with less accuracy
  142.  * since roundoff error in 1.0/x will be amplified.
  143.  * The precise demarcation should be the gradual underflow threshold.
  144.  */
  145. if( s < (-MAXLOGL+2.0L) )
  146.         {
  147.         x = 1.0L/x;
  148.         sign = -sign;
  149.         }
  150.  
  151. /* First bit of the power */
  152. if( n & 1 )
  153.         y = x;
  154.                
  155. else
  156.         {
  157.         y = 1.0L;
  158.         asign = 0;
  159.         }
  160.  
  161. w = x;
  162. n >>= 1;
  163. while( n )
  164.         {
  165.         w = w * w;      /* arg to the 2-to-the-kth power */
  166.         if( n & 1 )     /* if that bit is set, then include in product */
  167.                 y *= w;
  168.         n >>= 1;
  169.         }
  170.  
  171.  
  172. done:
  173.  
  174. if( asign )
  175.         y = -y; /* odd power of negative number */
  176. if( sign < 0 )
  177.         y = 1.0L/y;
  178. return(y);
  179. }
  180.