Line data    Source code

       1             : 
       2             : /* @(#)e_log.c 1.3 95/01/18 */
       3             : /*
       4             :  * ====================================================
       5             :  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       6             :  *
       7             :  * Developed at SunSoft, a Sun Microsystems, Inc. business.
       8             :  * Permission to use, copy, modify, and distribute this
       9             :  * software is freely granted, provided that this notice 
      10             :  * is preserved.
      11             :  * ====================================================
      12             :  */
      13             : 
      14             : //#include <sys/cdefs.h>
      15             : //__FBSDID("$FreeBSD$");
      16             : 
      17             : /*
      18             :  * k_log1p(f):
      19             :  * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)].
      20             :  *
      21             :  * The following describes the overall strategy for computing
      22             :  * logarithms in base e.  The argument reduction and adding the final
      23             :  * term of the polynomial are done by the caller for increased accuracy
      24             :  * when different bases are used.
      25             :  *
      26             :  * Method :                  
      27             :  *   1. Argument Reduction: find k and f such that 
      28             :  *                      x = 2^k * (1+f), 
      29             :  *         where  sqrt(2)/2 < 1+f < sqrt(2) .
      30             :  *
      31             :  *   2. Approximation of log(1+f).
      32             :  *      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
      33             :  *               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
      34             :  *               = 2s + s*R
      35             :  *      We use a special Reme algorithm on [0,0.1716] to generate 
      36             :  *      a polynomial of degree 14 to approximate R The maximum error 
      37             :  *      of this polynomial approximation is bounded by 2**-58.45. In
      38             :  *      other words,
      39             :  *                      2      4      6      8      10      12      14
      40             :  *          R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
      41             :  *      (the values of Lg1 to Lg7 are listed in the program)
      42             :  *      and
      43             :  *          |      2          14          |     -58.45
      44             :  *          | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
      45             :  *          |                             |
      46             :  *      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
      47             :  *      In order to guarantee error in log below 1ulp, we compute log
      48             :  *      by
      49             :  *              log(1+f) = f - s*(f - R)        (if f is not too large)
      50             :  *              log(1+f) = f - (hfsq - s*(hfsq+R)).     (better accuracy)
      51             :  *      
      52             :  *      3. Finally,  log(x) = k*ln2 + log(1+f).  
      53             :  *                          = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
      54             :  *         Here ln2 is split into two floating point number: 
      55             :  *                      ln2_hi + ln2_lo,
      56             :  *         where n*ln2_hi is always exact for |n| < 2000.
      57             :  *
      58             :  * Special cases:
      59             :  *      log(x) is NaN with signal if x < 0 (including -INF) ; 
      60             :  *      log(+INF) is +INF; log(0) is -INF with signal;
      61             :  *      log(NaN) is that NaN with no signal.
      62             :  *
      63             :  * Accuracy:
      64             :  *      according to an error analysis, the error is always less than
      65             :  *      1 ulp (unit in the last place).
      66             :  *
      67             :  * Constants:
      68             :  * The hexadecimal values are the intended ones for the following 
      69             :  * constants. The decimal values may be used, provided that the 
      70             :  * compiler will convert from decimal to binary accurately enough 
      71             :  * to produce the hexadecimal values shown.
      72             :  */
      73             : 
      74             : static const double
      75             : Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
      76             : Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
      77             : Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
      78             : Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
      79             : Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
      80             : Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
      81             : Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
      82             : 
      83             : /*
      84             :  * We always inline k_log1p(), since doing so produces a
      85             :  * substantial performance improvement (~40% on amd64).
      86             :  */
      87             : static inline double
      88           0 : k_log1p(double f)
      89             : {
      90             :         double hfsq,s,z,R,w,t1,t2;
      91             : 
      92           0 :         s = f/(2.0+f);
      93           0 :         z = s*s;
      94           0 :         w = z*z;
      95           0 :         t1= w*(Lg2+w*(Lg4+w*Lg6));
      96           0 :         t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
      97           0 :         R = t2+t1;
      98           0 :         hfsq=0.5*f*f;
      99           0 :         return s*(hfsq+R);
     100             : }