// 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 $(MREF 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; 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 @nogcin{    // 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.");}