| OLD | NEW | 
|    1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.c */ |    1 /* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.c */ | 
|    2 /* |    2 /* | 
|    3  * ==================================================== |    3  * ==================================================== | 
|    4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |    4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 
|    5  * |    5  * | 
|    6  * Developed at SunPro, a Sun Microsystems, Inc. business. |    6  * Developed at SunPro, a Sun Microsystems, Inc. business. | 
|    7  * Permission to use, copy, modify, and distribute this |    7  * Permission to use, copy, modify, and distribute this | 
|    8  * software is freely granted, provided that this notice |    8  * software is freely granted, provided that this notice | 
|    9  * is preserved. |    9  * is preserved. | 
|   10  * ==================================================== |   10  * ==================================================== | 
| (...skipping 37 matching lines...) Expand 10 before | Expand all | Expand 10 after  Loading... | 
|   48  * |   48  * | 
|   49  *              u = 1+x; |   49  *              u = 1+x; | 
|   50  *              if(u==1.0) return x ; else |   50  *              if(u==1.0) return x ; else | 
|   51  *                         return log(u)*(x/(u-1.0)); |   51  *                         return log(u)*(x/(u-1.0)); | 
|   52  * |   52  * | 
|   53  *       See HP-15C Advanced Functions Handbook, p.193. |   53  *       See HP-15C Advanced Functions Handbook, p.193. | 
|   54  */ |   54  */ | 
|   55  |   55  | 
|   56 #include "libm.h" |   56 #include "libm.h" | 
|   57  |   57  | 
|   58 static const double |   58 static const double ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ | 
|   59 ln2_hi = 6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */ |   59     ln2_lo = 1.90821492927058770002e-10,                 /* 3dea39ef 35793c76 */ | 
|   60 ln2_lo = 1.90821492927058770002e-10,  /* 3dea39ef 35793c76 */ |   60     Lg1 = 6.666666666666735130e-01,                      /* 3FE55555 55555593 */ | 
|   61 Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */ |   61     Lg2 = 3.999999999940941908e-01,                      /* 3FD99999 9997FA04 */ | 
|   62 Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */ |   62     Lg3 = 2.857142874366239149e-01,                      /* 3FD24924 94229359 */ | 
|   63 Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */ |   63     Lg4 = 2.222219843214978396e-01,                      /* 3FCC71C5 1D8E78AF */ | 
|   64 Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */ |   64     Lg5 = 1.818357216161805012e-01,                      /* 3FC74664 96CB03DE */ | 
|   65 Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */ |   65     Lg6 = 1.531383769920937332e-01,                      /* 3FC39A09 D078C69F */ | 
|   66 Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */ |   66     Lg7 = 1.479819860511658591e-01;                      /* 3FC2F112 DF3E5244 */ | 
|   67 Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */ |  | 
|   68  |   67  | 
|   69 double log1p(double x) |   68 double log1p(double x) { | 
|   70 { |   69   union { | 
|   71 »       union {double f; uint64_t i;} u = {x}; |   70     double f; | 
|   72 »       double_t hfsq,f,c,s,z,R,w,t1,t2,dk; |   71     uint64_t i; | 
|   73 »       uint32_t hx,hu; |   72   } u = {x}; | 
|   74 »       int k; |   73   double_t hfsq, f, c, s, z, R, w, t1, t2, dk; | 
 |   74   uint32_t hx, hu; | 
 |   75   int k; | 
|   75  |   76  | 
|   76 »       hx = u.i>>32; |   77   hx = u.i >> 32; | 
|   77 »       k = 1; |   78   k = 1; | 
|   78 »       if (hx < 0x3fda827a || hx>>31) {  /* 1+x < sqrt(2)+ */ |   79   if (hx < 0x3fda827a || hx >> 31) { /* 1+x < sqrt(2)+ */ | 
|   79 »       »       if (hx >= 0xbff00000) {  /* x <= -1.0 */ |   80     if (hx >= 0xbff00000) {          /* x <= -1.0 */ | 
|   80 »       »       »       if (x == -1) |   81       if (x == -1) | 
|   81 »       »       »       »       return x/0.0; /* log1p(-1) = -inf */ |   82         return x / 0.0;     /* log1p(-1) = -inf */ | 
|   82 »       »       »       return (x-x)/0.0;     /* log1p(x<-1) = NaN */ |   83       return (x - x) / 0.0; /* log1p(x<-1) = NaN */ | 
|   83 »       »       } |   84     } | 
|   84 »       »       if (hx<<1 < 0x3ca00000<<1) {  /* |x| < 2**-53 */ |   85     if (hx << 1 < 0x3ca00000 << 1) { /* |x| < 2**-53 */ | 
|   85 »       »       »       /* underflow if subnormal */ |   86       /* underflow if subnormal */ | 
|   86 »       »       »       if ((hx&0x7ff00000) == 0) |   87       if ((hx & 0x7ff00000) == 0) | 
|   87 »       »       »       »       FORCE_EVAL((float)x); |   88         FORCE_EVAL((float)x); | 
|   88 »       »       »       return x; |   89       return x; | 
|   89 »       »       } |   90     } | 
|   90 »       »       if (hx <= 0xbfd2bec4) {  /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ |   91     if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ | 
|   91 »       »       »       k = 0; |   92       k = 0; | 
|   92 »       »       »       c = 0; |   93       c = 0; | 
|   93 »       »       »       f = x; |   94       f = x; | 
|   94 »       »       } |   95     } | 
|   95 »       } else if (hx >= 0x7ff00000) |   96   } else if (hx >= 0x7ff00000) | 
|   96 »       »       return x; |   97     return x; | 
|   97 »       if (k) { |   98   if (k) { | 
|   98 »       »       u.f = 1 + x; |   99     u.f = 1 + x; | 
|   99 »       »       hu = u.i>>32; |  100     hu = u.i >> 32; | 
|  100 »       »       hu += 0x3ff00000 - 0x3fe6a09e; |  101     hu += 0x3ff00000 - 0x3fe6a09e; | 
|  101 »       »       k = (int)(hu>>20) - 0x3ff; |  102     k = (int)(hu >> 20) - 0x3ff; | 
|  102 »       »       /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ |  103     /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */ | 
|  103 »       »       if (k < 54) { |  104     if (k < 54) { | 
|  104 »       »       »       c = k >= 2 ? 1-(u.f-x) : x-(u.f-1); |  105       c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1); | 
|  105 »       »       »       c /= u.f; |  106       c /= u.f; | 
|  106 »       »       } else |  107     } else | 
|  107 »       »       »       c = 0; |  108       c = 0; | 
|  108 »       »       /* reduce u into [sqrt(2)/2, sqrt(2)] */ |  109     /* reduce u into [sqrt(2)/2, sqrt(2)] */ | 
|  109 »       »       hu = (hu&0x000fffff) + 0x3fe6a09e; |  110     hu = (hu & 0x000fffff) + 0x3fe6a09e; | 
|  110 »       »       u.i = (uint64_t)hu<<32 | (u.i&0xffffffff); |  111     u.i = (uint64_t)hu << 32 | (u.i & 0xffffffff); | 
|  111 »       »       f = u.f - 1; |  112     f = u.f - 1; | 
|  112 »       } |  113   } | 
|  113 »       hfsq = 0.5*f*f; |  114   hfsq = 0.5 * f * f; | 
|  114 »       s = f/(2.0+f); |  115   s = f / (2.0 + f); | 
|  115 »       z = s*s; |  116   z = s * s; | 
|  116 »       w = z*z; |  117   w = z * z; | 
|  117 »       t1 = w*(Lg2+w*(Lg4+w*Lg6)); |  118   t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); | 
|  118 »       t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); |  119   t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); | 
|  119 »       R = t2 + t1; |  120   R = t2 + t1; | 
|  120 »       dk = k; |  121   dk = k; | 
|  121 »       return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi; |  122   return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi; | 
|  122 } |  123 } | 
| OLD | NEW |