Changeset 19520ba9 in mainline
- Timestamp:
- 2021-06-18T16:05:49Z (3 years ago)
- Children:
- f250c5a
- Parents:
- ec6081c
- git-author:
- Maurizio Lombardi <mlombard@…> (2018-11-25 16:43:14)
- git-committer:
- Maurizio Lombardi <mlombard@…> (2021-06-18 16:05:49)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/math/generic/log.c
rec6081c r19520ba9 69 69 #include <math.h> 70 70 #include <stdint.h> 71 72 #include "internal.h" 73 71 74 static const double 72 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 73 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 75 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 76 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 77 two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ 74 78 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 75 79 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ … … 80 84 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 81 85 86 static const double zero = 0.0; 87 static volatile double vzero = 0.0; 88 89 82 90 double log(double x) 83 91 { 84 union {double f; uint64_t i;} u = {x}; 85 double_t hfsq,f,s,z,R,w,t1,t2,dk; 86 uint32_t hx; 87 int k; 88 hx = u.i>>32; 89 k = 0; 90 if (hx < 0x00100000 || hx>>31) { 91 if (u.i<<1 == 0) 92 return -1/(x*x); /* log(+-0)=-inf */ 93 if (hx>>31) 94 return (x-x)/0.0; /* log(-#) = NaN */ 95 /* subnormal number, scale x up */ 96 k -= 54; 97 x *= 0x1p54; 98 u.f = x; 99 hx = u.i>>32; 100 } else if (hx >= 0x7ff00000) { 101 return x; 102 } else if (hx == 0x3ff00000 && u.i<<32 == 0) 103 return 0; 92 double hfsq,f,s,z,R,w,t1,t2,dk; 93 int32_t k,hx,i,j; 94 uint32_t lx; 104 95 105 /* reduce x into [sqrt(2)/2, sqrt(2)] */ 106 hx += 0x3ff00000 - 0x3fe6a09e; 107 k += (int)(hx>>20) - 0x3ff; 108 hx = (hx&0x000fffff) + 0x3fe6a09e; 109 u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); 110 x = u.f; 111 f = x - 1.0; 112 hfsq = 0.5*f*f; 113 s = f/(2.0+f); 96 EXTRACT_WORDS(hx,lx,x); 97 98 k=0; 99 if (hx < 0x00100000) { /* x < 2**-1022 */ 100 if (((hx&0x7fffffff)|lx)==0) 101 return -two54/vzero; /* log(+-0)=-inf */ 102 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 103 k -= 54; x *= two54; /* subnormal number, scale up x */ 104 GET_HIGH_WORD(hx,x); 105 } 106 if (hx >= 0x7ff00000) return x+x; 107 k += (hx>>20)-1023; 108 hx &= 0x000fffff; 109 i = (hx+0x95f64)&0x100000; 110 SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ 111 k += (i>>20); 112 f = x-1.0; 113 if((0x000fffff&(2+hx))<3) { /* -2**-20 <= f < 2**-20 */ 114 if(f==zero) { 115 if(k==0) { 116 return zero; 117 } else { 118 dk=(double)k; 119 return dk*ln2_hi+dk*ln2_lo; 120 } 121 } 122 R = f*f*(0.5-0.33333333333333333*f); 123 if(k==0) return f-R; else {dk=(double)k; 124 return dk*ln2_hi-((R-dk*ln2_lo)-f);} 125 } 126 s = f/(2.0+f); 127 dk = (double)k; 114 128 z = s*s; 129 i = hx-0x6147a; 115 130 w = z*z; 116 t1 = w*(Lg2+w*(Lg4+w*Lg6)); 117 t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 118 R = t2 + t1; 119 dk = k; 120 121 return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi; 131 j = 0x6b851-hx; 132 t1= w*(Lg2+w*(Lg4+w*Lg6)); 133 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 134 i |= j; 135 R = t2+t1; 136 if(i>0) { 137 hfsq=0.5*f*f; 138 if(k==0) return f-(hfsq-s*(hfsq+R)); else 139 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); 140 } else { 141 if(k==0) return f-s*(f-R); else 142 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); 143 } 122 144 } 123 145
Note:
See TracChangeset
for help on using the changeset viewer.