1 | module math |
2 | |
3 | // log_n returns log base b of x |
4 | pub fn log_n(x f64, b f64) f64 { |
5 | y := log(x) |
6 | z := log(b) |
7 | return y / z |
8 | } |
9 | |
10 | // log10 returns the decimal logarithm of x. |
11 | // The special cases are the same as for log. |
12 | pub fn log10(x f64) f64 { |
13 | return log(x) * (1.0 / ln10) |
14 | } |
15 | |
16 | // log2 returns the binary logarithm of x. |
17 | // The special cases are the same as for log. |
18 | pub fn log2(x f64) f64 { |
19 | frac, exp := frexp(x) |
20 | // Make sure exact powers of two give an exact answer. |
21 | // Don't depend on log(0.5)*(1/ln2)+exp being exactly exp-1. |
22 | if frac == 0.5 { |
23 | return f64(exp - 1) |
24 | } |
25 | return log(frac) * (1.0 / ln2) + f64(exp) |
26 | } |
27 | |
28 | // log1p returns log(1+x) |
29 | pub fn log1p(x f64) f64 { |
30 | y := 1.0 + x |
31 | z := y - 1.0 |
32 | return log(y) - (z - x) / y // cancels errors with IEEE arithmetic |
33 | } |
34 | |
35 | // log_b returns the binary exponent of x. |
36 | // |
37 | // special cases are: |
38 | // log_b(±inf) = +inf |
39 | // log_b(0) = -inf |
40 | // log_b(nan) = nan |
41 | pub fn log_b(x f64) f64 { |
42 | if x == 0 { |
43 | return inf(-1) |
44 | } |
45 | if is_inf(x, 0) { |
46 | return inf(1) |
47 | } |
48 | if is_nan(x) { |
49 | return x |
50 | } |
51 | return f64(ilog_b_(x)) |
52 | } |
53 | |
54 | // ilog_b returns the binary exponent of x as an integer. |
55 | // |
56 | // special cases are: |
57 | // ilog_b(±inf) = max_i32 |
58 | // ilog_b(0) = min_i32 |
59 | // ilog_b(nan) = max_i32 |
60 | pub fn ilog_b(x f64) int { |
61 | if x == 0 { |
62 | return min_i32 |
63 | } |
64 | if is_nan(x) { |
65 | return max_i32 |
66 | } |
67 | if is_inf(x, 0) { |
68 | return max_i32 |
69 | } |
70 | return ilog_b_(x) |
71 | } |
72 | |
73 | // ilog_b returns the binary exponent of x. It assumes x is finite and |
74 | // non-zero. |
75 | fn ilog_b_(x_ f64) int { |
76 | x, exp := normalize(x_) |
77 | return int((f64_bits(x) >> shift) & mask) - bias + exp |
78 | } |
79 | |
80 | // log returns the natural logarithm of x |
81 | // |
82 | // Method : |
83 | // 1. Argument Reduction: find k and f such that |
84 | // x = 2^k * (1+f), |
85 | // where sqrt(2)/2 < 1+f < sqrt(2) . |
86 | // |
87 | // 2. Approximation of log(1+f). |
88 | // Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) |
89 | // = 2s + 2/3 s**3 + 2/5 s**5 + ....., |
90 | // = 2s + s*R |
91 | // We use a special Remez algorithm on [0,0.1716] to generate |
92 | // a polynomial of degree 14 to approximate R The maximum error |
93 | // of this polynomial approximation is bounded by 2**-58.45. In |
94 | // other words, |
95 | // 2 4 6 8 10 12 14 |
96 | // R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s |
97 | // (the values of Lg1 to Lg7 are listed in the program) |
98 | // and |
99 | // | 2 14 | -58.45 |
100 | // | Lg1*s +...+Lg7*s - R(z) | <= 2 |
101 | // | | |
102 | // Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. |
103 | // In order to guarantee error in log below 1ulp, we compute log |
104 | // by |
105 | // log(1+f) = f - s*(f - R) (if f is not too large) |
106 | // log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) |
107 | // |
108 | // 3. Finally, log(x) = k*ln2 + log(1+f). |
109 | // = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) |
110 | // Here ln2 is split into two floating point number: |
111 | // ln2_hi + ln2_lo, |
112 | // where n*ln2_hi is always exact for |n| < 2000. |
113 | // |
114 | // Special cases: |
115 | // log(x) is NaN with signal if x < 0 (including -inf) ; |
116 | // log(+inf) is +inf; log(0) is -inf with signal; |
117 | // log(NaN) is that NaN with no signal. |
118 | // |
119 | // Accuracy: |
120 | // according to an error analysis, the error is always less than |
121 | // 1 ulp (unit in the last place). |
122 | pub fn log(a f64) f64 { |
123 | ln2_hi := 6.93147180369123816490e-01 // 3fe62e42 fee00000 |
124 | ln2_lo := 1.90821492927058770002e-10 // 3dea39ef 35793c76 |
125 | l1 := 6.666666666666735130e-01 // 3FE55555 55555593 |
126 | l2 := 3.999999999940941908e-01 // 3FD99999 9997FA04 |
127 | l3 := 2.857142874366239149e-01 // 3FD24924 94229359 |
128 | l4 := 2.222219843214978396e-01 // 3FCC71C5 1D8E78AF |
129 | l5 := 1.818357216161805012e-01 // 3FC74664 96CB03DE |
130 | l6 := 1.531383769920937332e-01 // 3FC39A09 D078C69F |
131 | l7 := 1.479819860511658591e-01 // 3FC2F112 DF3E5244 |
132 | |
133 | x := a |
134 | if is_nan(x) || is_inf(x, 1) { |
135 | return x |
136 | } else if x < 0 { |
137 | return nan() |
138 | } else if x == 0 { |
139 | return inf(-1) |
140 | } |
141 | |
142 | mut f1, mut ki := frexp(x) |
143 | if f1 < sqrt2 / 2 { |
144 | f1 *= 2 |
145 | ki-- |
146 | } |
147 | |
148 | f := f1 - 1 |
149 | k := f64(ki) |
150 | |
151 | // compute |
152 | s := f / (2 + f) |
153 | s2 := s * s |
154 | s4 := s2 * s2 |
155 | t1 := s2 * (l1 + s4 * (l3 + s4 * (l5 + s4 * l7))) |
156 | t2 := s4 * (l2 + s4 * (l4 + s4 * l6)) |
157 | r := t1 + t2 |
158 | hfsq := 0.5 * f * f |
159 | return k * ln2_hi - ((hfsq - (s * (hfsq + r) + k * ln2_lo)) - f) |
160 | } |