1 | module math |
2 | |
3 | import math.internal |
4 | |
5 | const ( |
6 | f64_max_exp = f64(1024) |
7 | f64_min_exp = f64(-1021) |
8 | threshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF |
9 | ln2_x56 = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1 |
10 | ln2_halfx3 = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73 |
11 | ln2_half = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef |
12 | ln2hi = 6.93147180369123816490e-01 // 0x3fe62e42fee00000 |
13 | ln2lo = 1.90821492927058770002e-10 // 0x3dea39ef35793c76 |
14 | inv_ln2 = 1.44269504088896338700e+00 // 0x3ff71547652b82fe |
15 | // scaled coefficients related to expm1 |
16 | expm1_q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4 |
17 | expm1_q2 = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585 |
18 | expm1_q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7 |
19 | expm1_q4 = 4.00821782732936239552e-06 // 0x3ED0CFCA86E65239 |
20 | expm1_q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D |
21 | ) |
22 | |
23 | // exp returns e**x, the base-e exponential of x. |
24 | // |
25 | // special cases are: |
26 | // exp(+inf) = +inf |
27 | // exp(nan) = nan |
28 | // Very large values overflow to 0 or +inf. |
29 | // Very small values underflow to 1. |
30 | pub fn exp(x f64) f64 { |
31 | log2e := 1.44269504088896338700e+00 |
32 | overflow := 7.09782712893383973096e+02 |
33 | underflow := -7.45133219101941108420e+02 |
34 | near_zero := 1.0 / (1 << 28) // 2**-28 |
35 | // special cases |
36 | if is_nan(x) || is_inf(x, 1) { |
37 | return x |
38 | } |
39 | if is_inf(x, -1) { |
40 | return 0.0 |
41 | } |
42 | if x > overflow { |
43 | return inf(1) |
44 | } |
45 | if x < underflow { |
46 | return 0.0 |
47 | } |
48 | if -near_zero < x && x < near_zero { |
49 | return 1.0 + x |
50 | } |
51 | // reduce; computed as r = hi - lo for extra precision. |
52 | mut k := 0 |
53 | if x < 0 { |
54 | k = int(log2e * x - 0.5) |
55 | } |
56 | if x > 0 { |
57 | k = int(log2e * x + 0.5) |
58 | } |
59 | hi := x - f64(k) * math.ln2hi |
60 | lo := f64(k) * math.ln2lo |
61 | // compute |
62 | return expmulti(hi, lo, k) |
63 | } |
64 | |
65 | // exp2 returns 2**x, the base-2 exponential of x. |
66 | // |
67 | // special cases are the same as exp. |
68 | pub fn exp2(x f64) f64 { |
69 | overflow := 1.0239999999999999e+03 |
70 | underflow := -1.0740e+03 |
71 | if is_nan(x) || is_inf(x, 1) { |
72 | return x |
73 | } |
74 | if is_inf(x, -1) { |
75 | return 0 |
76 | } |
77 | if x > overflow { |
78 | return inf(1) |
79 | } |
80 | if x < underflow { |
81 | return 0 |
82 | } |
83 | // argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2. |
84 | // computed as r = hi - lo for extra precision. |
85 | mut k := 0 |
86 | if x > 0 { |
87 | k = int(x + 0.5) |
88 | } |
89 | if x < 0 { |
90 | k = int(x - 0.5) |
91 | } |
92 | mut t := x - f64(k) |
93 | hi := t * math.ln2hi |
94 | lo := -t * math.ln2lo |
95 | // compute |
96 | return expmulti(hi, lo, k) |
97 | } |
98 | |
99 | // ldexp calculates frac*(2**exp) |
100 | pub fn ldexp(frac f64, exp int) f64 { |
101 | return scalbn(frac, exp) |
102 | } |
103 | |
104 | // frexp breaks f into a normalized fraction |
105 | // and an integral power of two. |
106 | // It returns frac and exp satisfying f == frac × 2**exp, |
107 | // with the absolute value of frac in the interval [½, 1). |
108 | // |
109 | // special cases are: |
110 | // frexp(±0) = ±0, 0 |
111 | // frexp(±inf) = ±inf, 0 |
112 | // frexp(nan) = nan, 0 |
113 | // pub fn frexp(f f64) (f64, int) { |
114 | // mut y := f64_bits(x) |
115 | // ee := int((y >> 52) & 0x7ff) |
116 | // // special cases |
117 | // if ee == 0 { |
118 | // if x != 0.0 { |
119 | // x1p64 := f64_from_bits(0x43f0000000000000) |
120 | // z,e_ := frexp(x * x1p64) |
121 | // return z,e_ - 64 |
122 | // } |
123 | // return x,0 |
124 | // } else if ee == 0x7ff { |
125 | // return x,0 |
126 | // } |
127 | // e_ := ee - 0x3fe |
128 | // y &= 0x800fffffffffffff |
129 | // y |= 0x3fe0000000000000 |
130 | // return f64_from_bits(y),e_ |
131 | pub fn frexp(x f64) (f64, int) { |
132 | mut y := f64_bits(x) |
133 | ee := int((y >> 52) & 0x7ff) |
134 | if ee == 0 { |
135 | if x != 0.0 { |
136 | x1p64 := f64_from_bits(u64(0x43f0000000000000)) |
137 | z, e_ := frexp(x * x1p64) |
138 | return z, e_ - 64 |
139 | } |
140 | return x, 0 |
141 | } else if ee == 0x7ff { |
142 | return x, 0 |
143 | } |
144 | e_ := ee - 0x3fe |
145 | y &= u64(0x800fffffffffffff) |
146 | y |= u64(0x3fe0000000000000) |
147 | return f64_from_bits(y), e_ |
148 | } |
149 | |
150 | // expm1 calculates e**x - 1 |
151 | // special cases are: |
152 | // expm1(+inf) = +inf |
153 | // expm1(-inf) = -1 |
154 | // expm1(nan) = nan |
155 | pub fn expm1(x f64) f64 { |
156 | if is_inf(x, 1) || is_nan(x) { |
157 | return x |
158 | } |
159 | if is_inf(x, -1) { |
160 | return f64(-1) |
161 | } |
162 | // FIXME: this should be improved |
163 | if abs(x) < ln2 { // Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ... |
164 | mut i := 1.0 |
165 | mut sum := x |
166 | mut term := x / 1.0 |
167 | i++ |
168 | term *= x / f64(i) |
169 | sum += term |
170 | for abs(term) > abs(sum) * internal.f64_epsilon { |
171 | i++ |
172 | term *= x / f64(i) |
173 | sum += term |
174 | } |
175 | return sum |
176 | } else { |
177 | return exp(x) - 1 |
178 | } |
179 | } |
180 | |
181 | fn expmulti(hi f64, lo f64, k int) f64 { |
182 | exp_p1 := 1.66666666666666657415e-01 // 0x3FC55555; 0x55555555 |
183 | exp_p2 := -2.77777777770155933842e-03 // 0xBF66C16C; 0x16BEBD93 |
184 | exp_p3 := 6.61375632143793436117e-05 // 0x3F11566A; 0xAF25DE2C |
185 | exp_p4 := -1.65339022054652515390e-06 // 0xBEBBBD41; 0xC5D26BF1 |
186 | exp_p5 := 4.13813679705723846039e-08 // 0x3E663769; 0x72BEA4D0 |
187 | r := hi - lo |
188 | t := r * r |
189 | c := r - t * (exp_p1 + t * (exp_p2 + t * (exp_p3 + t * (exp_p4 + t * exp_p5)))) |
190 | y := 1 - ((lo - (r * c) / (2 - c)) - hi) |
191 | // TODO(rsc): make sure ldexp can handle boundary k |
192 | return ldexp(y, k) |
193 | } |