Go to most recent revision | Details | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
4349 | Serge | 1 | |
2 | /* |
||
3 | * ==================================================== |
||
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
||
5 | * |
||
6 | * Developed at SunPro, 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 |
||
15 | * |
||
16 | * Method : |
||
17 | * Let log10_2hi = leading 40 bits of log10(2) and |
||
18 | * log10_2lo = log10(2) - log10_2hi, |
||
19 | * ivln10 = 1/log(10) rounded. |
||
20 | * Then |
||
21 | * n = ilogb(x), |
||
22 | * if(n<0) n = n+1; |
||
23 | * x = scalbn(x,-n); |
||
24 | * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) |
||
25 | * |
||
26 | * Note 1: |
||
27 | * To guarantee log10(10**n)=n, where 10**n is normal, the rounding |
||
28 | * mode must set to Round-to-Nearest. |
||
29 | * Note 2: |
||
30 | * [1/log(10)] rounded to 53 bits has error .198 ulps; |
||
31 | * log10 is monotonic at all binary break points. |
||
32 | * |
||
33 | * Special cases: |
||
34 | * log10(x) is NaN with signal if x < 0; |
||
35 | * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; |
||
36 | * log10(NaN) is that NaN with no signal; |
||
37 | * log10(10**N) = N for N=0,1,...,22. |
||
38 | * |
||
39 | * Constants: |
||
40 | * The hexadecimal values are the intended ones for the following constants. |
||
41 | * The decimal values may be used, provided that the compiler will convert |
||
42 | * from decimal to binary accurately enough to produce the hexadecimal values |
||
43 | * shown. |
||
44 | */ |
||
45 | |||
46 | |||
47 | |||
48 | |||
49 | |||
50 | |||
51 | static const double |
||
52 | #else |
||
53 | static double |
||
54 | #endif |
||
55 | two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ |
||
56 | ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */ |
||
57 | log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ |
||
58 | log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ |
||
59 | |||
60 | |||
61 | static const double zero = 0.0; |
||
62 | #else |
||
63 | static double zero = 0.0; |
||
64 | #endif |
||
65 | |||
66 | |||
67 | double __ieee754_log10(double x) |
||
68 | #else |
||
69 | double __ieee754_log10(x) |
||
70 | double x; |
||
71 | #endif |
||
72 | { |
||
73 | double y,z; |
||
74 | __int32_t i,k,hx; |
||
75 | __uint32_t lx; |
||
76 | |||
77 | |||
78 | |||
79 | |||
80 | if (hx < 0x00100000) { /* x < 2**-1022 */ |
||
81 | if (((hx&0x7fffffff)|lx)==0) |
||
82 | return -two54/zero; /* log(+-0)=-inf */ |
||
83 | if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
||
84 | k -= 54; x *= two54; /* subnormal number, scale up x */ |
||
85 | GET_HIGH_WORD(hx,x); |
||
86 | } |
||
87 | if (hx >= 0x7ff00000) return x+x; |
||
88 | k += (hx>>20)-1023; |
||
89 | i = ((__uint32_t)k&0x80000000)>>31; |
||
90 | hx = (hx&0x000fffff)|((0x3ff-i)<<20); |
||
91 | y = (double)(k+i); |
||
92 | SET_HIGH_WORD(x,hx); |
||
93 | z = y*log10_2lo + ivln10*__ieee754_log(x); |
||
94 | return z+y*log10_2hi; |
||
95 | } |
||
96 | |||
97 | |||
98 | ><20); |