home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ldouble / powil.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.5 KB  |  162 lines
  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 MAXNUM 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. #include "mconf.h"
  48. extern long double MAXNUML, MAXLOGL, MINLOGL;
  49. extern double LOGE2;
  50.  
  51. long double powil( x, nn )
  52. long double x;
  53. int nn;
  54. {
  55. long double w, y;
  56. double s;
  57. int n, e, sign, asign, lx;
  58. double frexp();
  59.  
  60. if( x == 0.0L )
  61.     {
  62.     if( nn == 0 )
  63.         return( 1.0 );
  64.     else if( nn < 0 )
  65.         return( MAXNUML );
  66.     else
  67.         return( 0.0L );
  68.     }
  69.  
  70. if( nn == 0 )
  71.     return( 1.0L );
  72.  
  73.  
  74. if( x < 0.0L )
  75.     {
  76.     asign = -1;
  77.     x = -x;
  78.     }
  79. else
  80.     asign = 0;
  81.  
  82.  
  83. if( nn < 0 )
  84.     {
  85.     sign = -1;
  86.     n = -nn;
  87. /*
  88.     x = 1.0/x;
  89. */
  90.     }
  91. else
  92.     {
  93.     sign = 0;
  94.     n = nn;
  95.     }
  96.  
  97. /* Overflow detection */
  98.  
  99. /* Calculate approximate logarithm of answer */
  100. s = (double )x;
  101. s = frexp( s, &lx );
  102. e = (lx - 1)*n;
  103. if( (e == 0) || (e > 64) || (e < -64) )
  104.     {
  105.     s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
  106.     s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
  107.     }
  108. else
  109.     {
  110.     s = LOGE2 * e;
  111.     }
  112.  
  113. if( s > MAXLOGL )
  114.     {
  115.     mtherr( "powil", OVERFLOW );
  116.     y = MAXNUML;
  117.     goto done;
  118.     }
  119.  
  120. if( s < MINLOGL )
  121.     return(0.0);
  122.  
  123. /* Handle tiny denormal answer, but with less accuracy
  124.  * since roundoff error in 1.0/x will be amplified.
  125.  * The precise demarcation should be the gradual underflow threshold.
  126.  */
  127. if( s < (-MAXLOGL+2.0) )
  128.     {
  129.     x = 1.0/x;
  130.     sign = 0;
  131.     }
  132.  
  133. /* First bit of the power */
  134. if( n & 1 )
  135.     y = x;
  136.         
  137. else
  138.     {
  139.     y = 1.0;
  140.     asign = 0;
  141.     }
  142.  
  143. w = x;
  144. n >>= 1;
  145. while( n )
  146.     {
  147.     w = w * w;    /* arg to the 2-to-the-kth power */
  148.     if( n & 1 )    /* if that bit is set, then include in product */
  149.         y *= w;
  150.     n >>= 1;
  151.     }
  152.  
  153.  
  154. done:
  155.  
  156. if( asign )
  157.     y = -y; /* odd power of negative number */
  158. if( sign )
  159.     y = 1.0/y;
  160. return(y);
  161. }
  162.