1 module ctstdmath;
2 
3 import std.math : LN2, LOG2E, PI, PI_2, abs, fabs, sqrt, hypot, poly;
4 unittest
5 {
6     enum a = real(1).abs, b = real(1).fabs, c = real(1).sqrt, d = real(1).hypot(real(1)), e = real(1).poly([real(0)]);
7 }
8 unittest
9 {
10     enum a = real(1).cbrt, b = real(1).exp2, c = real(1).exp, d = real(1).lg, e = real(1).log, f = real(1).expm1;
11 }
12 
13 ///
14 real cbrt(in real a)
15 {
16     if (a < 0)
17         return -(-a).cbrt;
18     if (a.isNaN || a.isInfinity)
19         return a;
20     if (a == 0)
21         return 0;
22     // safe initial guess.
23     real x = a.sqrt, y = real.infinity;
24     while (x < y)
25     {
26         y = x;
27         x = (x * 2 + (a / (x * x))) / 3;
28     }
29     return x;
30 }
31 ///
32 unittest
33 {
34     static assert(8.cbrt == 2);
35     static assert(0.cbrt == 0);
36     static assert((-8).cbrt == -2);
37     static assert(real.nan.cbrt.isNaN);
38     static assert(real.infinity.cbrt.isInfinity);
39 }
40 
41 /// exponential function.
42 real exp(real x)
43 {
44     return (x * LOG2E).exp2;
45 }
46 
47 /// power of 2.
48 real exp2(real x)
49 {
50     if (x.isNaN) return x;
51     if (x < -(14).exp2small) return 0;
52     if (14.exp2small < x) return real.infinity;
53     real ret = 1;
54     if (0 < x)
55     {
56         foreach_reverse (byte i; 0..15)
57         {
58             if (x <= i.exp2small)
59                 continue;
60             x -= i.exp2small;
61             ret *= i.exp2exp2;
62         }
63     }
64     else
65     {
66         foreach_reverse (byte i; 0..15)
67         {
68             if (-(i.exp2small) <= x)
69                 continue;
70             x += i.exp2small;
71             ret /= i.exp2exp2;
72         }
73     }
74     if (ret == 0 || ret == real.infinity)
75         return ret;
76     if (+0.5 < x)
77     {
78         x -= 1;
79         ret *= 2;
80     }
81     if (x < -0.5)
82     {
83         x += 1;
84         ret /= 2;
85     }
86     return (x * LN2).expsmall * ret;
87 }
88 ///
89 unittest
90 {
91     static assert(real.nan.exp2.isNaN);
92     static assert(real.infinity.exp2.isInfinity);
93     static assert((-real.infinity).exp2 == 0);
94     static assert(0.exp2 == 1);
95     static assert(2.exp2 == 4);
96     static assert(8.exp2 == 256);
97     static assert((-1).exp2 == 0.5);
98 }
99 
100 private auto expsmall(in real x)
101 in
102 {
103     assert (-1 < x);
104     assert (x < +1);
105 }
106 body
107 {
108     real ret = 0, nextTerm = 1, old = real.nan;
109     size_t i;
110     while (ret != old)
111     {
112         old = ret;
113         ret += nextTerm;
114         i += 1;
115         nextTerm *= x / i;
116     }
117     return ret;
118 }
119 
120 private auto exp2small(byte n)
121 {
122     real ret = 1;
123     if (0 <= n)
124         foreach (i; 0..n)
125             ret *= 2;
126     else
127         foreach (i; 0..-n)
128             ret *= real(0.5);
129     return ret;
130 }
131 private auto doubling(real x, in ubyte n)
132 {
133     foreach (i; 0..n)
134         x *= x;
135     return x;
136 }
137 private auto exp2exp2(byte n)
138 {
139     return (0 <= n) ? 2.doubling(n) : 0.5.doubling(-n);
140 }
141 unittest
142 {
143     static assert(2.exp2small == 4);
144     static assert(2.exp2exp2 == 16);
145     static assert((-2).exp2small * 4 == 1);
146     static assert((-2).exp2exp2 * 16 == 1);
147 }
148 
149 /// equivalent to exp(x) - 1, but more accurate when |x| is small.
150 real expm1(real x)
151 {
152     if (!(-1 < x && x < 1))
153         return x.exp - 1;
154     real ret = 0, nextTerm = 1, old = real.nan;
155     size_t i;
156     while (ret != old)
157     {
158         i += 1;
159         nextTerm *= x / i;
160         old = ret;
161         ret += nextTerm;
162     }
163     return ret;
164 }
165 
166 /// Natural logarithm.
167 real log(real x)
168 {
169     return x.lg * LN2;
170 }
171 alias ln = log; /// ditto
172 
173 /// Base-2 logarithm.
174 real lg(real x)
175 {
176     real ret = 0;
177     if (x < 0) return real.nan;
178     if (x == 0) return -real.infinity;
179     if (x.isNaN) return x;
180     if (x.isInfinity) return real.infinity; // x = +inf
181     enum left = real(2) / 3, right = real(4) / 3;
182     if (right < x)
183         foreach_reverse (byte i; 0..15)
184         {
185             if (x / i.exp2exp2 < left)
186                 continue;
187             x /= i.exp2exp2;
188             ret += i.exp2small;
189         }
190     else if (x < left)
191         foreach_reverse (byte i; 0..15)
192         {
193             if (right < i.exp2exp2 * x)
194                 continue;
195             x *= i.exp2exp2;
196             ret -= i.exp2small;
197         }
198     return x.logsmall / LN2 + ret;
199 }
200 /// ditto
201 alias log2 = lg;
202 unittest
203 {
204     static assert((-1).lg.isNaN);
205     static assert(real.infinity.lg.isInfinity);
206     static assert(real.nan.lg.isNaN);
207     static assert(0.lg < 0 && 0.lg.isInfinity);
208     static assert(0.5.lg == -1);
209     static assert(1.lg == 0);
210     static assert(2.lg == 1);
211 }
212 
213 private real logsmall(real x)
214 in
215 {
216     assert (0.5 < x);
217     assert (x < 1.5);
218 }
219 body
220 {
221     return log1p(x - 1);
222 }
223 
224 /// equivalent to log(x + 1), but more accurate when |x| is small.
225 real log1p(real x)
226 {
227     real ret = 0, power = -1, old = real.nan;
228     size_t i;
229     while (ret != old)
230     {
231         i += 1;
232         power *= -x;
233         old = ret;
234         ret += power / i;
235     }
236     return ret;
237 }
238 
239 /// floor function.
240 real floor(real x)
241 {
242     if (x.isNaN)
243         return x;
244     if (x <= -0x1p+63 || +0x1p+63 <= x)
245         return x;
246     if (x < 0)
247         return -(-x).ceil;
248     ulong left = 0, right = 1UL << 63;
249     // binary search [)
250     foreach (i; 0..63)
251     {
252         auto mid = (left & right) + ((left ^ right) >> 1);
253         if (x < mid)
254             right = mid;
255         else
256             left = mid;
257     }
258     return left;
259 }
260 ///
261 unittest
262 {
263     static assert((-0.01).floor == -1);
264     static assert(real(0).floor == 0);
265     static assert(0.01.floor == 0);
266 }
267 
268 /// ceiling function.
269 real ceil(real x)
270 {
271     if (x.isNaN)
272         return x;
273     if (x <= -0x1p+63 || +0x1p+63 <= x)
274         return x;
275     if (x <= 0)
276         return -(-x).floor;
277     ulong left = 0, right = 1UL << 63;
278     // binary search (]
279     foreach (i; 0..63)
280     {
281         auto mid = (left & right) + ((left ^ right) >> 1);
282         if (x <= mid)
283             right = mid;
284         else
285             left = mid;
286     }
287     return right;
288 }
289 /// ditto
290 alias ceiling = ceil;
291 ///
292 unittest
293 {
294     static assert((-0.01).ceil == 0);
295     static assert(real(0).ceil == 0);
296     static assert(0.01.ceil == 1);
297 }
298 
299 /// reduce x into [-y .. y] modulo 2y.
300 real modHalf(real x, real y)
301 {
302     if (x.isNaN || y.isNaN || y == 0)
303         return real.nan;
304     if (y.isInfinity)
305         return x;
306     if (y < 0)
307         return x.modHalf(-y);
308     if (x < 0)
309         return -(-x).modHalf(y);
310     void dfs(real mod)
311     {
312         if (mod * 2 <= x)
313             return dfs(mod * 2);
314         if (mod <= x)
315             x -= mod;
316     }
317     dfs(y * 2);
318     if (y < x)
319         return x - 2 * y;
320     return x;
321 }
322 ///
323 unittest
324 {
325     static assert(0.modHalf(1) == 0);
326     static assert(0.5.modHalf(1) == 0.5);
327     static assert(1.modHalf(1) == 1);
328     static assert(1.5.modHalf(1) == -0.5);
329     static assert(2.modHalf(1) == 0);
330     static assert(2.5.modHalf(1) == 0.5);
331 }
332 
333 /// cosine.
334 real cos(real x)
335 {
336     return (x.modHalf(PI) + PI_2).sin;
337 }
338 ///
339 unittest
340 {
341     static assert(PI_2.cos == 0);
342 }
343 
344 /// sine.
345 real sin(real x)
346 {
347     x = x.modHalf(PI);
348     if (x < 0)
349         return -(-x).sin;
350     if (PI < x * 2)
351         return (PI - x).sin;
352     real ret = 0, nextTerm = x, old = real.nan;
353     size_t i;
354     x = -(x * x);
355     while (ret != old)
356     {
357         old = ret;
358         ret += nextTerm;
359         i += 2;
360         nextTerm *= x / (i * (i + 1));
361     }
362     return ret;
363 }
364 ///
365 unittest
366 {
367     static assert(0.sin == 0);
368 }
369 
370 /// tangent.
371 real tan(real x)
372 {
373     return x.sin / x.cos;
374 }
375 /// cotangent.
376 real cot(real x)
377 {
378     return x.cos / x.sin;
379 }
380 /// secant.
381 real sec(real x)
382 {
383     return 1 / x.cos;
384 }
385 /// cosecant.
386 real csc(real x)
387 {
388     return 1 / x.sin;
389 }
390 /// ditto
391 alias cosec = csc;
392 
393 /// arctangent.
394 real atan(real x)
395 {
396     if (x < 0)
397         return -(-x).atan;
398     if (x.isInfinity)
399         return PI_2;
400     if (x == 1)
401         return PI_2 / 2;
402     if (1 < x)
403         return PI_2 - (1 / x).atan;
404     if (1 < 3 * x * x) // without this, very slow when |x| ~ 1.
405         return PI_2 / 3 + ((x - 1 / real(3).sqrt) / (1 + x / real(3).sqrt)).atan;
406     real ret = 0, numeratorX = x, old = real.nan;
407     x *= x;
408     real i = 2;
409     while (ret != old)
410     {
411         old = ret;
412         ret += numeratorX * (i + 1 - (i - 1) * x) / (i * i - 1);
413         i += 4;
414         numeratorX *= x * x;
415     }
416     import std.experimental.logger;
417     i.trace;
418     return ret;
419 }
420 
421 ///
422 bool isInfinity(real x)
423 {
424     return x < -real.max || real.max < x;
425 }
426 ///
427 unittest
428 {
429     static assert(!isInfinity(real.init));
430     static assert(!isInfinity(-real.init));
431     static assert(!isInfinity(real.nan));
432     static assert(!isInfinity(-real.nan));
433     static assert(isInfinity(real.infinity));
434     static assert(isInfinity(-real.infinity));
435     static assert(isInfinity(-1.0L / 0.0L));
436     static assert(!isInfinity(real(0)));
437 }
438 
439 ///
440 bool isNaN(real x)
441 {
442     return !(0 <= x) && !(x <= 0);
443 }
444 ///
445 unittest
446 {
447     static assert(isNaN(real.nan));
448     static assert(isNaN(-real.nan));
449     static assert(!isNaN(real.infinity));
450     static assert(!isNaN(-real.infinity));
451     static assert(!isNaN(-1.0L / 0.0L));
452     static assert(!isNaN(real(0)));
453 }