1 | module math |
2 | |
3 | const ( |
4 | pow10tab = [f64(1e+00), 1e+01, 1e+02, 1e+03, 1e+04, 1e+05, 1e+06, 1e+07, 1e+08, 1e+09, |
5 | 1e+10, 1e+11, 1e+12, 1e+13, 1e+14, 1e+15, 1e+16, 1e+17, 1e+18, 1e+19, 1e+20, 1e+21, 1e+22, |
6 | 1e+23, 1e+24, 1e+25, 1e+26, 1e+27, 1e+28, 1e+29, 1e+30, 1e+31] |
7 | pow10postab32 = [f64(1e+00), 1e+32, 1e+64, 1e+96, 1e+128, 1e+160, 1e+192, 1e+224, 1e+256, 1e+288] |
8 | pow10negtab32 = [f64(1e-00), 1e-32, 1e-64, 1e-96, 1e-128, 1e-160, 1e-192, 1e-224, 1e-256, 1e-288, |
9 | 1e-320] |
10 | ) |
11 | |
12 | // powf returns base raised to the provided power. (float32) |
13 | [inline] |
14 | pub fn powf(a f32, b f32) f32 { |
15 | return f32(pow(a, b)) |
16 | } |
17 | |
18 | // pow10 returns 10**n, the base-10 exponential of n. |
19 | // |
20 | // special cases are: |
21 | // pow10(n) = 0 for n < -323 |
22 | // pow10(n) = +inf for n > 308 |
23 | pub fn pow10(n int) f64 { |
24 | if 0 <= n && n <= 308 { |
25 | return math.pow10postab32[u32(n) / 32] * math.pow10tab[u32(n) % 32] |
26 | } |
27 | if -323 <= n && n <= 0 { |
28 | return math.pow10negtab32[u32(-n) / 32] / math.pow10tab[u32(-n) % 32] |
29 | } |
30 | // n < -323 || 308 < n |
31 | if n > 0 { |
32 | return inf(1) |
33 | } |
34 | // n < -323 |
35 | return 0.0 |
36 | } |
37 | |
38 | // powi returns base raised to power (a**b) as an integer (i64) |
39 | // |
40 | // special case: |
41 | // powi(a, b) = -1 for a = 0 and b < 0 |
42 | pub fn powi(a i64, b i64) i64 { |
43 | mut b_ := b |
44 | mut p := a |
45 | mut v := i64(1) |
46 | |
47 | if b_ < 0 { // exponent < 0 |
48 | if a == 0 { |
49 | return -1 // division by 0 |
50 | } |
51 | return if a * a != 1 { |
52 | 0 |
53 | } else { |
54 | if (b_ & 1) > 0 { |
55 | a |
56 | } else { |
57 | 1 |
58 | } |
59 | } |
60 | } |
61 | |
62 | for ; b_ > 0; { |
63 | if b_ & 1 > 0 { |
64 | v *= p |
65 | } |
66 | p *= p |
67 | b_ >>= 1 |
68 | } |
69 | |
70 | return v |
71 | } |
72 | |
73 | // pow returns the base x, raised to the provided power y. (float64) |
74 | // |
75 | // todo(playXE): make this function work on JS backend, probably problem of JS codegen that it does not work. |
76 | pub fn pow(x f64, y f64) f64 { |
77 | if y == 0 || x == 1 { |
78 | return 1 |
79 | } else if y == 1 { |
80 | return x |
81 | } else if is_nan(x) || is_nan(y) { |
82 | return nan() |
83 | } else if y == 2 { |
84 | return x * x |
85 | } else if y == 3 { |
86 | return x * x * x |
87 | } else if x == 0 { |
88 | if y < 0 { |
89 | if is_odd_int(y) { |
90 | return copysign(inf(1), x) |
91 | } |
92 | return inf(1) |
93 | } else if y > 0 { |
94 | if is_odd_int(y) { |
95 | return x |
96 | } |
97 | return 0 |
98 | } |
99 | } else if is_inf(y, 0) { |
100 | if x == -1 { |
101 | return 1 |
102 | } else if (abs(x) < 1) == is_inf(y, 1) { |
103 | return 0 |
104 | } else { |
105 | return inf(1) |
106 | } |
107 | } else if is_inf(x, 0) { |
108 | if is_inf(x, -1) { |
109 | return pow(1 / x, -y) |
110 | } |
111 | |
112 | if y < 0 { |
113 | return 0 |
114 | } else if y > 0 { |
115 | return inf(1) |
116 | } |
117 | } else if y == 0.5 { |
118 | return sqrt(x) |
119 | } else if y == -0.5 { |
120 | return 1 / sqrt(x) |
121 | } |
122 | mut yi, mut yf := modf(abs(y)) |
123 | |
124 | if yf != 0 && x < 0 { |
125 | return nan() |
126 | } |
127 | if yi >= (u64(1) << 63) { |
128 | // yi is a large even int that will lead to overflow (or underflow to 0) |
129 | // for all x except -1 (x == 1 was handled earlier) |
130 | |
131 | if x == -1 { |
132 | return 1 |
133 | } else if (abs(x) < 1) == (y > 0) { |
134 | return 0 |
135 | } else { |
136 | return inf(1) |
137 | } |
138 | } |
139 | |
140 | if yf == 0.0 { |
141 | mut result := x |
142 | for _ in 1 .. i64(yi) { |
143 | result *= x |
144 | } |
145 | if y > 0 { |
146 | return result |
147 | } |
148 | return 1 / result |
149 | } |
150 | |
151 | // ans = a1 * 2**ae (= 1 for now). |
152 | mut a1 := 1.0 |
153 | mut ae := 0 |
154 | |
155 | // ans *= x**yf |
156 | if yf != 0 { |
157 | if yf > 0.5 { |
158 | yf-- |
159 | yi++ |
160 | } |
161 | |
162 | a1 = exp(yf * log(x)) |
163 | } |
164 | |
165 | // ans *= x**yi |
166 | // by multiplying in successive squarings |
167 | // of x according to bits of yi. |
168 | // accumulate powers of two into exp. |
169 | mut x1, mut xe := frexp(x) |
170 | |
171 | for i := i64(yi); i != 0; i >>= 1 { |
172 | // these series of casts is a little weird but we have to do them to prevent left shift of negative error |
173 | if xe < int(u32(u32(-1) << 12)) || 1 << 12 < xe { |
174 | // catch xe before it overflows the left shift below |
175 | // Since i !=0 it has at least one bit still set, so ae will accumulate xe |
176 | // on at least one more iteration, ae += xe is a lower bound on ae |
177 | // the lower bound on ae exceeds the size of a float64 exp |
178 | // so the final call to Ldexp will produce under/overflow (0/Inf) |
179 | ae += xe |
180 | break |
181 | } |
182 | if i & 1 == 1 { |
183 | a1 *= x1 |
184 | ae += xe |
185 | } |
186 | x1 *= x1 |
187 | xe <<= 1 |
188 | if x1 < .5 { |
189 | x1 += x1 |
190 | xe-- |
191 | } |
192 | } |
193 | |
194 | // ans = a1*2**ae |
195 | // if y < 0 { ans = 1 / ans } |
196 | // but in the opposite order |
197 | if y < 0 { |
198 | a1 = 1 / a1 |
199 | ae = -ae |
200 | } |
201 | return ldexp(a1, ae) |
202 | } |