home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / qfloat / qlog.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.1 KB  |  81 lines
  1. /*                        qlog.c    */
  2. /* natural logarithm */
  3.  
  4. #include "qhead.h"
  5.  
  6. extern short qone[], qtwo[], qlog2[], qsqrt2[];
  7.  
  8. qlog( x, y )
  9. short *x, *y;
  10. {
  11. short xx[NQ], z[NQ], a[NQ], b[NQ], t[NQ], qj[NQ];
  12. long i, j, n, nsq;
  13. long ex;
  14.  
  15. if( x[0] < 0 )
  16.     {
  17.     qclear(y);
  18.     printf( "qlog domain error\n" );
  19.     return;
  20.     }
  21. if( x[1] == 0 )
  22.     {
  23.     qinfin( y );
  24.     y[0] = -1;
  25.     printf( "qlog singularity\n" );
  26.     return(0);
  27.     }
  28.  
  29. /* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */
  30. qmov(x, xx );
  31. ex = *(xx+1);
  32. ex -= 040000;
  33. xx[1] = 040000;
  34.  
  35. /* Adjust range to 1/sqrt(2), sqrt(2) */
  36. qsqrt2[1] -= 1;
  37. if( qcmp( xx, qsqrt2 ) < 0 )
  38.     {
  39.     ex -= 1;
  40.     xx[1] += 1;
  41.     }
  42. qsqrt2[1] += 1;
  43.  
  44. qadd( qone, xx, b );
  45. qsub( qone, xx, a );
  46. if( a[1] == 0 )
  47.     {
  48.     qclear(y);
  49.     goto bdone;
  50.     }
  51. qdiv( b, a, y );    /* store (x-1)/(x+1) in y */
  52.  
  53. qmul( y, y, z );
  54.  
  55. qmov( qone, a );
  56. qmov( qone, b );
  57. qmov( qone, qj );
  58. do
  59.     {
  60.     qadd( qtwo, qj, qj );    /* 2 * i + 1        */
  61.     qmul( z, a, a );
  62.     qdiv( qj, a, t );
  63.     qadd( t, b, b );
  64.     }
  65. while( (b[1] - t[1]) < NBITS );
  66.  
  67.  
  68. qmul( b, y, y );
  69. y[1] += 1;
  70.  
  71. bdone:
  72.  
  73. /* now add log of 2**ex */
  74. if( ex != 0 )
  75.     {
  76.     ltoq( &ex, b );
  77.     qmul( qlog2, b, b );
  78.     qadd( b, y, y );
  79.     }
  80. }
  81.