home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ieee / elog.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.3 KB  |  91 lines
  1. /*                        xlog.c    */
  2. /* natural logarithm */
  3. /* by Stephen L. Moshier. */
  4.  
  5. #include "mconf.h"
  6. #include "ehead.h"
  7.  
  8.  
  9.  
  10. void elog( x, y )
  11. unsigned short *x, *y;
  12. {
  13. unsigned short xx[NE], z[NE], a[NE], b[NE], t[NE], qj[NE];
  14. long ex;
  15.  
  16.  
  17. if( x[NE-1] & (unsigned short )0x8000 )
  18.     {
  19.     eclear(y);
  20.     mtherr( "elog", DOMAIN );
  21.     return;
  22.     }
  23. if( ecmp( x, ezero ) == 0 )
  24.     {
  25.     einfin( y );
  26.     eneg(y);
  27.     mtherr( "elog", SING );
  28.     return;
  29.     }
  30. if( ecmp( x, eone ) == 0 )
  31.     {
  32.     eclear( y );
  33.     return;
  34.     }
  35.  
  36. /* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */
  37. efrexp( x, &ex, xx );
  38. /*
  39. emov(x, xx );
  40. ex = xx[NX-1] & 0x7fff;
  41. ex -= 0x3ffe;
  42. xx[NX-1] = 0x3ffe;
  43. */
  44.  
  45. /* Adjust range to 1/sqrt(2), sqrt(2) */
  46. esqrt2[NE-1] -= 1;
  47. if( ecmp( xx, esqrt2 ) < 0 )
  48.     {
  49.     ex -= 1;
  50.     emul( xx, etwo, xx );
  51.     }
  52. esqrt2[NE-1] += 1;
  53.  
  54. esub( eone, xx, a );
  55. if( a[NE-1] == 0 )
  56.     {
  57.     eclear( y );
  58.     goto logdon;
  59.     }
  60. eadd( eone, xx, b );
  61. ediv( b, a, y );    /* store (x-1)/(x+1) in y */
  62.  
  63. emul( y, y, z );
  64.  
  65. emov( eone, a );
  66. emov( eone, b );
  67. emov( eone, qj );
  68. do
  69.     {
  70.     eadd( etwo, qj, qj );    /* 2 * i + 1        */
  71.     emul( z, a, a );
  72.     ediv( qj, a, t );
  73.     eadd( t, b, b );
  74.     }
  75. while( ((b[NE-1] & 0x7fff) - (t[NE-1] & 0x7fff)) < NBITS );
  76.  
  77.  
  78. emul( b, y, y );
  79. emul( y, etwo, y );
  80.  
  81. logdon:
  82.  
  83. /* now add log of 2**ex */
  84. if( ex != 0 )
  85.     {
  86.     ltoe( &ex, b );
  87.     emul( elog2, b, b );
  88.     eadd( b, y, y );
  89.     }
  90. }
  91.