home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ieee / epow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.1 KB  |  170 lines
  1. /*                        epow.c    */
  2. /*  power function: z = x**y */
  3. /*  by Stephen L. Moshier. */
  4.  
  5.  
  6. #include "ehead.h"
  7.  
  8. extern int rndprc;
  9. void epowi();
  10.  
  11. void epow( x, y, z )
  12. unsigned short *x, *y, *z;
  13. {
  14. unsigned short w[NE];
  15. int rndsav;
  16. long li;
  17.  
  18. efloor( y, w );
  19. if( ecmp(y,w) == 0 )
  20.     {
  21.     eifrac( y, &li, w );
  22.     if( li < 0 )
  23.         li = -li;
  24.     if( li < 0x7fffffff )
  25.         {
  26.         epowi( x, y, z );
  27.         return;
  28.         }
  29.     }
  30. /* z = exp( y * log(x) ) */
  31. rndsav = rndprc;
  32. rndprc = NBITS;
  33. elog( x, w );
  34. emul( y, w, w );
  35. eexp( w, z );
  36. rndprc = rndsav;
  37. emul( eone, z, z );
  38. }
  39.  
  40.  
  41. /* y is integer valued. */
  42.  
  43. void epowi( x, y, z )
  44. unsigned short x[], y[], z[];
  45. {
  46. unsigned short w[NE];
  47. long li, lx;
  48. unsigned long lu;
  49. int rndsav;
  50. unsigned short signx;
  51. /* unsigned short signy; */
  52.  
  53. rndsav = rndprc;
  54. eifrac( y, &li, w );
  55. if( li < 0 )
  56.     lx = -li;
  57. else
  58.     lx = li;
  59.  
  60. if( lx == 0x7fffffff )
  61.     {
  62.     epow( x, y, z );
  63.     goto done;
  64.     }
  65.  
  66. if( (x[NE-1] & (unsigned short )0x7fff) == 0.0 )
  67.     {
  68.     if( li == 0 )
  69.         {
  70.         emov( eone, z );
  71.         return;
  72.         }
  73.     else if( li < 0 )
  74.         {
  75.         einfin( z );
  76.         return;
  77.         }
  78.     else
  79.         {
  80.         eclear( z );
  81.         return;
  82.         }
  83.     }
  84.  
  85. if( li == 0L )
  86.     {
  87.     emov( eone, z );
  88.     return;
  89.     }
  90.  
  91. emov( x, w );
  92. signx = w[NE-1] & (unsigned short )0x8000;
  93. w[NE-1] &= (unsigned short )0x7fff;
  94.  
  95. /* Overflow detection */
  96. /*
  97. lx = li * (w[NE-1] - 0x3fff);
  98. if( lx > 16385L )
  99.     {
  100.     einfin( z );
  101.     mtherr( "epowi", OVERFLOW );
  102.     goto done;
  103.     }
  104. if( lx < -16450L )
  105.     {
  106.     eclear( z );
  107.     return;
  108.     }
  109. */
  110. rndprc = NBITS;
  111.  
  112. if( li < 0 )
  113.     {
  114.     lu = (unsigned int )( -li );
  115. /*    signy = 0xffff;*/
  116.     ediv( w, eone, w );
  117.     }
  118. else
  119.     {
  120.     lu = (unsigned int )li;
  121. /*    signy = 0;*/
  122.     }
  123.  
  124. /* First bit of the power */
  125. if( lu & 1 )
  126.     {
  127.     emov( w, z );
  128.     }    
  129. else
  130.     {
  131.     emov( eone, z );
  132.     signx = 0;
  133.     }
  134.  
  135.  
  136. lu >>= 1;
  137. while( lu != 0L )
  138.     {
  139.     emul( w, w, w );    /* arg to the 2-to-the-kth power */
  140.     if( lu & 1L )    /* if that bit is set, then include in product */
  141.         emul( w, z, z );
  142.     lu >>= 1;
  143.     }
  144.  
  145.  
  146. done:
  147.  
  148. if( signx )
  149.     eneg( z ); /* odd power of negative number */
  150.  
  151. /*
  152. if( signy )
  153.       {
  154.       if( ecmp( z, ezero ) != 0 )
  155.          {
  156.         ediv( z, eone, z );
  157.         }
  158.     else
  159.         {
  160.         einfin( z );
  161.         printf( "epowi OVERFLOW\n" );
  162.         }
  163.     }
  164. */
  165. rndprc = rndsav;
  166. emul( eone, z, z );
  167. }
  168.  
  169.  
  170.