// 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, 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: * TABLE_SV = * * $0
Special Values
* SVH = $(TR $(TH $1) $(TH $2)) * SV = $(TR $(TD $1) $(TD $2)) * TH3 = $(TR $(TH $1) $(TH $2) $(TH $3)) * TD3 = $(TR $(TD $1) $(TD $2) $(TD $3)) * TABLE_DOMRG = * $(SVH Domain X, Range Y) $(SV $1, $2) *
* DOMAIN=$1 * RANGE=$1 * NAN = $(RED NAN) * SUP = $0 * GAMMA = Γ * THETA = θ * INTEGRAL = ∫ * INTEGRATE = $(BIG ∫$(SMALL $1)$2) * POWER = $1$2 * SUB = $1$2 * BIGSUM = $(BIG Σ $2$(SMALL $1)) * CHOOSE = $(BIG () $(SMALL $1)$(SMALL $2) $(BIG )) * PLUSMN = ± * INFIN = ∞ * PLUSMNINF = ±∞ * PI = π * LT = < * GT = > * SQRT = √ * HALF = ½ * * 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; static import core.math; static import core.stdc.math; static import core.stdc.fenv; import std.traits : CommonType, isFloatingPoint, isIntegral, isNumeric, isSigned, isUnsigned, Largest, Unqual; // @@@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; version (DigitalMars) { version = INLINE_YL2X; // x87 has opcodes for these } version (X86) version = X86_Any; version (X86_64) version = X86_Any; version (PPC) version = PPC_Any; version (PPC64) version = PPC_Any; version (MIPS32) version = MIPS_Any; version (MIPS64) version = MIPS_Any; version (AArch64) version = ARM_Any; version (ARM) version = ARM_Any; version (S390) version = IBMZ_Any; version (SPARC) version = SPARC_Any; version (SPARC64) version = SPARC_Any; version (SystemZ) version = IBMZ_Any; version (RISCV32) version = RISCV_Any; version (RISCV64) version = RISCV_Any; version (D_InlineAsm_X86) version = InlineAsm_X86_Any; version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any; version (InlineAsm_X86_Any) version = InlineAsm_X87; version (InlineAsm_X87) { static assert(real.mant_dig == 64); version (CRuntime_Microsoft) version = InlineAsm_X87_MSVC; } version (X86_64) version = StaticallyHaveSSE; version (X86) version (OSX) version = StaticallyHaveSSE; version (StaticallyHaveSSE) { private enum bool haveSSE = true; } else version (X86) { static import core.cpuid; private alias haveSSE = core.cpuid.sse; } /** Computes the value of a positive integer `x`, raised to the power `n`, modulo `m`. * * Params: * x = base * n = exponent * m = modulus * * Returns: * `x` to the power `n`, modulo `m`. * The return type is the largest of `x`'s and `m`'s type. * * The function requires that all values have unsigned types. */ Unqual!(Largest!(F, H)) powmod(F, G, H)(F x, G n, H m) if (isUnsigned!F && isUnsigned!G && isUnsigned!H) { import std.meta : AliasSeq; alias T = Unqual!(Largest!(F, H)); static if (T.sizeof <= 4) { alias DoubleT = AliasSeq!(void, ushort, uint, void, ulong)[T.sizeof]; } static T mulmod(T a, T b, T c) { static if (T.sizeof == 8) { static T addmod(T a, T b, T c) { b = c - b; if (a >= b) return a - b; else return c - b + a; } T result = 0, tmp; b %= c; while (a > 0) { if (a & 1) result = addmod(result, b, c); a >>= 1; b = addmod(b, b, c); } return result; } else { DoubleT result = cast(DoubleT) (cast(DoubleT) a * cast(DoubleT) b); return result % c; } } T base = x, result = 1, modulus = m; Unqual!G exponent = n; while (exponent > 0) { if (exponent & 1) result = mulmod(result, base, modulus); base = mulmod(base, base, modulus); exponent >>= 1; } return result; } /// @safe pure nothrow @nogc unittest { assert(powmod(1U, 10U, 3U) == 1); assert(powmod(3U, 2U, 6U) == 3); assert(powmod(5U, 5U, 15U) == 5); assert(powmod(2U, 3U, 5U) == 3); assert(powmod(2U, 4U, 5U) == 1); assert(powmod(2U, 5U, 5U) == 2); } @safe pure nothrow @nogc unittest { ulong a = 18446744073709551615u, b = 20u, c = 18446744073709551610u; assert(powmod(a, b, c) == 95367431640625u); a = 100; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 18223853583554725198u); a = 117; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 11493139548346411394u); a = 134; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 10979163786734356774u); a = 151; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 7023018419737782840u); a = 168; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 58082701842386811u); a = 185; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 17423478386299876798u); a = 202; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 5522733478579799075u); a = 219; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 15230218982491623487u); a = 236; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 5198328724976436000u); a = 0; b = 7919; c = 18446744073709551557u; assert(powmod(a, b, c) == 0); a = 123; b = 0; c = 18446744073709551557u; assert(powmod(a, b, c) == 1); immutable ulong a1 = 253, b1 = 7919, c1 = 18446744073709551557u; assert(powmod(a1, b1, c1) == 3883707345459248860u); uint x = 100 ,y = 7919, z = 1844674407u; assert(powmod(x, y, z) == 1613100340u); x = 134; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 734956622u); x = 151; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 1738696945u); x = 168; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 1247580927u); x = 185; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 1293855176u); x = 202; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 1566963682u); x = 219; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 181227807u); x = 236; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 217988321u); x = 253; y = 7919; z = 1844674407u; assert(powmod(x, y, z) == 1588843243u); x = 0; y = 7919; z = 184467u; assert(powmod(x, y, z) == 0); x = 123; y = 0; z = 1844674u; assert(powmod(x, y, z) == 1); immutable ubyte x1 = 117; immutable uint y1 = 7919; immutable uint z1 = 1844674407u; auto res = powmod(x1, y1, z1); assert(is(typeof(res) == uint)); assert(res == 9479781u); immutable ushort x2 = 123; immutable uint y2 = 203; immutable ubyte z2 = 113; auto res2 = powmod(x2, y2, z2); assert(is(typeof(res2) == ushort)); assert(res2 == 42u); } 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) { // 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; }