The OpenD Programming Language

1 // Written in the D programming language.
2 
3 /**
4  * Contains the elementary mathematical functions (powers, roots,
5  * and trigonometric functions), and low-level floating-point operations.
6  * Mathematical special functions are available in $(MREF std, mathspecial).
7  *
8 $(SCRIPT inhibitQuickIndex = 1;)
9 
10 $(DIVC quickindex,
11 $(BOOKTABLE ,
12 $(TR $(TH Category) $(TH Members) )
13 $(TR $(TDNW $(SUBMODULE Constants, constants)) $(TD
14     $(SUBREF constants, E)
15     $(SUBREF constants, PI)
16     $(SUBREF constants, PI_2)
17     $(SUBREF constants, PI_4)
18     $(SUBREF constants, M_1_PI)
19     $(SUBREF constants, M_2_PI)
20     $(SUBREF constants, M_2_SQRTPI)
21     $(SUBREF constants, LN10)
22     $(SUBREF constants, LN2)
23     $(SUBREF constants, LOG2)
24     $(SUBREF constants, LOG2E)
25     $(SUBREF constants, LOG2T)
26     $(SUBREF constants, LOG10E)
27     $(SUBREF constants, SQRT2)
28     $(SUBREF constants, SQRT1_2)
29 ))
30 $(TR $(TDNW $(SUBMODULE Algebraic, algebraic)) $(TD
31     $(SUBREF algebraic, abs)
32     $(SUBREF algebraic, fabs)
33     $(SUBREF algebraic, sqrt)
34     $(SUBREF algebraic, cbrt)
35     $(SUBREF algebraic, hypot)
36     $(SUBREF algebraic, poly)
37     $(SUBREF algebraic, nextPow2)
38     $(SUBREF algebraic, truncPow2)
39 ))
40 $(TR $(TDNW $(SUBMODULE Trigonometry, trigonometry)) $(TD
41     $(SUBREF trigonometry, sin)
42     $(SUBREF trigonometry, cos)
43     $(SUBREF trigonometry, tan)
44     $(SUBREF trigonometry, asin)
45     $(SUBREF trigonometry, acos)
46     $(SUBREF trigonometry, atan)
47     $(SUBREF trigonometry, atan2)
48     $(SUBREF trigonometry, sinh)
49     $(SUBREF trigonometry, cosh)
50     $(SUBREF trigonometry, tanh)
51     $(SUBREF trigonometry, asinh)
52     $(SUBREF trigonometry, acosh)
53     $(SUBREF trigonometry, atanh)
54 ))
55 $(TR $(TDNW $(SUBMODULE Rounding, rounding)) $(TD
56     $(SUBREF rounding, ceil)
57     $(SUBREF rounding, floor)
58     $(SUBREF rounding, round)
59     $(SUBREF rounding, lround)
60     $(SUBREF rounding, trunc)
61     $(SUBREF rounding, rint)
62     $(SUBREF rounding, lrint)
63     $(SUBREF rounding, nearbyint)
64     $(SUBREF rounding, rndtol)
65     $(SUBREF rounding, quantize)
66 ))
67 $(TR $(TDNW $(SUBMODULE Exponentiation & Logarithms, exponential)) $(TD
68     $(SUBREF exponential, pow)
69     $(SUBREF exponential, powmod)
70     $(SUBREF exponential, exp)
71     $(SUBREF exponential, exp2)
72     $(SUBREF exponential, expm1)
73     $(SUBREF exponential, ldexp)
74     $(SUBREF exponential, frexp)
75     $(SUBREF exponential, log)
76     $(SUBREF exponential, log2)
77     $(SUBREF exponential, log10)
78     $(SUBREF exponential, logb)
79     $(SUBREF exponential, ilogb)
80     $(SUBREF exponential, log1p)
81     $(SUBREF exponential, scalbn)
82 ))
83 $(TR $(TDNW $(SUBMODULE Remainder, remainder)) $(TD
84     $(SUBREF remainder, fmod)
85     $(SUBREF remainder, modf)
86     $(SUBREF remainder, remainder)
87     $(SUBREF remainder, remquo)
88 ))
89 $(TR $(TDNW $(SUBMODULE Floating-point operations, operations)) $(TD
90     $(SUBREF operations, approxEqual)
91     $(SUBREF operations, feqrel)
92     $(SUBREF operations, fdim)
93     $(SUBREF operations, fmax)
94     $(SUBREF operations, fmin)
95     $(SUBREF operations, fma)
96     $(SUBREF operations, isClose)
97     $(SUBREF operations, nextDown)
98     $(SUBREF operations, nextUp)
99     $(SUBREF operations, nextafter)
100     $(SUBREF operations, NaN)
101     $(SUBREF operations, getNaNPayload)
102     $(SUBREF operations, cmp)
103 ))
104 $(TR $(TDNW $(SUBMODULE Introspection, traits)) $(TD
105     $(SUBREF traits, isFinite)
106     $(SUBREF traits, isIdentical)
107     $(SUBREF traits, isInfinity)
108     $(SUBREF traits, isNaN)
109     $(SUBREF traits, isNormal)
110     $(SUBREF traits, isSubnormal)
111     $(SUBREF traits, signbit)
112     $(SUBREF traits, sgn)
113     $(SUBREF traits, copysign)
114     $(SUBREF traits, isPowerOf2)
115 ))
116 $(TR $(TDNW $(SUBMODULE Hardware Control, hardware)) $(TD
117     $(SUBREF hardware, IeeeFlags)
118     $(SUBREF hardware, ieeeFlags)
119     $(SUBREF hardware, resetIeeeFlags)
120     $(SUBREF hardware, FloatingPointControl)
121 ))
122 )
123 )
124 
125  * The functionality closely follows the IEEE754-2008 standard for
126  * floating-point arithmetic, including the use of camelCase names rather
127  * than C99-style lower case names. All of these functions behave correctly
128  * when presented with an infinity or NaN.
129  *
130  * The following IEEE 'real' formats are currently supported:
131  * $(UL
132  * $(LI 64 bit Big-endian  'double' (eg PowerPC))
133  * $(LI 128 bit Big-endian 'quadruple' (eg SPARC))
134  * $(LI 64 bit Little-endian 'double' (eg x86-SSE2))
135  * $(LI 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium))
136  * $(LI 128 bit Little-endian 'quadruple' (not implemented on any known processor!))
137  * $(LI Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support)
138  * )
139  * Unlike C, there is no global 'errno' variable. Consequently, almost all of
140  * these functions are pure nothrow.
141  *
142  * Macros:
143  *      SUBMODULE = $(MREF_ALTTEXT $1, std, math, $2)
144  *      SUBREF = $(REF_ALTTEXT $(TT $2), $2, std, math, $1)$(NBSP)
145  *
146  * Copyright: Copyright The D Language Foundation 2000 - 2011.
147  *            D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p,
148  *            log2, floor, ceil and lrint functions are based on the CEPHES math library,
149  *            which is Copyright (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT)
150  *            and are incorporated herein by permission of the author.  The author
151  *            reserves the right to distribute this material elsewhere under different
152  *            copying permissions.  These modifications are distributed here under
153  *            the following terms:
154  * License:   $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
155  * Authors:   $(HTTP digitalmars.com, Walter Bright), Don Clugston,
156  *            Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger
157  * Source: $(PHOBOSSRC std/math/package.d)
158  */
159 module std.math;
160 
161 public import std.math.algebraic;
162 public import std.math.constants;
163 public import std.math.exponential;
164 public import std.math.operations;
165 public import std.math.hardware;
166 public import std.math.remainder;
167 public import std.math.rounding;
168 public import std.math.traits;
169 public import std.math.trigonometry;
170 
171 package(std): // Not public yet
172 /* Return the value that lies halfway between x and y on the IEEE number line.
173  *
174  * Formally, the result is the arithmetic mean of the binary significands of x
175  * and y, multiplied by the geometric mean of the binary exponents of x and y.
176  * x and y must have the same sign, and must not be NaN.
177  * Note: this function is useful for ensuring O(log n) behaviour in algorithms
178  * involving a 'binary chop'.
179  *
180  * Special cases:
181  * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
182  * is the arithmetic mean (x + y) / 2.
183  * If x and y are even powers of 2, the return value is the geometric mean,
184  *   ieeeMean(x, y) = sqrt(x * y).
185  *
186  */
187 T ieeeMean(T)(const T x, const T y)  @trusted pure nothrow @nogc
188 in
189 {
190     // both x and y must have the same sign, and must not be NaN.
191     assert(signbit(x) == signbit(y));
192     assert(x == x && y == y);
193 }
194 do
195 {
196     // Runtime behaviour for contract violation:
197     // If signs are opposite, or one is a NaN, return 0.
198     if (!((x >= 0 && y >= 0) || (x <= 0 && y <= 0))) return 0.0;
199 
200     // The implementation is simple: cast x and y to integers,
201     // average them (avoiding overflow), and cast the result back to a floating-point number.
202 
203     alias F = floatTraits!(T);
204     T u;
205     static if (F.realFormat == RealFormat.ieeeExtended ||
206                F.realFormat == RealFormat.ieeeExtended53)
207     {
208         // There's slight additional complexity because they are actually
209         // 79-bit reals...
210         ushort *ue = cast(ushort *)&u;
211         ulong *ul = cast(ulong *)&u;
212         ushort *xe = cast(ushort *)&x;
213         ulong *xl = cast(ulong *)&x;
214         ushort *ye = cast(ushort *)&y;
215         ulong *yl = cast(ulong *)&y;
216 
217         // Ignore the useless implicit bit. (Bonus: this prevents overflows)
218         ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL);
219 
220         // @@@ BUG? @@@
221         // Cast shouldn't be here
222         ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK)
223                                  + (ye[F.EXPPOS_SHORT] & F.EXPMASK));
224         if (m & 0x8000_0000_0000_0000L)
225         {
226             ++e;
227             m &= 0x7FFF_FFFF_FFFF_FFFFL;
228         }
229         // Now do a multi-byte right shift
230         const uint c = e & 1; // carry
231         e >>= 1;
232         m >>>= 1;
233         if (c)
234             m |= 0x4000_0000_0000_0000L; // shift carry into significand
235         if (e)
236             *ul = m | 0x8000_0000_0000_0000L; // set implicit bit...
237         else
238             *ul = m; // ... unless exponent is 0 (subnormal or zero).
239 
240         ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit
241     }
242     else static if (F.realFormat == RealFormat.ieeeQuadruple)
243     {
244         // This would be trivial if 'ucent' were implemented...
245         ulong *ul = cast(ulong *)&u;
246         ulong *xl = cast(ulong *)&x;
247         ulong *yl = cast(ulong *)&y;
248 
249         // Multi-byte add, then multi-byte right shift.
250         import core.checkedint : addu;
251         bool carry;
252         ulong ml = addu(xl[MANTISSA_LSB], yl[MANTISSA_LSB], carry);
253 
254         ulong mh = carry + (xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) +
255             (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL);
256 
257         ul[MANTISSA_MSB] = (mh >>> 1) | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
258         ul[MANTISSA_LSB] = (ml >>> 1) | (mh & 1) << 63;
259     }
260     else static if (F.realFormat == RealFormat.ieeeDouble)
261     {
262         ulong *ul = cast(ulong *)&u;
263         ulong *xl = cast(ulong *)&x;
264         ulong *yl = cast(ulong *)&y;
265         ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL)
266                    + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
267         m |= ((*xl) & 0x8000_0000_0000_0000L);
268         *ul = m;
269     }
270     else static if (F.realFormat == RealFormat.ieeeSingle)
271     {
272         uint *ul = cast(uint *)&u;
273         uint *xl = cast(uint *)&x;
274         uint *yl = cast(uint *)&y;
275         uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1;
276         m |= ((*xl) & 0x8000_0000);
277         *ul = m;
278     }
279     else
280     {
281         assert(0, "Not implemented");
282     }
283     return u;
284 }
285 
286 @safe pure nothrow @nogc unittest
287 {
288     assert(ieeeMean(-0.0,-1e-20)<0);
289     assert(ieeeMean(0.0,1e-20)>0);
290 
291     assert(ieeeMean(1.0L,4.0L)==2L);
292     assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013);
293     assert(ieeeMean(-1.0L,-4.0L)==-2L);
294     assert(ieeeMean(-1.0,-4.0)==-2);
295     assert(ieeeMean(-1.0f,-4.0f)==-2f);
296     assert(ieeeMean(-1.0,-2.0)==-1.5);
297     assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon))
298                  ==-1.5*(1+5*real.epsilon));
299     assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
300 
301     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
302     {
303       assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
304       assert(ieeeMean(0.0L,real.infinity)==1.5);
305     }
306     assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal)
307            == 0.5*real.min_normal*(1-2*real.epsilon));
308 }
309 
310 
311 // The following IEEE 'real' formats are currently supported.
312 version (LittleEndian)
313 {
314     static assert(real.mant_dig == 53 || real.mant_dig == 64
315                || real.mant_dig == 113,
316       "Only 64-bit, 80-bit, and 128-bit reals"~
317       " are supported for LittleEndian CPUs");
318 }
319 else
320 {
321     static assert(real.mant_dig == 53 || real.mant_dig == 113,
322     "Only 64-bit and 128-bit reals are supported for BigEndian CPUs.");
323 }
324 
325 // Underlying format exposed through floatTraits
326 enum RealFormat
327 {
328     ieeeHalf,
329     ieeeSingle,
330     ieeeDouble,
331     ieeeExtended,   // x87 80-bit real
332     ieeeExtended53, // x87 real rounded to precision of double.
333     ibmExtended,    // IBM 128-bit extended
334     ieeeQuadruple,
335 }
336 
337 // Constants used for extracting the components of the representation.
338 // They supplement the built-in floating point properties.
339 template floatTraits(T)
340 {
341     import std.traits : Unqual;
342 
343     // EXPMASK is a ushort mask to select the exponent portion (without sign)
344     // EXPSHIFT is the number of bits the exponent is left-shifted by in its ushort
345     // EXPBIAS is the exponent bias - 1 (exp == EXPBIAS yields ×2^-1).
346     // EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
347     // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
348     // RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal
349     enum Unqual!T RECIP_EPSILON = (1/T.epsilon);
350     static if (T.mant_dig == 24)
351     {
352         // Single precision float
353         enum ushort EXPMASK = 0x7F80;
354         enum ushort EXPSHIFT = 7;
355         enum ushort EXPBIAS = 0x3F00;
356         enum uint EXPMASK_INT = 0x7F80_0000;
357         enum uint MANTISSAMASK_INT = 0x007F_FFFF;
358         enum realFormat = RealFormat.ieeeSingle;
359         version (LittleEndian)
360         {
361             enum EXPPOS_SHORT = 1;
362             enum SIGNPOS_BYTE = 3;
363         }
364         else
365         {
366             enum EXPPOS_SHORT = 0;
367             enum SIGNPOS_BYTE = 0;
368         }
369     }
370     else static if (T.mant_dig == 53)
371     {
372         static if (T.sizeof == 8)
373         {
374             // Double precision float, or real == double
375             enum ushort EXPMASK = 0x7FF0;
376             enum ushort EXPSHIFT = 4;
377             enum ushort EXPBIAS = 0x3FE0;
378             enum uint EXPMASK_INT = 0x7FF0_0000;
379             enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
380             enum ulong MANTISSAMASK_LONG = 0x000F_FFFF_FFFF_FFFF;
381             enum realFormat = RealFormat.ieeeDouble;
382             version (LittleEndian)
383             {
384                 enum EXPPOS_SHORT = 3;
385                 enum SIGNPOS_BYTE = 7;
386             }
387             else
388             {
389                 enum EXPPOS_SHORT = 0;
390                 enum SIGNPOS_BYTE = 0;
391             }
392         }
393         else static if (T.sizeof == 12)
394         {
395             // Intel extended real80 rounded to double
396             enum ushort EXPMASK = 0x7FFF;
397             enum ushort EXPSHIFT = 0;
398             enum ushort EXPBIAS = 0x3FFE;
399             enum realFormat = RealFormat.ieeeExtended53;
400             version (LittleEndian)
401             {
402                 enum EXPPOS_SHORT = 4;
403                 enum SIGNPOS_BYTE = 9;
404             }
405             else
406             {
407                 enum EXPPOS_SHORT = 0;
408                 enum SIGNPOS_BYTE = 0;
409             }
410         }
411         else
412             static assert(false, "No traits support for " ~ T.stringof);
413     }
414     else static if (T.mant_dig == 64)
415     {
416         // Intel extended real80
417         enum ushort EXPMASK = 0x7FFF;
418         enum ushort EXPSHIFT = 0;
419         enum ushort EXPBIAS = 0x3FFE;
420         enum realFormat = RealFormat.ieeeExtended;
421         version (LittleEndian)
422         {
423             enum EXPPOS_SHORT = 4;
424             enum SIGNPOS_BYTE = 9;
425         }
426         else
427         {
428             enum EXPPOS_SHORT = 0;
429             enum SIGNPOS_BYTE = 0;
430         }
431     }
432     else static if (T.mant_dig == 113)
433     {
434         // Quadruple precision float
435         enum ushort EXPMASK = 0x7FFF;
436         enum ushort EXPSHIFT = 0;
437         enum ushort EXPBIAS = 0x3FFE;
438         enum realFormat = RealFormat.ieeeQuadruple;
439         version (LittleEndian)
440         {
441             enum EXPPOS_SHORT = 7;
442             enum SIGNPOS_BYTE = 15;
443         }
444         else
445         {
446             enum EXPPOS_SHORT = 0;
447             enum SIGNPOS_BYTE = 0;
448         }
449     }
450     else static if (T.mant_dig == 106)
451     {
452         // IBM Extended doubledouble
453         enum ushort EXPMASK = 0x7FF0;
454         enum ushort EXPSHIFT = 4;
455         enum realFormat = RealFormat.ibmExtended;
456 
457         // For IBM doubledouble the larger magnitude double comes first.
458         // It's really a double[2] and arrays don't index differently
459         // between little and big-endian targets.
460         enum DOUBLEPAIR_MSB = 0;
461         enum DOUBLEPAIR_LSB = 1;
462 
463         // The exponent/sign byte is for most significant part.
464         version (LittleEndian)
465         {
466             enum EXPPOS_SHORT = 3;
467             enum SIGNPOS_BYTE = 7;
468         }
469         else
470         {
471             enum EXPPOS_SHORT = 0;
472             enum SIGNPOS_BYTE = 0;
473         }
474     }
475     else
476         static assert(false, "No traits support for " ~ T.stringof);
477 }
478 
479 // These apply to all floating-point types
480 version (LittleEndian)
481 {
482     enum MANTISSA_LSB = 0;
483     enum MANTISSA_MSB = 1;
484 }
485 else
486 {
487     enum MANTISSA_LSB = 1;
488     enum MANTISSA_MSB = 0;
489 }