Changeset 19520ba9 in mainline
 Timestamp:
 20210618T16:05:49Z (2 years ago)
 Children:
 f250c5a
 Parents:
 ec6081c
 gitauthor:
 Maurizio Lombardi <mlombard@…> (20181125 16:43:14)
 gitcommitter:
 Maurizio Lombardi <mlombard@…> (20210618 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.93147180369123816490e01, /* 3fe62e42 fee00000 */ 73 ln2_lo = 1.90821492927058770002e10, /* 3dea39ef 35793c76 */ 75 ln2_hi = 6.93147180369123816490e01, /* 3fe62e42 fee00000 */ 76 ln2_lo = 1.90821492927058770002e10, /* 3dea39ef 35793c76 */ 77 two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ 74 78 Lg1 = 6.666666666666735130e01, /* 3FE55555 55555593 */ 75 79 Lg2 = 3.999999999940941908e01, /* 3FD99999 9997FA04 */ … … 80 84 Lg7 = 1.479819860511658591e01; /* 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 (xx)/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 (xx)/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 = x1.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.50.33333333333333333*f); 123 if(k==0) return fR; else {dk=(double)k; 124 return dk*ln2_hi((Rdk*ln2_lo)f);} 125 } 126 s = f/(2.0+f); 127 dk = (double)k; 114 128 z = s*s; 129 i = hx0x6147a; 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 = 0x6b851hx; 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(hfsqs*(hfsq+R)); else 139 return dk*ln2_hi((hfsq(s*(hfsq+R)+dk*ln2_lo))f); 140 } else { 141 if(k==0) return fs*(fR); else 142 return dk*ln2_hi((s*(fR)dk*ln2_lo)f); 143 } 122 144 } 123 145
Note:
See TracChangeset
for help on using the changeset viewer.