1 /* This Source Code Form is subject to the terms of the Mozilla Public
2  * License, v. 2.0. If a copy of the MPL was not distributed with this
3  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
4 
5 /**
6  * This package provides mathematical functions.
7  *
8  * The $(D_PSYMBOL tanya.math) package itself provides only representation
9  * functions for built-in types, such as functions that provide information
10  * about internal representation of floating-point numbers and low-level
11  * operatons on these. Actual mathematical functions and additional types can
12  * be found in its submodules. $(D_PSYMBOL tanya.math) doesn't import any
13  * submodules publically, they should be imported explicitly.
14  *
15  * Copyright: Eugene Wissner 2016-2020.
16  * License: $(LINK2 https://www.mozilla.org/en-US/MPL/2.0/,
17  *                  Mozilla Public License, v. 2.0).
18  * Authors: $(LINK2 mailto:info@caraus.de, Eugene Wissner)
19  * Source: $(LINK2 https://github.com/caraus-ecms/tanya/blob/master/source/tanya/math/package.d,
20  *                 tanya/math/package.d)
21  */
22 module tanya.math;
23 
24 import tanya.math.nbtheory;
25 import tanya.meta.trait;
26 import tanya.meta.transform;
27 
28 /// Floating-point number precisions according to IEEE-754.
29 enum IEEEPrecision : ubyte
30 {
31     single = 4, /// Single precision: 64-bit.
32     double_ = 8, /// Single precision: 64-bit.
33     doubleExtended = 10, /// Double extended precision: 80-bit.
34 }
35 
36 /**
37  * Tests the precision of floating-point type $(D_PARAM F).
38  *
39  * For $(D_KEYWORD float) $(D_PSYMBOL ieeePrecision) always evaluates to
40  * $(D_INLINECODE IEEEPrecision.single); for $(D_KEYWORD double) - to
41  * $(D_INLINECODE IEEEPrecision.double). It returns different values only
42  * for $(D_KEYWORD real), since $(D_KEYWORD real) is a platform-dependent type.
43  *
44  * If $(D_PARAM F) is a $(D_KEYWORD real) and the target platform isn't
45  * currently supported, static assertion error will be raised (you can use
46  * $(D_INLINECODE is(typeof(ieeePrecision!F))) for testing the platform support
47  * without a compilation error).
48  *
49  * Params:
50  *  F = Type to be tested.
51  *
52  * Returns: Precision according to IEEE-754.
53  *
54  * See_Also: $(D_PSYMBOL IEEEPrecision).
55  */
56 template ieeePrecision(F)
57 if (isFloatingPoint!F)
58 {
59     static if (F.sizeof == float.sizeof)
60     {
61         enum IEEEPrecision ieeePrecision = IEEEPrecision.single;
62     }
63     else static if (F.sizeof == double.sizeof)
64     {
65         enum IEEEPrecision ieeePrecision = IEEEPrecision.double_;
66     }
67     else version (X86)
68     {
69         enum IEEEPrecision ieeePrecision = IEEEPrecision.doubleExtended;
70     }
71     else version (X86_64)
72     {
73         enum IEEEPrecision ieeePrecision = IEEEPrecision.doubleExtended;
74     }
75     else
76     {
77         static assert(false, "Unsupported IEEE 754 floating point precision");
78     }
79 }
80 
81 ///
82 @nogc nothrow pure @safe unittest
83 {
84     static assert(ieeePrecision!float == IEEEPrecision.single);
85     static assert(ieeePrecision!double == IEEEPrecision.double_);
86 }
87 
88 package(tanya) union FloatBits(F)
89 {
90     Unqual!F floating;
91     static if (ieeePrecision!F == IEEEPrecision.single)
92     {
93         uint integral;
94         enum uint expMask = 0x7f800000;
95     }
96     else static if (ieeePrecision!F == IEEEPrecision.double_)
97     {
98         ulong integral;
99         enum ulong expMask = 0x7ff0000000000000;
100     }
101     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
102     {
103         struct // Little-endian.
104         {
105             ulong mantissa;
106             ushort exp;
107         }
108         enum ulong mantissaMask = 0x7fffffffffffffff;
109         enum uint expMask = 0x7fff;
110     }
111     else
112     {
113         static assert(false, "Unsupported IEEE 754 floating point precision");
114     }
115 }
116 
117 /**
118  * Floating-point number classifications.
119  */
120 enum FloatingPointClass : ubyte
121 {
122     /**
123      * Not a Number.
124      *
125      * See_Also: $(D_PSYMBOL isNaN).
126      */
127     nan,
128 
129     /// Zero.
130     zero,
131 
132     /**
133      * Infinity.
134      *
135      * See_Also: $(D_PSYMBOL isInfinity).
136      */
137     infinite,
138 
139     /**
140      * Denormalized number.
141      *
142      * See_Also: $(D_PSYMBOL isSubnormal).
143      */
144     subnormal,
145 
146     /**
147      * Normalized number.
148      *
149      * See_Also: $(D_PSYMBOL isNormal).
150      */
151     normal,
152 }
153 
154 /**
155  * Returns whether $(D_PARAM x) is a NaN, zero, infinity, subnormal or
156  * normalized number.
157  *
158  * This function doesn't distinguish between negative and positive infinity,
159  * negative and positive NaN or negative and positive zero.
160  *
161  * Params:
162  *  F = Type of the floating point number.
163  *  x = Floating point number.
164  *
165  * Returns: Classification of $(D_PARAM x).
166  */
167 FloatingPointClass classify(F)(F x)
168 if (isFloatingPoint!F)
169 {
170     if (x == 0)
171     {
172         return FloatingPointClass.zero;
173     }
174     FloatBits!F bits;
175     bits.floating = abs(x);
176 
177     static if (ieeePrecision!F == IEEEPrecision.single)
178     {
179         if (bits.integral > bits.expMask)
180         {
181             return FloatingPointClass.nan;
182         }
183         else if (bits.integral == bits.expMask)
184         {
185             return FloatingPointClass.infinite;
186         }
187         else if (bits.integral < (1 << 23))
188         {
189             return FloatingPointClass.subnormal;
190         }
191     }
192     else static if (ieeePrecision!F == IEEEPrecision.double_)
193     {
194         if (bits.integral > bits.expMask)
195         {
196             return FloatingPointClass.nan;
197         }
198         else if (bits.integral == bits.expMask)
199         {
200             return FloatingPointClass.infinite;
201         }
202         else if (bits.integral < (1L << 52))
203         {
204             return FloatingPointClass.subnormal;
205         }
206     }
207     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
208     {
209         if (bits.exp == bits.expMask)
210         {
211             if ((bits.mantissa & bits.mantissaMask) == 0)
212             {
213                 return FloatingPointClass.infinite;
214             }
215             else
216             {
217                 return FloatingPointClass.nan;
218             }
219         }
220         else if (bits.exp == 0)
221         {
222             return FloatingPointClass.subnormal;
223         }
224         else if (bits.mantissa < (1L << 63)) // "Unnormal".
225         {
226             return FloatingPointClass.nan;
227         }
228     }
229 
230     return FloatingPointClass.normal;
231 }
232 
233 ///
234 @nogc nothrow pure @safe unittest
235 {
236     assert(classify(0.0) == FloatingPointClass.zero);
237     assert(classify(double.nan) == FloatingPointClass.nan);
238     assert(classify(double.infinity) == FloatingPointClass.infinite);
239     assert(classify(-double.infinity) == FloatingPointClass.infinite);
240     assert(classify(1.4) == FloatingPointClass.normal);
241     assert(classify(1.11254e-307 / 10) == FloatingPointClass.subnormal);
242 
243     assert(classify(0.0f) == FloatingPointClass.zero);
244     assert(classify(float.nan) == FloatingPointClass.nan);
245     assert(classify(float.infinity) == FloatingPointClass.infinite);
246     assert(classify(-float.infinity) == FloatingPointClass.infinite);
247     assert(classify(0.3) == FloatingPointClass.normal);
248     assert(classify(5.87747e-38f / 10) == FloatingPointClass.subnormal);
249 
250     assert(classify(0.0L) == FloatingPointClass.zero);
251     assert(classify(real.nan) == FloatingPointClass.nan);
252     assert(classify(real.infinity) == FloatingPointClass.infinite);
253     assert(classify(-real.infinity) == FloatingPointClass.infinite);
254 }
255 
256 /**
257  * Determines whether $(D_PARAM x) is a finite number.
258  *
259  * Params:
260  *  F = Type of the floating point number.
261  *  x = Floating point number.
262  *
263  * Returns: $(D_KEYWORD true) if $(D_PARAM x) is a finite number,
264  *          $(D_KEYWORD false) otherwise.
265  *
266  * See_Also: $(D_PSYMBOL isInfinity).
267  */
268 bool isFinite(F)(F x)
269 if (isFloatingPoint!F)
270 {
271     FloatBits!F bits;
272     static if (ieeePrecision!F == IEEEPrecision.single
273             || ieeePrecision!F == IEEEPrecision.double_)
274     {
275         bits.floating = x;
276         bits.integral &= bits.expMask;
277         return bits.integral != bits.expMask;
278     }
279     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
280     {
281         bits.floating = abs(x);
282         return (bits.exp != bits.expMask)
283             && (bits.exp == 0 || bits.mantissa >= (1L << 63));
284     }
285 }
286 
287 ///
288 @nogc nothrow pure @safe unittest
289 {
290     assert(!isFinite(float.infinity));
291     assert(!isFinite(-double.infinity));
292     assert(isFinite(0.0));
293     assert(!isFinite(float.nan));
294     assert(isFinite(5.87747e-38f / 10));
295     assert(isFinite(1.11254e-307 / 10));
296     assert(isFinite(0.5));
297 }
298 
299 /**
300  * Determines whether $(D_PARAM x) is $(B n)ot $(B a) $(B n)umber (NaN).
301  *
302  * Params:
303  *  F = Type of the floating point number.
304  *  x = Floating point number.
305  *
306  * Returns: $(D_KEYWORD true) if $(D_PARAM x) is not a number,
307  *          $(D_KEYWORD false) otherwise.
308  */
309 bool isNaN(F)(F x)
310 if (isFloatingPoint!F)
311 {
312     FloatBits!F bits;
313     bits.floating = abs(x);
314 
315     static if (ieeePrecision!F == IEEEPrecision.single
316             || ieeePrecision!F == IEEEPrecision.double_)
317     {
318         return bits.integral > bits.expMask;
319     }
320     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
321     {
322         const maskedMantissa = bits.mantissa & bits.mantissaMask;
323         if ((bits.exp == bits.expMask && maskedMantissa != 0)
324          || ((bits.exp != 0) && (bits.mantissa < (1L << 63))))
325         {
326             return true;
327         }
328         return false;
329     }
330 }
331 
332 ///
333 @nogc nothrow pure @safe unittest
334 {
335     assert(isNaN(float.init));
336     assert(isNaN(double.init));
337     assert(isNaN(real.init));
338 }
339 
340 /**
341  * Determines whether $(D_PARAM x) is a positive or negative infinity.
342  *
343  * Params:
344  *  F = Type of the floating point number.
345  *  x = Floating point number.
346  *
347  * Returns: $(D_KEYWORD true) if $(D_PARAM x) is infinity, $(D_KEYWORD false)
348  *          otherwise.
349  *
350  * See_Also: $(D_PSYMBOL isFinite).
351  */
352 bool isInfinity(F)(F x)
353 if (isFloatingPoint!F)
354 {
355     FloatBits!F bits;
356     bits.floating = abs(x);
357     static if (ieeePrecision!F == IEEEPrecision.single
358             || ieeePrecision!F == IEEEPrecision.double_)
359     {
360         return bits.integral == bits.expMask;
361     }
362     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
363     {
364         return (bits.exp == bits.expMask)
365             && ((bits.mantissa & bits.mantissaMask) == 0);
366     }
367 }
368 
369 ///
370 @nogc nothrow pure @safe unittest
371 {
372     assert(isInfinity(float.infinity));
373     assert(isInfinity(-float.infinity));
374     assert(isInfinity(double.infinity));
375     assert(isInfinity(-double.infinity));
376     assert(isInfinity(real.infinity));
377     assert(isInfinity(-real.infinity));
378 }
379 
380 /**
381  * Determines whether $(D_PARAM x) is a denormilized number or not.
382  *
383  * Denormalized number is a number between `0` and `1` that cannot be
384  * represented as
385  *
386  * <pre>
387  * m*2<sup>e</sup>
388  * </pre>
389  *
390  * where $(I m) is the mantissa and $(I e) is an exponent that fits into the
391  * exponent field of the type $(D_PARAM F).
392  *
393  * `0` is neither normalized nor denormalized.
394  *
395  * Params:
396  *  F = Type of the floating point number.
397  *  x = Floating point number.
398  *
399  * Returns: $(D_KEYWORD true) if $(D_PARAM x) is a denormilized number,
400  *          $(D_KEYWORD false) otherwise.
401  *
402  * See_Also: $(D_PSYMBOL isNormal).
403  */
404 bool isSubnormal(F)(F x)
405 if (isFloatingPoint!F)
406 {
407     FloatBits!F bits;
408     bits.floating = abs(x);
409     static if (ieeePrecision!F == IEEEPrecision.single)
410     {
411         return bits.integral < (1 << 23) && bits.integral > 0;
412     }
413     else static if (ieeePrecision!F == IEEEPrecision.double_)
414     {
415         return bits.integral < (1L << 52) && bits.integral > 0;
416     }
417     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
418     {
419         return bits.exp == 0 && bits.mantissa != 0;
420     }
421 }
422 
423 ///
424 @nogc nothrow pure @safe unittest
425 {
426     assert(!isSubnormal(0.0f));
427     assert(!isSubnormal(float.nan));
428     assert(!isSubnormal(float.infinity));
429     assert(!isSubnormal(0.3f));
430     assert(isSubnormal(5.87747e-38f / 10));
431 
432     assert(!isSubnormal(0.0));
433     assert(!isSubnormal(double.nan));
434     assert(!isSubnormal(double.infinity));
435     assert(!isSubnormal(1.4));
436     assert(isSubnormal(1.11254e-307 / 10));
437 
438     assert(!isSubnormal(0.0L));
439     assert(!isSubnormal(real.nan));
440     assert(!isSubnormal(real.infinity));
441 }
442 
443 /**
444  * Determines whether $(D_PARAM x) is a normilized number or not.
445  *
446  * Normalized number is a number that can be represented as
447  *
448  * <pre>
449  * m*2<sup>e</sup>
450  * </pre>
451  *
452  * where $(I m) is the mantissa and $(I e) is an exponent that fits into the
453  * exponent field of the type $(D_PARAM F).
454  *
455  * `0` is neither normalized nor denormalized.
456  *
457  * Params:
458  *  F = Type of the floating point number.
459  *  x = Floating point number.
460  *
461  * Returns: $(D_KEYWORD true) if $(D_PARAM x) is a normilized number,
462  *          $(D_KEYWORD false) otherwise.
463  *
464  * See_Also: $(D_PSYMBOL isSubnormal).
465  */
466 bool isNormal(F)(F x)
467 if (isFloatingPoint!F)
468 {
469     static if (ieeePrecision!F == IEEEPrecision.single
470             || ieeePrecision!F == IEEEPrecision.double_)
471     {
472         FloatBits!F bits;
473         bits.floating = x;
474         bits.integral &= bits.expMask;
475         return bits.integral != 0 && bits.integral != bits.expMask;
476     }
477     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
478     {
479         return classify(x) == FloatingPointClass.normal;
480     }
481 }
482 
483 ///
484 @nogc nothrow pure @safe unittest
485 {
486     assert(!isNormal(0.0f));
487     assert(!isNormal(float.nan));
488     assert(!isNormal(float.infinity));
489     assert(isNormal(0.3f));
490     assert(!isNormal(5.87747e-38f / 10));
491 
492     assert(!isNormal(0.0));
493     assert(!isNormal(double.nan));
494     assert(!isNormal(double.infinity));
495     assert(isNormal(1.4));
496     assert(!isNormal(1.11254e-307 / 10));
497 
498     assert(!isNormal(0.0L));
499     assert(!isNormal(real.nan));
500     assert(!isNormal(real.infinity));
501 }
502 
503 /**
504  * Determines whether the sign bit of $(D_PARAM x) is set or not.
505  *
506  * If the sign bit, $(D_PARAM x) is a negative number, otherwise positive.
507  *
508  * Params:
509  *  F = Type of the floating point number.
510  *  x = Floating point number.
511  *
512  * Returns: $(D_KEYWORD true) if the sign bit of $(D_PARAM x) is set,
513  *          $(D_KEYWORD false) otherwise.
514  */
515 bool signBit(F)(F x)
516 if (isFloatingPoint!F)
517 {
518     FloatBits!F bits;
519     bits.floating = x;
520     static if (ieeePrecision!F == IEEEPrecision.single)
521     {
522         return (bits.integral & (1 << 31)) != 0;
523     }
524     else static if (ieeePrecision!F == IEEEPrecision.double_)
525     {
526         return (bits.integral & (1L << 63)) != 0;
527     }
528     else static if (ieeePrecision!F == IEEEPrecision.doubleExtended)
529     {
530         return (bits.exp & (1 << 15)) != 0;
531     }
532 }
533 
534 ///
535 @nogc nothrow pure @safe unittest
536 {
537     assert(signBit(-1.0f));
538     assert(!signBit(1.0f));
539 
540     assert(signBit(-1.0));
541     assert(!signBit(1.0));
542 
543     assert(signBit(-1.0L));
544     assert(!signBit(1.0L));
545 }
546 
547 /**
548  * Computes $(D_PARAM x) to the power $(D_PARAM y) modulo $(D_PARAM z).
549  *
550  * Params:
551  *  I = Base type.
552  *  G = Exponent type.
553  *  H = Divisor type:
554  *  x = Base.
555  *  y = Exponent.
556  *  z = Divisor.
557  *
558  * Returns: Reminder of the division of $(D_PARAM x) to the power $(D_PARAM y)
559  *          by $(D_PARAM z).
560  *
561  * Precondition: $(D_INLINECODE z > 0)
562  */
563 H pow(I, G, H)(in auto ref I x, in auto ref G y, in auto ref H z)
564 if (isIntegral!I && isIntegral!G && isIntegral!H)
565 in
566 {
567     assert(z > 0, "Division by zero");
568 }
569 do
570 {
571     G mask = G.max / 2 + 1;
572     H result;
573 
574     if (y == 0)
575     {
576         return 1 % z;
577     }
578     else if (y == 1)
579     {
580         return x % z;
581     }
582     do
583     {
584         immutable bit = y & mask;
585         if (!result && bit)
586         {
587             result = x;
588             continue;
589         }
590 
591         result *= result;
592         if (bit)
593         {
594             result *= x;
595         }
596         result %= z;
597     }
598     while (mask >>= 1);
599 
600     return result;
601 }
602 
603 ///
604 @nogc nothrow pure @safe unittest
605 {
606     assert(pow(3, 5, 7) == 5);
607     assert(pow(2, 2, 1) == 0);
608     assert(pow(3, 3, 3) == 0);
609     assert(pow(7, 4, 2) == 1);
610     assert(pow(53, 0, 2) == 1);
611     assert(pow(53, 1, 3) == 2);
612     assert(pow(53, 2, 5) == 4);
613     assert(pow(0, 0, 5) == 1);
614     assert(pow(0, 5, 5) == 0);
615 }