/src/moddable/xs/tools/fdlibm/e_log10.c
Line | Count | Source (jump to first uncovered line) |
1 | | |
2 | | /* |
3 | | * ==================================================== |
4 | | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 | | * |
6 | | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 | | * Permission to use, copy, modify, and distribute this |
8 | | * software is freely granted, provided that this notice |
9 | | * is preserved. |
10 | | * ==================================================== |
11 | | */ |
12 | | |
13 | | /* |
14 | | * Return the base 10 logarithm of x. See e_log.c and k_log.h for most |
15 | | * comments. |
16 | | * |
17 | | * log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) |
18 | | * in not-quite-routine extra precision. |
19 | | */ |
20 | | |
21 | | #include "math_private.h" |
22 | | #include "k_log.h" |
23 | | |
24 | | static const double |
25 | | two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ |
26 | | ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ |
27 | | ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ |
28 | | log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ |
29 | | log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ |
30 | | |
31 | | static const double zero = 0.0; |
32 | | static volatile double vzero = 0.0; |
33 | | |
34 | | double |
35 | | __ieee754_log10(double x) |
36 | 97 | { |
37 | 97 | double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; |
38 | 97 | int32_t i,k,hx; |
39 | 97 | u_int32_t lx; |
40 | | |
41 | 97 | EXTRACT_WORDS(hx,lx,x); |
42 | | |
43 | 97 | k=0; |
44 | 97 | if (hx < 0x00100000) { /* x < 2**-1022 */ |
45 | 54 | if (((hx&0x7fffffff)|lx)==0) |
46 | 34 | return -two54/vzero; /* log(+-0)=-inf */ |
47 | 20 | if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
48 | 0 | k -= 54; x *= two54; /* subnormal number, scale up x */ |
49 | 0 | GET_HIGH_WORD(hx,x); |
50 | 0 | } |
51 | 43 | if (hx >= 0x7ff00000) return x+x; |
52 | 21 | if (hx == 0x3ff00000 && lx == 0) |
53 | 4 | return zero; /* log(1) = +0 */ |
54 | 17 | k += (hx>>20)-1023; |
55 | 17 | hx &= 0x000fffff; |
56 | 17 | i = (hx+0x95f64)&0x100000; |
57 | 17 | SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ |
58 | 17 | k += (i>>20); |
59 | 17 | y = (double)k; |
60 | 17 | f = x - 1.0; |
61 | 17 | hfsq = 0.5*f*f; |
62 | 17 | r = k_log1p(f); |
63 | | |
64 | | /* See e_log2.c for most details. */ |
65 | 17 | hi = f - hfsq; |
66 | 17 | SET_LOW_WORD(hi,0); |
67 | 17 | lo = (f - hi) - hfsq + r; |
68 | 17 | val_hi = hi*ivln10hi; |
69 | 17 | y2 = y*log10_2hi; |
70 | 17 | val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi; |
71 | | |
72 | | /* |
73 | | * Extra precision in for adding y*log10_2hi is not strictly needed |
74 | | * since there is no very large cancellation near x = sqrt(2) or |
75 | | * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs |
76 | | * with some parallelism and it reduces the error for many args. |
77 | | */ |
78 | 17 | w = y2 + val_hi; |
79 | 17 | val_lo += (y2 - w) + val_hi; |
80 | 17 | val_hi = w; |
81 | | |
82 | 17 | return val_lo + val_hi; |
83 | 21 | } |