| // Written in the D programming language. |
| |
| /** |
| * Contains the elementary mathematical functions (powers, roots, |
| * and trigonometric functions), and low-level floating-point operations. |
| * Mathematical special functions are available in `std.mathspecial`. |
| * |
| $(SCRIPT inhibitQuickIndex = 1;) |
| |
| $(DIVC quickindex, |
| $(BOOKTABLE , |
| $(TR $(TH Category) $(TH Members) ) |
| $(TR $(TDNW $(SUBMODULE Constants, constants)) $(TD |
| $(SUBREF constants, E) |
| $(SUBREF constants, PI) |
| $(SUBREF constants, PI_2) |
| $(SUBREF constants, PI_4) |
| $(SUBREF constants, M_1_PI) |
| $(SUBREF constants, M_2_PI) |
| $(SUBREF constants, M_2_SQRTPI) |
| $(SUBREF constants, LN10) |
| $(SUBREF constants, LN2) |
| $(SUBREF constants, LOG2) |
| $(SUBREF constants, LOG2E) |
| $(SUBREF constants, LOG2T) |
| $(SUBREF constants, LOG10E) |
| $(SUBREF constants, SQRT2) |
| $(SUBREF constants, SQRT1_2) |
| )) |
| $(TR $(TDNW $(SUBMODULE Algebraic, algebraic)) $(TD |
| $(SUBREF algebraic, abs) |
| $(SUBREF algebraic, fabs) |
| $(SUBREF algebraic, sqrt) |
| $(SUBREF algebraic, cbrt) |
| $(SUBREF algebraic, hypot) |
| $(SUBREF algebraic, poly) |
| $(SUBREF algebraic, nextPow2) |
| $(SUBREF algebraic, truncPow2) |
| )) |
| $(TR $(TDNW $(SUBMODULE Trigonometry, trigonometry)) $(TD |
| $(SUBREF trigonometry, sin) |
| $(SUBREF trigonometry, cos) |
| $(SUBREF trigonometry, tan) |
| $(SUBREF trigonometry, asin) |
| $(SUBREF trigonometry, acos) |
| $(SUBREF trigonometry, atan) |
| $(SUBREF trigonometry, atan2) |
| $(SUBREF trigonometry, sinh) |
| $(SUBREF trigonometry, cosh) |
| $(SUBREF trigonometry, tanh) |
| $(SUBREF trigonometry, asinh) |
| $(SUBREF trigonometry, acosh) |
| $(SUBREF trigonometry, atanh) |
| )) |
| $(TR $(TDNW $(SUBMODULE Rounding, rounding)) $(TD |
| $(SUBREF rounding, ceil) |
| $(SUBREF rounding, floor) |
| $(SUBREF rounding, round) |
| $(SUBREF rounding, lround) |
| $(SUBREF rounding, trunc) |
| $(SUBREF rounding, rint) |
| $(SUBREF rounding, lrint) |
| $(SUBREF rounding, nearbyint) |
| $(SUBREF rounding, rndtol) |
| $(SUBREF rounding, quantize) |
| )) |
| $(TR $(TDNW $(SUBMODULE Exponentiation & Logarithms, exponential)) $(TD |
| $(SUBREF exponential, pow) |
| $(SUBREF exponential, powmod) |
| $(SUBREF exponential, exp) |
| $(SUBREF exponential, exp2) |
| $(SUBREF exponential, expm1) |
| $(SUBREF exponential, ldexp) |
| $(SUBREF exponential, frexp) |
| $(SUBREF exponential, log) |
| $(SUBREF exponential, log2) |
| $(SUBREF exponential, log10) |
| $(SUBREF exponential, logb) |
| $(SUBREF exponential, ilogb) |
| $(SUBREF exponential, log1p) |
| $(SUBREF exponential, scalbn) |
| )) |
| $(TR $(TDNW $(SUBMODULE Remainder, remainder)) $(TD |
| $(SUBREF remainder, fmod) |
| $(SUBREF remainder, modf) |
| $(SUBREF remainder, remainder) |
| $(SUBREF remainder, remquo) |
| )) |
| $(TR $(TDNW $(SUBMODULE Floating-point operations, operations)) $(TD |
| $(SUBREF operations, approxEqual) |
| $(SUBREF operations, feqrel) |
| $(SUBREF operations, fdim) |
| $(SUBREF operations, fmax) |
| $(SUBREF operations, fmin) |
| $(SUBREF operations, fma) |
| $(SUBREF operations, isClose) |
| $(SUBREF operations, nextDown) |
| $(SUBREF operations, nextUp) |
| $(SUBREF operations, nextafter) |
| $(SUBREF operations, NaN) |
| $(SUBREF operations, getNaNPayload) |
| $(SUBREF operations, cmp) |
| )) |
| $(TR $(TDNW $(SUBMODULE Introspection, traits)) $(TD |
| $(SUBREF traits, isFinite) |
| $(SUBREF traits, isIdentical) |
| $(SUBREF traits, isInfinity) |
| $(SUBREF traits, isNaN) |
| $(SUBREF traits, isNormal) |
| $(SUBREF traits, isSubnormal) |
| $(SUBREF traits, signbit) |
| $(SUBREF traits, sgn) |
| $(SUBREF traits, copysign) |
| $(SUBREF traits, isPowerOf2) |
| )) |
| $(TR $(TDNW $(SUBMODULE Hardware Control, hardware)) $(TD |
| $(SUBREF hardware, IeeeFlags) |
| $(SUBREF hardware, ieeeFlags) |
| $(SUBREF hardware, resetIeeeFlags) |
| $(SUBREF hardware, FloatingPointControl) |
| )) |
| ) |
| ) |
| |
| * The functionality closely follows the IEEE754-2008 standard for |
| * floating-point arithmetic, including the use of camelCase names rather |
| * than C99-style lower case names. All of these functions behave correctly |
| * when presented with an infinity or NaN. |
| * |
| * The following IEEE 'real' formats are currently supported: |
| * $(UL |
| * $(LI 64 bit Big-endian 'double' (eg PowerPC)) |
| * $(LI 128 bit Big-endian 'quadruple' (eg SPARC)) |
| * $(LI 64 bit Little-endian 'double' (eg x86-SSE2)) |
| * $(LI 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium)) |
| * $(LI 128 bit Little-endian 'quadruple' (not implemented on any known processor!)) |
| * $(LI Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support) |
| * ) |
| * Unlike C, there is no global 'errno' variable. Consequently, almost all of |
| * these functions are pure nothrow. |
| * |
| * Macros: |
| * SUBMODULE = $(MREF_ALTTEXT $1, std, math, $2) |
| * SUBREF = $(REF_ALTTEXT $(TT $2), $2, std, math, $1)$(NBSP) |
| * |
| * Copyright: Copyright The D Language Foundation 2000 - 2011. |
| * D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p, |
| * log2, floor, ceil and lrint functions are based on the CEPHES math library, |
| * which is Copyright (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT) |
| * and are incorporated herein by permission of the author. The author |
| * reserves the right to distribute this material elsewhere under different |
| * copying permissions. These modifications are distributed here under |
| * the following terms: |
| * License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). |
| * Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, |
| * Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger |
| * Source: $(PHOBOSSRC std/math/package.d) |
| */ |
| module std.math; |
| |
| public import std.math.algebraic; |
| public import std.math.constants; |
| public import std.math.exponential; |
| public import std.math.operations; |
| public import std.math.hardware; |
| public import std.math.remainder; |
| public import std.math.rounding; |
| public import std.math.traits; |
| public import std.math.trigonometry; |
| |
| // @@@DEPRECATED_2.102@@@ |
| // Note: Exposed accidentally, should be deprecated / removed |
| deprecated("std.meta.AliasSeq was unintentionally available from std.math " |
| ~ "and will be removed after 2.102. Please import std.meta instead") |
| public import std.meta : AliasSeq; |
| |
| package(std): // Not public yet |
| /* Return the value that lies halfway between x and y on the IEEE number line. |
| * |
| * Formally, the result is the arithmetic mean of the binary significands of x |
| * and y, multiplied by the geometric mean of the binary exponents of x and y. |
| * x and y must have the same sign, and must not be NaN. |
| * Note: this function is useful for ensuring O(log n) behaviour in algorithms |
| * involving a 'binary chop'. |
| * |
| * Special cases: |
| * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value |
| * is the arithmetic mean (x + y) / 2. |
| * If x and y are even powers of 2, the return value is the geometric mean, |
| * ieeeMean(x, y) = sqrt(x * y). |
| * |
| */ |
| T ieeeMean(T)(const T x, const T y) @trusted pure nothrow @nogc |
| in |
| { |
| // both x and y must have the same sign, and must not be NaN. |
| assert(signbit(x) == signbit(y)); |
| assert(x == x && y == y); |
| } |
| do |
| { |
| // Runtime behaviour for contract violation: |
| // If signs are opposite, or one is a NaN, return 0. |
| if (!((x >= 0 && y >= 0) || (x <= 0 && y <= 0))) return 0.0; |
| |
| // The implementation is simple: cast x and y to integers, |
| // average them (avoiding overflow), and cast the result back to a floating-point number. |
| |
| alias F = floatTraits!(T); |
| T u; |
| static if (F.realFormat == RealFormat.ieeeExtended || |
| F.realFormat == RealFormat.ieeeExtended53) |
| { |
| // There's slight additional complexity because they are actually |
| // 79-bit reals... |
| ushort *ue = cast(ushort *)&u; |
| ulong *ul = cast(ulong *)&u; |
| ushort *xe = cast(ushort *)&x; |
| ulong *xl = cast(ulong *)&x; |
| ushort *ye = cast(ushort *)&y; |
| ulong *yl = cast(ulong *)&y; |
| |
| // Ignore the useless implicit bit. (Bonus: this prevents overflows) |
| ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL); |
| |
| // @@@ BUG? @@@ |
| // Cast shouldn't be here |
| ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK) |
| + (ye[F.EXPPOS_SHORT] & F.EXPMASK)); |
| if (m & 0x8000_0000_0000_0000L) |
| { |
| ++e; |
| m &= 0x7FFF_FFFF_FFFF_FFFFL; |
| } |
| // Now do a multi-byte right shift |
| const uint c = e & 1; // carry |
| e >>= 1; |
| m >>>= 1; |
| if (c) |
| m |= 0x4000_0000_0000_0000L; // shift carry into significand |
| if (e) |
| *ul = m | 0x8000_0000_0000_0000L; // set implicit bit... |
| else |
| *ul = m; // ... unless exponent is 0 (subnormal or zero). |
| |
| ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit |
| } |
| else static if (F.realFormat == RealFormat.ieeeQuadruple) |
| { |
| // This would be trivial if 'ucent' were implemented... |
| ulong *ul = cast(ulong *)&u; |
| ulong *xl = cast(ulong *)&x; |
| ulong *yl = cast(ulong *)&y; |
| |
| // Multi-byte add, then multi-byte right shift. |
| import core.checkedint : addu; |
| bool carry; |
| ulong ml = addu(xl[MANTISSA_LSB], yl[MANTISSA_LSB], carry); |
| |
| ulong mh = carry + (xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) + |
| (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL); |
| |
| ul[MANTISSA_MSB] = (mh >>> 1) | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000); |
| ul[MANTISSA_LSB] = (ml >>> 1) | (mh & 1) << 63; |
| } |
| else static if (F.realFormat == RealFormat.ieeeDouble) |
| { |
| ulong *ul = cast(ulong *)&u; |
| ulong *xl = cast(ulong *)&x; |
| ulong *yl = cast(ulong *)&y; |
| ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) |
| + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1; |
| m |= ((*xl) & 0x8000_0000_0000_0000L); |
| *ul = m; |
| } |
| else static if (F.realFormat == RealFormat.ieeeSingle) |
| { |
| uint *ul = cast(uint *)&u; |
| uint *xl = cast(uint *)&x; |
| uint *yl = cast(uint *)&y; |
| uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1; |
| m |= ((*xl) & 0x8000_0000); |
| *ul = m; |
| } |
| else |
| { |
| assert(0, "Not implemented"); |
| } |
| return u; |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| assert(ieeeMean(-0.0,-1e-20)<0); |
| assert(ieeeMean(0.0,1e-20)>0); |
| |
| assert(ieeeMean(1.0L,4.0L)==2L); |
| assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013); |
| assert(ieeeMean(-1.0L,-4.0L)==-2L); |
| assert(ieeeMean(-1.0,-4.0)==-2); |
| assert(ieeeMean(-1.0f,-4.0f)==-2f); |
| assert(ieeeMean(-1.0,-2.0)==-1.5); |
| assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon)) |
| ==-1.5*(1+5*real.epsilon)); |
| assert(ieeeMean(0x1p60,0x1p-10)==0x1p25); |
| |
| static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) |
| { |
| assert(ieeeMean(1.0L,real.infinity)==0x1p8192L); |
| assert(ieeeMean(0.0L,real.infinity)==1.5); |
| } |
| assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal) |
| == 0.5*real.min_normal*(1-2*real.epsilon)); |
| } |
| |
| |
| // The following IEEE 'real' formats are currently supported. |
| version (LittleEndian) |
| { |
| static assert(real.mant_dig == 53 || real.mant_dig == 64 |
| || real.mant_dig == 113, |
| "Only 64-bit, 80-bit, and 128-bit reals"~ |
| " are supported for LittleEndian CPUs"); |
| } |
| else |
| { |
| static assert(real.mant_dig == 53 || real.mant_dig == 113, |
| "Only 64-bit and 128-bit reals are supported for BigEndian CPUs."); |
| } |
| |
| // Underlying format exposed through floatTraits |
| enum RealFormat |
| { |
| ieeeHalf, |
| ieeeSingle, |
| ieeeDouble, |
| ieeeExtended, // x87 80-bit real |
| ieeeExtended53, // x87 real rounded to precision of double. |
| ibmExtended, // IBM 128-bit extended |
| ieeeQuadruple, |
| } |
| |
| // Constants used for extracting the components of the representation. |
| // They supplement the built-in floating point properties. |
| template floatTraits(T) |
| { |
| import std.traits : Unqual; |
| |
| // EXPMASK is a ushort mask to select the exponent portion (without sign) |
| // EXPSHIFT is the number of bits the exponent is left-shifted by in its ushort |
| // EXPBIAS is the exponent bias - 1 (exp == EXPBIAS yields ×2^-1). |
| // EXPPOS_SHORT is the index of the exponent when represented as a ushort array. |
| // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array. |
| // RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal |
| enum Unqual!T RECIP_EPSILON = (1/T.epsilon); |
| static if (T.mant_dig == 24) |
| { |
| // Single precision float |
| enum ushort EXPMASK = 0x7F80; |
| enum ushort EXPSHIFT = 7; |
| enum ushort EXPBIAS = 0x3F00; |
| enum uint EXPMASK_INT = 0x7F80_0000; |
| enum uint MANTISSAMASK_INT = 0x007F_FFFF; |
| enum realFormat = RealFormat.ieeeSingle; |
| version (LittleEndian) |
| { |
| enum EXPPOS_SHORT = 1; |
| enum SIGNPOS_BYTE = 3; |
| } |
| else |
| { |
| enum EXPPOS_SHORT = 0; |
| enum SIGNPOS_BYTE = 0; |
| } |
| } |
| else static if (T.mant_dig == 53) |
| { |
| static if (T.sizeof == 8) |
| { |
| // Double precision float, or real == double |
| enum ushort EXPMASK = 0x7FF0; |
| enum ushort EXPSHIFT = 4; |
| enum ushort EXPBIAS = 0x3FE0; |
| enum uint EXPMASK_INT = 0x7FF0_0000; |
| enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only |
| enum realFormat = RealFormat.ieeeDouble; |
| version (LittleEndian) |
| { |
| enum EXPPOS_SHORT = 3; |
| enum SIGNPOS_BYTE = 7; |
| } |
| else |
| { |
| enum EXPPOS_SHORT = 0; |
| enum SIGNPOS_BYTE = 0; |
| } |
| } |
| else static if (T.sizeof == 12) |
| { |
| // Intel extended real80 rounded to double |
| enum ushort EXPMASK = 0x7FFF; |
| enum ushort EXPSHIFT = 0; |
| enum ushort EXPBIAS = 0x3FFE; |
| enum realFormat = RealFormat.ieeeExtended53; |
| version (LittleEndian) |
| { |
| enum EXPPOS_SHORT = 4; |
| enum SIGNPOS_BYTE = 9; |
| } |
| else |
| { |
| enum EXPPOS_SHORT = 0; |
| enum SIGNPOS_BYTE = 0; |
| } |
| } |
| else |
| static assert(false, "No traits support for " ~ T.stringof); |
| } |
| else static if (T.mant_dig == 64) |
| { |
| // Intel extended real80 |
| enum ushort EXPMASK = 0x7FFF; |
| enum ushort EXPSHIFT = 0; |
| enum ushort EXPBIAS = 0x3FFE; |
| enum realFormat = RealFormat.ieeeExtended; |
| version (LittleEndian) |
| { |
| enum EXPPOS_SHORT = 4; |
| enum SIGNPOS_BYTE = 9; |
| } |
| else |
| { |
| enum EXPPOS_SHORT = 0; |
| enum SIGNPOS_BYTE = 0; |
| } |
| } |
| else static if (T.mant_dig == 113) |
| { |
| // Quadruple precision float |
| enum ushort EXPMASK = 0x7FFF; |
| enum ushort EXPSHIFT = 0; |
| enum ushort EXPBIAS = 0x3FFE; |
| enum realFormat = RealFormat.ieeeQuadruple; |
| version (LittleEndian) |
| { |
| enum EXPPOS_SHORT = 7; |
| enum SIGNPOS_BYTE = 15; |
| } |
| else |
| { |
| enum EXPPOS_SHORT = 0; |
| enum SIGNPOS_BYTE = 0; |
| } |
| } |
| else static if (T.mant_dig == 106) |
| { |
| // IBM Extended doubledouble |
| enum ushort EXPMASK = 0x7FF0; |
| enum ushort EXPSHIFT = 4; |
| enum realFormat = RealFormat.ibmExtended; |
| |
| // For IBM doubledouble the larger magnitude double comes first. |
| // It's really a double[2] and arrays don't index differently |
| // between little and big-endian targets. |
| enum DOUBLEPAIR_MSB = 0; |
| enum DOUBLEPAIR_LSB = 1; |
| |
| // The exponent/sign byte is for most significant part. |
| version (LittleEndian) |
| { |
| enum EXPPOS_SHORT = 3; |
| enum SIGNPOS_BYTE = 7; |
| } |
| else |
| { |
| enum EXPPOS_SHORT = 0; |
| enum SIGNPOS_BYTE = 0; |
| } |
| } |
| else |
| static assert(false, "No traits support for " ~ T.stringof); |
| } |
| |
| // These apply to all floating-point types |
| version (LittleEndian) |
| { |
| enum MANTISSA_LSB = 0; |
| enum MANTISSA_MSB = 1; |
| } |
| else |
| { |
| enum MANTISSA_LSB = 1; |
| enum MANTISSA_MSB = 0; |
| } |