home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / POWERI.C < prev    next >
Encoding:
Text File  |  1987-01-03  |  1.9 KB  |  124 lines
  1. /*                            powi.c
  2.  *
  3.  *    Real raised to integer power
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, powi();
  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.  * Approximately 1e-16 relative.
  29.  *
  30.  * Returns MAXNUM or zero on overflow, which is detected
  31.  * by integer estimation of the exponent of the result.
  32.  *
  33.  */
  34.  
  35. /*                            powi.c    */
  36.  
  37. /* Cephes Math Library Release 1.1:  March, 1985
  38.  * Copyright 1985 by Stephen L. Moshier
  39.  * Contributed to BIX for personal, noncommercial use only.
  40.  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
  41.  
  42. extern double MAXNUM, MAXLOG;
  43.  
  44. double powi( x, n )
  45. double x;
  46. int n;
  47. {
  48. int e, sign, asign;
  49. short *p;
  50. double w, y;
  51. double log();
  52.  
  53. if( x == 0.0 )
  54.     {
  55.     if( n == 0 )
  56.         return( 1.0 );
  57.     else if( n < 0 )
  58.         return( MAXNUM );
  59.     else
  60.         return( 0.0 );
  61.     }
  62.  
  63. if( n == 0 )
  64.     return( 1.0 );
  65.  
  66.  
  67. if( x < 0.0 )
  68.     {
  69.     asign = -1;
  70.     x = -x;
  71.     }
  72. else
  73.     asign = 0;
  74.  
  75.  
  76. if( n < 0 )
  77.     {
  78.     sign = -1;
  79.     n = -n;
  80.     x = 1.0/x;
  81.     }
  82. else
  83.     sign = 0;
  84.  
  85. /*                            powi.c    */
  86.  
  87. /* First bit of the power */
  88. if( n & 1 )
  89.     y = x;
  90.         
  91. else
  92.     {
  93.     y = 1.0;
  94.     asign = 0;
  95.     }
  96.  
  97. w = n * log( x );    /* Overflow detection */
  98. if( w >= MAXLOG )
  99.     {
  100.     y = MAXNUM;
  101.     goto done;
  102.     }
  103. if( w < -MAXLOG )
  104.     return(0.0);
  105.  
  106. w = x;
  107. n >>= 1;
  108. while( n )
  109.     {
  110.     w = w * w;    /* arg to the 2-to-the-kth power */
  111.     if( n & 1 )    /* if that bit is set, then include in product */
  112.         y *= w;
  113.     n >>= 1;
  114.     }
  115.  
  116.  
  117. done:
  118.  
  119. if( asign )
  120.     return( -y ); /* odd power of negative number */
  121.  
  122. return(y);
  123. }
  124.