home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / qfloat / qpow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.6 KB  |  153 lines
  1. /*                        qpow.c    */
  2. /*  power function: z = x**y */
  3.  
  4. #include "qhead.h"
  5. #include "mconf.h"
  6.  
  7. extern short qone[];
  8.  
  9.  
  10. qpow( x, y, z )
  11. short *x, *y, *z;
  12. {
  13. short w[NQ];
  14. long li;
  15.  
  16. qfloor( y, w );
  17. if( qcmp(y,w) == 0 )
  18.     {
  19.     qifrac( y, &li, w );
  20.     if( li < 0 )
  21.         li = -li;
  22.     if( li < 0x7fffffff )
  23.         {
  24.         qpowi( x, y, z );
  25.         return;
  26.         }
  27.     }
  28. /* z = exp( y * log(x) ) */
  29.  
  30. qlog( x, w );
  31. qmul( y, w, w );
  32. qexp( w, z );
  33. }
  34.  
  35.  
  36. /* y is integer valued. */
  37.  
  38. qpowi( x, y, z )
  39. short x[], y[], z[];
  40. {
  41. short w[NQ];
  42. long li, lx;
  43. int signx, signy;
  44.  
  45. qifrac( y, &li, w );
  46. if( li < 0 )
  47.     lx = -li;
  48. else
  49.     lx = li;
  50.  
  51. if( lx == 0x7fffffff )
  52.     {
  53.     qpow( x, y, z );
  54.     return;
  55.     }
  56.  
  57.  
  58. if( x[1] == 0.0 )
  59.     {
  60.     if( li == 0 )
  61.         {
  62.         qmov( qone, z );
  63.         return;
  64.         }
  65.     else if( li < 0 )
  66.         {
  67.         qinfin( z );
  68.         return;
  69.         }
  70.     else
  71.         {
  72.         qclear( z );
  73.         return;
  74.         }
  75.     }
  76.  
  77. if( li == 0L )
  78.     {
  79.     qmov( qone, z );
  80.     return;
  81.     }
  82.  
  83. qmov( x, w );
  84. signx = w[0];
  85. w[0] = 0;
  86.  
  87.  
  88. if( li < 0 )
  89.     {
  90.     li = -li;
  91.     signy = -1;
  92.     }
  93. else
  94.     signy = 0;
  95.  
  96.  
  97. /* First bit of the power */
  98. if( li & 1 )
  99.     {
  100.     qmov( w, z );
  101.     }    
  102. else
  103.     {
  104.     qmov( qone, z );
  105.     signx = 0;
  106.     }
  107.  
  108. /* Overflow detection */
  109.  
  110. lx = signy * li * (w[1] - 0x4001);
  111. if( lx > 16383L )
  112.     {
  113.     qinfin( z );
  114.     mtherr( "qpowi", OVERFLOW );
  115.     goto done;
  116.     }
  117. if( lx < -16383L )
  118.     {
  119.     qclear( z );
  120.     return;
  121.     }
  122.  
  123.  
  124. li >>= 1;
  125. while( li != 0L )
  126.     {
  127.     qmul( w, w, w );    /* arg to the 2-to-the-kth power */
  128.     if( li & 1L )    /* if that bit is set, then include in product */
  129.         qmul( w, z, z );
  130.     li >>= 1;
  131.     }
  132.  
  133.  
  134. done:
  135.  
  136. if( signx )
  137.     qneg( z ); /* odd power of negative number */
  138.  
  139. if( signy )
  140.     {
  141.     if( z[1] != 0 )
  142.         {
  143.         qdiv( z, qone, z );
  144.         }
  145.     else
  146.         {
  147.         qinfin( z );
  148.         mtherr( "qpowi", OVERFLOW );
  149.         }
  150.     }
  151. }
  152.  
  153.