| // 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 $(D std.mathspecial). |
| * |
| $(SCRIPT inhibitQuickIndex = 1;) |
| |
| $(DIVC quickindex, |
| $(BOOKTABLE , |
| $(TR $(TH Category) $(TH Members) ) |
| $(TR $(TDNW Constants) $(TD |
| $(MYREF E) $(MYREF PI) $(MYREF PI_2) $(MYREF PI_4) $(MYREF M_1_PI) |
| $(MYREF M_2_PI) $(MYREF M_2_SQRTPI) $(MYREF LN10) $(MYREF LN2) |
| $(MYREF LOG2) $(MYREF LOG2E) $(MYREF LOG2T) $(MYREF LOG10E) |
| $(MYREF SQRT2) $(MYREF SQRT1_2) |
| )) |
| $(TR $(TDNW Classics) $(TD |
| $(MYREF abs) $(MYREF fabs) $(MYREF sqrt) $(MYREF cbrt) $(MYREF hypot) |
| $(MYREF poly) $(MYREF nextPow2) $(MYREF truncPow2) |
| )) |
| $(TR $(TDNW Trigonometry) $(TD |
| $(MYREF sin) $(MYREF cos) $(MYREF tan) $(MYREF asin) $(MYREF acos) |
| $(MYREF atan) $(MYREF atan2) $(MYREF sinh) $(MYREF cosh) $(MYREF tanh) |
| $(MYREF asinh) $(MYREF acosh) $(MYREF atanh) $(MYREF expi) |
| )) |
| $(TR $(TDNW Rounding) $(TD |
| $(MYREF ceil) $(MYREF floor) $(MYREF round) $(MYREF lround) |
| $(MYREF trunc) $(MYREF rint) $(MYREF lrint) $(MYREF nearbyint) |
| $(MYREF rndtol) $(MYREF quantize) |
| )) |
| $(TR $(TDNW Exponentiation & Logarithms) $(TD |
| $(MYREF pow) $(MYREF exp) $(MYREF exp2) $(MYREF expm1) $(MYREF ldexp) |
| $(MYREF frexp) $(MYREF log) $(MYREF log2) $(MYREF log10) $(MYREF logb) |
| $(MYREF ilogb) $(MYREF log1p) $(MYREF scalbn) |
| )) |
| $(TR $(TDNW Modulus) $(TD |
| $(MYREF fmod) $(MYREF modf) $(MYREF remainder) |
| )) |
| $(TR $(TDNW Floating-point operations) $(TD |
| $(MYREF approxEqual) $(MYREF feqrel) $(MYREF fdim) $(MYREF fmax) |
| $(MYREF fmin) $(MYREF fma) $(MYREF nextDown) $(MYREF nextUp) |
| $(MYREF nextafter) $(MYREF NaN) $(MYREF getNaNPayload) |
| $(MYREF cmp) |
| )) |
| $(TR $(TDNW Introspection) $(TD |
| $(MYREF isFinite) $(MYREF isIdentical) $(MYREF isInfinity) $(MYREF isNaN) |
| $(MYREF isNormal) $(MYREF isSubnormal) $(MYREF signbit) $(MYREF sgn) |
| $(MYREF copysign) $(MYREF isPowerOf2) |
| )) |
| $(TR $(TDNW Complex Numbers) $(TD |
| $(MYREF abs) $(MYREF conj) $(MYREF sin) $(MYREF cos) $(MYREF expi) |
| )) |
| $(TR $(TDNW Hardware Control) $(TD |
| $(MYREF IeeeFlags) $(MYREF 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. |
| * |
| * Status: |
| * The semantics and names of feqrel and approxEqual will be revised. |
| * |
| * Macros: |
| * TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> |
| * <caption>Special Values</caption> |
| * $0</table> |
| * 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 = <table border="1" cellpadding="4" cellspacing="0"> |
| * $(SVH Domain X, Range Y) |
| $(SV $1, $2) |
| * </table> |
| * DOMAIN=$1 |
| * RANGE=$1 |
| |
| * NAN = $(RED NAN) |
| * SUP = <span style="vertical-align:super;font-size:smaller">$0</span> |
| * GAMMA = Γ |
| * THETA = θ |
| * INTEGRAL = ∫ |
| * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) |
| * POWER = $1<sup>$2</sup> |
| * SUB = $1<sub>$2</sub> |
| * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>) |
| * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG )) |
| * PLUSMN = ± |
| * INFIN = ∞ |
| * PLUSMNINF = ±∞ |
| * PI = π |
| * LT = < |
| * GT = > |
| * SQRT = √ |
| * HALF = ½ |
| * |
| * Copyright: Copyright Digital Mars 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.d) |
| */ |
| |
| /* NOTE: This file has been patched from the original DMD distribution to |
| * work with the GDC compiler. |
| */ |
| module std.math; |
| |
| version (Win64) |
| { |
| version (D_InlineAsm_X86_64) |
| version = Win64_DMD_InlineAsm; |
| } |
| |
| static import core.math; |
| static import core.stdc.math; |
| static import core.stdc.fenv; |
| import std.traits; // CommonType, isFloatingPoint, isIntegral, isSigned, isUnsigned, Largest, Unqual |
| |
| version (LDC) |
| { |
| import ldc.intrinsics; |
| } |
| |
| 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; |
| } |
| |
| version (D_SoftFloat) |
| { |
| // Some soft float implementations may support IEEE floating flags. |
| // The implementation here supports hardware flags only and is so currently |
| // only available for supported targets. |
| } |
| else version (X86_Any) version = IeeeFlagsSupport; |
| else version (PPC_Any) version = IeeeFlagsSupport; |
| else version (RISCV_Any) version = IeeeFlagsSupport; |
| else version (MIPS_Any) version = IeeeFlagsSupport; |
| else version (ARM_Any) version = IeeeFlagsSupport; |
| |
| // Struct FloatingPointControl is only available if hardware FP units are available. |
| version (D_HardFloat) |
| { |
| // FloatingPointControl.clearExceptions() depends on version IeeeFlagsSupport |
| version (IeeeFlagsSupport) version = FloatingPointControlSupport; |
| } |
| |
| version (GNU) |
| { |
| // The compiler can unexpectedly rearrange floating point operations and |
| // access to the floating point status flags when optimizing. This means |
| // ieeeFlags tests cannot be reliably checked in optimized code. |
| // See https://github.com/ldc-developers/ldc/issues/888 |
| } |
| else |
| { |
| version = IeeeFlagsUnittest; |
| version = FloatingPointControlUnittest; |
| } |
| |
| version (unittest) |
| { |
| import core.stdc.stdio; // : sprintf; |
| |
| static if (real.sizeof > double.sizeof) |
| enum uint useDigits = 16; |
| else |
| enum uint useDigits = 15; |
| |
| /****************************************** |
| * Compare floating point numbers to n decimal digits of precision. |
| * Returns: |
| * 1 match |
| * 0 nomatch |
| */ |
| |
| private bool equalsDigit(real x, real y, uint ndigits) |
| { |
| if (signbit(x) != signbit(y)) |
| return 0; |
| |
| if (isInfinity(x) && isInfinity(y)) |
| return 1; |
| if (isInfinity(x) || isInfinity(y)) |
| return 0; |
| |
| if (isNaN(x) && isNaN(y)) |
| return 1; |
| if (isNaN(x) || isNaN(y)) |
| return 0; |
| |
| char[30] bufx; |
| char[30] bufy; |
| assert(ndigits < bufx.length); |
| |
| int ix; |
| int iy; |
| version (CRuntime_Microsoft) |
| alias real_t = double; |
| else |
| alias real_t = real; |
| ix = sprintf(bufx.ptr, is(real_t == real) ? "%.*Lg" : "%.*g", ndigits, cast(real_t) x); |
| iy = sprintf(bufy.ptr, is(real_t == real) ? "%.*Lg" : "%.*g", ndigits, cast(real_t) y); |
| assert(ix < bufx.length && ix > 0); |
| assert(ix < bufy.length && ix > 0); |
| |
| return bufx[0 .. ix] == bufy[0 .. iy]; |
| } |
| } |
| |
| |
| |
| package: |
| // 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 == 106 |
| || real.mant_dig == 113, |
| "Only 64-bit and 128-bit reals are supported for BigEndian CPUs."~ |
| " double-double reals have partial support"); |
| } |
| |
| // 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 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; |
| } |
| |
| // Common code for math implementations. |
| |
| // Helper for floor/ceil |
| T floorImpl(T)(const T x) @trusted pure nothrow @nogc |
| { |
| alias F = floatTraits!(T); |
| // Take care not to trigger library calls from the compiler, |
| // while ensuring that we don't get defeated by some optimizers. |
| union floatBits |
| { |
| T rv; |
| ushort[T.sizeof/2] vu; |
| |
| // Other kinds of extractors for real formats. |
| static if (F.realFormat == RealFormat.ieeeSingle) |
| int vi; |
| } |
| floatBits y = void; |
| y.rv = x; |
| |
| // Find the exponent (power of 2) |
| // Do this by shifting the raw value so that the exponent lies in the low bits, |
| // then mask out the sign bit, and subtract the bias. |
| static if (F.realFormat == RealFormat.ieeeSingle) |
| { |
| int exp = ((y.vi >> (T.mant_dig - 1)) & 0xff) - 0x7f; |
| } |
| else static if (F.realFormat == RealFormat.ieeeDouble) |
| { |
| int exp = ((y.vu[F.EXPPOS_SHORT] >> 4) & 0x7ff) - 0x3ff; |
| |
| version (LittleEndian) |
| int pos = 0; |
| else |
| int pos = 3; |
| } |
| else static if (F.realFormat == RealFormat.ieeeExtended || |
| F.realFormat == RealFormat.ieeeExtended53) |
| { |
| int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; |
| |
| version (LittleEndian) |
| int pos = 0; |
| else |
| int pos = 4; |
| } |
| else static if (F.realFormat == RealFormat.ieeeQuadruple) |
| { |
| int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; |
| |
| version (LittleEndian) |
| int pos = 0; |
| else |
| int pos = 7; |
| } |
| else |
| static assert(false, "Not implemented for this architecture"); |
| |
| if (exp < 0) |
| { |
| if (x < 0.0) |
| return -1.0; |
| else |
| return 0.0; |
| } |
| |
| static if (F.realFormat == RealFormat.ieeeSingle) |
| { |
| if (exp < (T.mant_dig - 1)) |
| { |
| // Clear all bits representing the fraction part. |
| const uint fraction_mask = F.MANTISSAMASK_INT >> exp; |
| |
| if ((y.vi & fraction_mask) != 0) |
| { |
| // If 'x' is negative, then first substract 1.0 from the value. |
| if (y.vi < 0) |
| y.vi += 0x00800000 >> exp; |
| y.vi &= ~fraction_mask; |
| } |
| } |
| } |
| else |
| { |
| static if (F.realFormat == RealFormat.ieeeExtended53) |
| exp = (T.mant_dig + 11 - 1) - exp; // mant_dig is really 64 |
| else |
| exp = (T.mant_dig - 1) - exp; |
| |
| // Zero 16 bits at a time. |
| while (exp >= 16) |
| { |
| version (LittleEndian) |
| y.vu[pos++] = 0; |
| else |
| y.vu[pos--] = 0; |
| exp -= 16; |
| } |
| |
| // Clear the remaining bits. |
| if (exp > 0) |
| y.vu[pos] &= 0xffff ^ ((1 << exp) - 1); |
| |
| if ((x < 0.0) && (x != y.rv)) |
| y.rv -= 1.0; |
| } |
| |
| return y.rv; |
| } |
| |
| public: |
| |
| // Values obtained from Wolfram Alpha. 116 bits ought to be enough for anybody. |
| // Wolfram Alpha LLC. 2011. Wolfram|Alpha. http://www.wolframalpha.com/input/?i=e+in+base+16 (access July 6, 2011). |
| enum real E = 0x1.5bf0a8b1457695355fb8ac404e7a8p+1L; /** e = 2.718281... */ |
| enum real LOG2T = 0x1.a934f0979a3715fc9257edfe9b5fbp+1L; /** $(SUB log, 2)10 = 3.321928... */ |
| enum real LOG2E = 0x1.71547652b82fe1777d0ffda0d23a8p+0L; /** $(SUB log, 2)e = 1.442695... */ |
| enum real LOG2 = 0x1.34413509f79fef311f12b35816f92p-2L; /** $(SUB log, 10)2 = 0.301029... */ |
| enum real LOG10E = 0x1.bcb7b1526e50e32a6ab7555f5a67cp-2L; /** $(SUB log, 10)e = 0.434294... */ |
| enum real LN2 = 0x1.62e42fefa39ef35793c7673007e5fp-1L; /** ln 2 = 0.693147... */ |
| enum real LN10 = 0x1.26bb1bbb5551582dd4adac5705a61p+1L; /** ln 10 = 2.302585... */ |
| enum real PI = 0x1.921fb54442d18469898cc51701b84p+1L; /** $(_PI) = 3.141592... */ |
| enum real PI_2 = PI/2; /** $(PI) / 2 = 1.570796... */ |
| enum real PI_4 = PI/4; /** $(PI) / 4 = 0.785398... */ |
| enum real M_1_PI = 0x1.45f306dc9c882a53f84eafa3ea69cp-2L; /** 1 / $(PI) = 0.318309... */ |
| enum real M_2_PI = 2*M_1_PI; /** 2 / $(PI) = 0.636619... */ |
| enum real M_2_SQRTPI = 0x1.20dd750429b6d11ae3a914fed7fd8p+0L; /** 2 / $(SQRT)$(PI) = 1.128379... */ |
| enum real SQRT2 = 0x1.6a09e667f3bcc908b2fb1366ea958p+0L; /** $(SQRT)2 = 1.414213... */ |
| enum real SQRT1_2 = SQRT2/2; /** $(SQRT)$(HALF) = 0.707106... */ |
| // Note: Make sure the magic numbers in compiler backend for x87 match these. |
| |
| |
| /*********************************** |
| * Calculates the absolute value of a number |
| * |
| * Params: |
| * Num = (template parameter) type of number |
| * x = real number value |
| * z = complex number value |
| * y = imaginary number value |
| * |
| * Returns: |
| * The absolute value of the number. If floating-point or integral, |
| * the return type will be the same as the input; if complex or |
| * imaginary, the returned value will be the corresponding floating |
| * point type. |
| * |
| * For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) ) |
| * = hypot(z.re, z.im). |
| */ |
| Num abs(Num)(Num x) @safe pure nothrow |
| if (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)) && |
| !(is(Num* : const(ifloat*)) || is(Num* : const(idouble*)) |
| || is(Num* : const(ireal*)))) |
| { |
| static if (isFloatingPoint!(Num)) |
| return fabs(x); |
| else |
| return x >= 0 ? x : -x; |
| } |
| |
| /// ditto |
| auto abs(Num)(Num z) @safe pure nothrow @nogc |
| if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*)) |
| || is(Num* : const(creal*))) |
| { |
| return hypot(z.re, z.im); |
| } |
| |
| /// ditto |
| auto abs(Num)(Num y) @safe pure nothrow @nogc |
| if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*)) |
| || is(Num* : const(ireal*))) |
| { |
| return fabs(y.im); |
| } |
| |
| /// ditto |
| @safe pure nothrow @nogc unittest |
| { |
| assert(isIdentical(abs(-0.0L), 0.0L)); |
| assert(isNaN(abs(real.nan))); |
| assert(abs(-real.infinity) == real.infinity); |
| assert(abs(-3.2Li) == 3.2L); |
| assert(abs(71.6Li) == 71.6L); |
| assert(abs(-56) == 56); |
| assert(abs(2321312L) == 2321312L); |
| assert(abs(-1L+1i) == sqrt(2.0L)); |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| import std.meta : AliasSeq; |
| foreach (T; AliasSeq!(float, double, real)) |
| { |
| T f = 3; |
| assert(abs(f) == f); |
| assert(abs(-f) == f); |
| } |
| foreach (T; AliasSeq!(cfloat, cdouble, creal)) |
| { |
| T f = -12+3i; |
| assert(abs(f) == hypot(f.re, f.im)); |
| assert(abs(-f) == hypot(f.re, f.im)); |
| } |
| } |
| |
| /*********************************** |
| * Complex conjugate |
| * |
| * conj(x + iy) = x - iy |
| * |
| * Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2) |
| * is always a real number |
| */ |
| auto conj(Num)(Num z) @safe pure nothrow @nogc |
| if (is(Num* : const(cfloat*)) || is(Num* : const(cdouble*)) |
| || is(Num* : const(creal*))) |
| { |
| //FIXME |
| //Issue 14206 |
| static if (is(Num* : const(cdouble*))) |
| return cast(cdouble) conj(cast(creal) z); |
| else |
| return z.re - z.im*1fi; |
| } |
| |
| /** ditto */ |
| auto conj(Num)(Num y) @safe pure nothrow @nogc |
| if (is(Num* : const(ifloat*)) || is(Num* : const(idouble*)) |
| || is(Num* : const(ireal*))) |
| { |
| return -y; |
| } |
| |
| /// |
| @safe pure nothrow @nogc unittest |
| { |
| creal c = 7 + 3Li; |
| assert(conj(c) == 7-3Li); |
| ireal z = -3.2Li; |
| assert(conj(z) == -z); |
| } |
| //Issue 14206 |
| @safe pure nothrow @nogc unittest |
| { |
| cdouble c = 7 + 3i; |
| assert(conj(c) == 7-3i); |
| idouble z = -3.2i; |
| assert(conj(z) == -z); |
| } |
| //Issue 14206 |
| @safe pure nothrow @nogc unittest |
| { |
| cfloat c = 7f + 3fi; |
| assert(conj(c) == 7f-3fi); |
| ifloat z = -3.2fi; |
| assert(conj(z) == -z); |
| } |
| |
| /*********************************** |
| * Returns cosine of x. x is in radians. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH cos(x)) $(TH invalid?)) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes) ) |
| * ) |
| * Bugs: |
| * Results are undefined if |x| >= $(POWER 2,64). |
| */ |
| |
| real cos(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.cos(x); } |
| //FIXME |
| ///ditto |
| double cos(double x) @safe pure nothrow @nogc { return cos(cast(real) x); } |
| //FIXME |
| ///ditto |
| float cos(float x) @safe pure nothrow @nogc { return cos(cast(real) x); } |
| |
| @safe unittest |
| { |
| real function(real) pcos = &cos; |
| assert(pcos != null); |
| } |
| |
| /*********************************** |
| * Returns $(HTTP en.wikipedia.org/wiki/Sine, sine) of x. x is in $(HTTP en.wikipedia.org/wiki/Radian, radians). |
| * |
| * $(TABLE_SV |
| * $(TH3 x , sin(x) , invalid?) |
| * $(TD3 $(NAN) , $(NAN) , yes ) |
| * $(TD3 $(PLUSMN)0.0, $(PLUSMN)0.0, no ) |
| * $(TD3 $(PLUSMNINF), $(NAN) , yes ) |
| * ) |
| * |
| * Params: |
| * x = angle in radians (not degrees) |
| * Returns: |
| * sine of x |
| * See_Also: |
| * $(MYREF cos), $(MYREF tan), $(MYREF asin) |
| * Bugs: |
| * Results are undefined if |x| >= $(POWER 2,64). |
| */ |
| |
| real sin(real x) @safe pure nothrow @nogc { pragma(inline, true); return core.math.sin(x); } |
| //FIXME |
| ///ditto |
| double sin(double x) @safe pure nothrow @nogc { return sin(cast(real) x); } |
| //FIXME |
| ///ditto |
| float sin(float x) @safe pure nothrow @nogc { return sin(cast(real) x); } |
| |
| /// |
| @safe unittest |
| { |
| import std.math : sin, PI; |
| import std.stdio : writefln; |
| |
| void someFunc() |
| { |
| real x = 30.0; |
| auto result = sin(x * (PI / 180)); // convert degrees to radians |
| writefln("The sine of %s degrees is %s", x, result); |
| } |
| } |
| |
| @safe unittest |
| { |
| real function(real) psin = &sin; |
| assert(psin != null); |
| } |
| |
| /*********************************** |
| * Returns sine for complex and imaginary arguments. |
| * |
| * sin(z) = sin(z.re)*cosh(z.im) + cos(z.re)*sinh(z.im)i |
| * |
| * If both sin($(THETA)) and cos($(THETA)) are required, |
| * it is most efficient to use expi($(THETA)). |
| */ |
| creal sin(creal z) @safe pure nothrow @nogc |
| { |
| const creal cs = expi(z.re); |
| const creal csh = coshisinh(z.im); |
| return cs.im * csh.re + cs.re * csh.im * 1i; |
| } |
| |
| /** ditto */ |
| ireal sin(ireal y) @safe pure nothrow @nogc |
| { |
| return cosh(y.im)*1i; |
| } |
| |
| /// |
| @safe pure nothrow @nogc unittest |
| { |
| assert(sin(0.0+0.0i) == 0.0); |
| assert(sin(2.0+0.0i) == sin(2.0L) ); |
| } |
| |
| /*********************************** |
| * cosine, complex and imaginary |
| * |
| * cos(z) = cos(z.re)*cosh(z.im) - sin(z.re)*sinh(z.im)i |
| */ |
| creal cos(creal z) @safe pure nothrow @nogc |
| { |
| const creal cs = expi(z.re); |
| const creal csh = coshisinh(z.im); |
| return cs.re * csh.re - cs.im * csh.im * 1i; |
| } |
| |
| /** ditto */ |
| real cos(ireal y) @safe pure nothrow @nogc |
| { |
| return cosh(y.im); |
| } |
| |
| /// |
| @safe pure nothrow @nogc unittest |
| { |
| assert(cos(0.0+0.0i)==1.0); |
| assert(cos(1.3L+0.0i)==cos(1.3L)); |
| assert(cos(5.2Li)== cosh(5.2L)); |
| } |
| |
| /**************************************************************************** |
| * Returns tangent of x. x is in radians. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH tan(x)) $(TH invalid?)) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) |
| * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD yes)) |
| * ) |
| */ |
| |
| real tan(real x) @trusted pure nothrow @nogc |
| { |
| version (D_InlineAsm_X86) |
| { |
| asm pure nothrow @nogc |
| { |
| fld x[EBP] ; // load theta |
| fxam ; // test for oddball values |
| fstsw AX ; |
| sahf ; |
| jc trigerr ; // x is NAN, infinity, or empty |
| // 387's can handle subnormals |
| SC18: fptan ; |
| fstsw AX ; |
| sahf ; |
| jnp Clear1 ; // C2 = 1 (x is out of range) |
| |
| // Do argument reduction to bring x into range |
| fldpi ; |
| fxch ; |
| SC17: fprem1 ; |
| fstsw AX ; |
| sahf ; |
| jp SC17 ; |
| fstp ST(1) ; // remove pi from stack |
| jmp SC18 ; |
| |
| trigerr: |
| jnp Lret ; // if theta is NAN, return theta |
| fstp ST(0) ; // dump theta |
| } |
| return real.nan; |
| |
| Clear1: asm pure nothrow @nogc{ |
| fstp ST(0) ; // dump X, which is always 1 |
| } |
| |
| Lret: {} |
| } |
| else version (D_InlineAsm_X86_64) |
| { |
| version (Win64) |
| { |
| asm pure nothrow @nogc |
| { |
| fld real ptr [RCX] ; // load theta |
| } |
| } |
| else |
| { |
| asm pure nothrow @nogc |
| { |
| fld x[RBP] ; // load theta |
| } |
| } |
| asm pure nothrow @nogc |
| { |
| fxam ; // test for oddball values |
| fstsw AX ; |
| test AH,1 ; |
| jnz trigerr ; // x is NAN, infinity, or empty |
| // 387's can handle subnormals |
| SC18: fptan ; |
| fstsw AX ; |
| test AH,4 ; |
| jz Clear1 ; // C2 = 1 (x is out of range) |
| |
| // Do argument reduction to bring x into range |
| fldpi ; |
| fxch ; |
| SC17: fprem1 ; |
| fstsw AX ; |
| test AH,4 ; |
| jnz SC17 ; |
| fstp ST(1) ; // remove pi from stack |
| jmp SC18 ; |
| |
| trigerr: |
| test AH,4 ; |
| jz Lret ; // if theta is NAN, return theta |
| fstp ST(0) ; // dump theta |
| } |
| return real.nan; |
| |
| Clear1: asm pure nothrow @nogc{ |
| fstp ST(0) ; // dump X, which is always 1 |
| } |
| |
| Lret: {} |
| } |
| else |
| { |
| // Coefficients for tan(x) and PI/4 split into three parts. |
| static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple) |
| { |
| static immutable real[6] P = [ |
| 2.883414728874239697964612246732416606301E10L, |
| -2.307030822693734879744223131873392503321E9L, |
| 5.160188250214037865511600561074819366815E7L, |
| -4.249691853501233575668486667664718192660E5L, |
| 1.272297782199996882828849455156962260810E3L, |
| -9.889929415807650724957118893791829849557E-1L |
| ]; |
| static immutable real[7] Q = [ |
| 8.650244186622719093893836740197250197602E10L, |
| -4.152206921457208101480801635640958361612E10L, |
| 2.758476078803232151774723646710890525496E9L, |
| -5.733709132766856723608447733926138506824E7L, |
| 4.529422062441341616231663543669583527923E5L, |
| -1.317243702830553658702531997959756728291E3L, |
| 1.0 |
| ]; |
| |
| enum real P1 = |
| 7.853981633974483067550664827649598009884357452392578125E-1L; |
| enum real P2 = |
| 2.8605943630549158983813312792950660807511260829685741796657E-18L; |
| enum real P3 = |
| 2.1679525325309452561992610065108379921905808E-35L; |
| } |
| else |
| { |
| static immutable real[3] P = [ |
| -1.7956525197648487798769E7L, |
| 1.1535166483858741613983E6L, |
| -1.3093693918138377764608E4L, |
| ]; |
| static immutable real[5] Q = [ |
| -5.3869575592945462988123E7L, |
| 2.5008380182335791583922E7L, |
| -1.3208923444021096744731E6L, |
| 1.3681296347069295467845E4L, |
| 1.0000000000000000000000E0L, |
| ]; |
| |
| enum real P1 = 7.853981554508209228515625E-1L; |
| enum real P2 = 7.946627356147928367136046290398E-9L; |
| enum real P3 = 3.061616997868382943065164830688E-17L; |
| } |
| |
| // Special cases. |
| if (x == 0.0 || isNaN(x)) |
| return x; |
| if (isInfinity(x)) |
| return real.nan; |
| |
| // Make argument positive but save the sign. |
| bool sign = false; |
| if (signbit(x)) |
| { |
| sign = true; |
| x = -x; |
| } |
| |
| // Compute x mod PI/4. |
| real y = floor(x / PI_4); |
| // Strip high bits of integer part. |
| real z = ldexp(y, -4); |
| // Compute y - 16 * (y / 16). |
| z = y - ldexp(floor(z), 4); |
| |
| // Integer and fraction part modulo one octant. |
| int j = cast(int)(z); |
| |
| // Map zeros and singularities to origin. |
| if (j & 1) |
| { |
| j += 1; |
| y += 1.0; |
| } |
| |
| z = ((x - y * P1) - y * P2) - y * P3; |
| const real zz = z * z; |
| |
| if (zz > 1.0e-20L) |
| y = z + z * (zz * poly(zz, P) / poly(zz, Q)); |
| else |
| y = z; |
| |
| if (j & 2) |
| y = -1.0 / y; |
| |
| return (sign) ? -y : y; |
| } |
| } |
| |
| @safe nothrow @nogc unittest |
| { |
| static real[2][] vals = // angle,tan |
| [ |
| [ 0, 0], |
| [ .5, .5463024898], |
| [ 1, 1.557407725], |
| [ 1.5, 14.10141995], |
| [ 2, -2.185039863], |
| [ 2.5,-.7470222972], |
| [ 3, -.1425465431], |
| [ 3.5, .3745856402], |
| [ 4, 1.157821282], |
| [ 4.5, 4.637332055], |
| [ 5, -3.380515006], |
| [ 5.5,-.9955840522], |
| [ 6, -.2910061914], |
| [ 6.5, .2202772003], |
| [ 10, .6483608275], |
| |
| // special angles |
| [ PI_4, 1], |
| //[ PI_2, real.infinity], // PI_2 is not _exactly_ pi/2. |
| [ 3*PI_4, -1], |
| [ PI, 0], |
| [ 5*PI_4, 1], |
| //[ 3*PI_2, -real.infinity], |
| [ 7*PI_4, -1], |
| [ 2*PI, 0], |
| ]; |
| int i; |
| |
| for (i = 0; i < vals.length; i++) |
| { |
| real x = vals[i][0]; |
| real r = vals[i][1]; |
| real t = tan(x); |
| |
| //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r); |
| assert(approxEqual(r, t)); |
| |
| x = -x; |
| r = -r; |
| t = tan(x); |
| //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r); |
| assert(approxEqual(r, t)); |
| } |
| // overflow |
| assert(isNaN(tan(real.infinity))); |
| assert(isNaN(tan(-real.infinity))); |
| // NaN propagation |
| assert(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) )); |
| } |
| |
| @system unittest |
| { |
| assert(equalsDigit(tan(PI / 3), std.math.sqrt(3.0), useDigits)); |
| } |
| |
| /*************** |
| * Calculates the arc cosine of x, |
| * returning a value ranging from 0 to $(PI). |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH acos(x)) $(TH invalid?)) |
| * $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes)) |
| * $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes)) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes)) |
| * ) |
| */ |
| real acos(real x) @safe pure nothrow @nogc |
| { |
| return atan2(sqrt(1-x*x), x); |
| } |
| |
| /// ditto |
| double acos(double x) @safe pure nothrow @nogc { return acos(cast(real) x); } |
| |
| /// ditto |
| float acos(float x) @safe pure nothrow @nogc { return acos(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(equalsDigit(acos(0.5), std.math.PI / 3, useDigits)); |
| } |
| |
| /*************** |
| * Calculates the arc sine of x, |
| * returning a value ranging from -$(PI)/2 to $(PI)/2. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH asin(x)) $(TH invalid?)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) |
| * $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes)) |
| * $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes)) |
| * ) |
| */ |
| real asin(real x) @safe pure nothrow @nogc |
| { |
| return atan2(x, sqrt(1-x*x)); |
| } |
| |
| /// ditto |
| double asin(double x) @safe pure nothrow @nogc { return asin(cast(real) x); } |
| |
| /// ditto |
| float asin(float x) @safe pure nothrow @nogc { return asin(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(asin(0.5).approxEqual(PI / 6)); |
| } |
| |
| /*************** |
| * Calculates the arc tangent of x, |
| * returning a value ranging from -$(PI)/2 to $(PI)/2. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH atan(x)) $(TH invalid?)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes)) |
| * ) |
| */ |
| real atan(real x) @safe pure nothrow @nogc |
| { |
| version (InlineAsm_X86_Any) |
| { |
| return atan2(x, 1.0L); |
| } |
| else |
| { |
| // Coefficients for atan(x) |
| static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple) |
| { |
| static immutable real[9] P = [ |
| -6.880597774405940432145577545328795037141E2L, |
| -2.514829758941713674909996882101723647996E3L, |
| -3.696264445691821235400930243493001671932E3L, |
| -2.792272753241044941703278827346430350236E3L, |
| -1.148164399808514330375280133523543970854E3L, |
| -2.497759878476618348858065206895055957104E2L, |
| -2.548067867495502632615671450650071218995E1L, |
| -8.768423468036849091777415076702113400070E-1L, |
| -6.635810778635296712545011270011752799963E-4L |
| ]; |
| static immutable real[9] Q = [ |
| 2.064179332321782129643673263598686441900E3L, |
| 8.782996876218210302516194604424986107121E3L, |
| 1.547394317752562611786521896296215170819E4L, |
| 1.458510242529987155225086911411015961174E4L, |
| 7.928572347062145288093560392463784743935E3L, |
| 2.494680540950601626662048893678584497900E3L, |
| 4.308348370818927353321556740027020068897E2L, |
| 3.566239794444800849656497338030115886153E1L, |
| 1.0 |
| ]; |
| } |
| else |
| { |
| static immutable real[5] P = [ |
| -5.0894116899623603312185E1L, |
| -9.9988763777265819915721E1L, |
| -6.3976888655834347413154E1L, |
| -1.4683508633175792446076E1L, |
| -8.6863818178092187535440E-1L, |
| ]; |
| static immutable real[6] Q = [ |
| 1.5268235069887081006606E2L, |
| 3.9157570175111990631099E2L, |
| 3.6144079386152023162701E2L, |
| 1.4399096122250781605352E2L, |
| 2.2981886733594175366172E1L, |
| 1.0000000000000000000000E0L, |
| ]; |
| } |
| |
| // tan(PI/8) |
| enum real TAN_PI_8 = 0.414213562373095048801688724209698078569672L; |
| // tan(3 * PI/8) |
| enum real TAN3_PI_8 = 2.414213562373095048801688724209698078569672L; |
| |
| // Special cases. |
| if (x == 0.0) |
| return x; |
| if (isInfinity(x)) |
| return copysign(PI_2, x); |
| |
| // Make argument positive but save the sign. |
| bool sign = false; |
| if (signbit(x)) |
| { |
| sign = true; |
| x = -x; |
| } |
| |
| // Range reduction. |
| real y; |
| if (x > TAN3_PI_8) |
| { |
| y = PI_2; |
| x = -(1.0 / x); |
| } |
| else if (x > TAN_PI_8) |
| { |
| y = PI_4; |
| x = (x - 1.0)/(x + 1.0); |
| } |
| else |
| y = 0.0; |
| |
| // Rational form in x^^2. |
| const real z = x * x; |
| y = y + (poly(z, P) / poly(z, Q)) * z * x + x; |
| |
| return (sign) ? -y : y; |
| } |
| } |
| |
| /// ditto |
| double atan(double x) @safe pure nothrow @nogc { return atan(cast(real) x); } |
| |
| /// ditto |
| float atan(float x) @safe pure nothrow @nogc { return atan(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(equalsDigit(atan(std.math.sqrt(3.0)), PI / 3, useDigits)); |
| } |
| |
| /*************** |
| * Calculates the arc tangent of y / x, |
| * returning a value ranging from -$(PI) to $(PI). |
| * |
| * $(TABLE_SV |
| * $(TR $(TH y) $(TH x) $(TH atan(y, x))) |
| * $(TR $(TD $(NAN)) $(TD anything) $(TD $(NAN)) ) |
| * $(TR $(TD anything) $(TD $(NAN)) $(TD $(NAN)) ) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) ) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) $(TD $(PLUSMN)0.0) ) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(LT)0.0) $(TD $(PLUSMN)$(PI))) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD -0.0) $(TD $(PLUSMN)$(PI))) |
| * $(TR $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) $(TD $(PI)/2) ) |
| * $(TR $(TD $(LT)0.0) $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) ) |
| * $(TR $(TD $(GT)0.0) $(TD $(INFIN)) $(TD $(PLUSMN)0.0) ) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything) $(TD $(PLUSMN)$(PI)/2)) |
| * $(TR $(TD $(GT)0.0) $(TD -$(INFIN)) $(TD $(PLUSMN)$(PI)) ) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN)) $(TD $(PLUSMN)$(PI)/4)) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN)) $(TD $(PLUSMN)3$(PI)/4)) |
| * ) |
| */ |
| real atan2(real y, real x) @trusted pure nothrow @nogc |
| { |
| version (InlineAsm_X86_Any) |
| { |
| version (Win64) |
| { |
| asm pure nothrow @nogc { |
| naked; |
| fld real ptr [RDX]; // y |
| fld real ptr [RCX]; // x |
| fpatan; |
| ret; |
| } |
| } |
| else |
| { |
| asm pure nothrow @nogc { |
| fld y; |
| fld x; |
| fpatan; |
| } |
| } |
| } |
| else |
| { |
| // Special cases. |
| if (isNaN(x) || isNaN(y)) |
| return real.nan; |
| if (y == 0.0) |
| { |
| if (x >= 0 && !signbit(x)) |
| return copysign(0, y); |
| else |
| return copysign(PI, y); |
| } |
| if (x == 0.0) |
| return copysign(PI_2, y); |
| if (isInfinity(x)) |
| { |
| if (signbit(x)) |
| { |
| if (isInfinity(y)) |
| return copysign(3*PI_4, y); |
| else |
| return copysign(PI, y); |
| } |
| else |
| { |
| if (isInfinity(y)) |
| return copysign(PI_4, y); |
| else |
| return copysign(0.0, y); |
| } |
| } |
| if (isInfinity(y)) |
| return copysign(PI_2, y); |
| |
| // Call atan and determine the quadrant. |
| real z = atan(y / x); |
| |
| if (signbit(x)) |
| { |
| if (signbit(y)) |
| z = z - PI; |
| else |
| z = z + PI; |
| } |
| |
| if (z == 0.0) |
| return copysign(z, y); |
| |
| return z; |
| } |
| } |
| |
| /// ditto |
| double atan2(double y, double x) @safe pure nothrow @nogc |
| { |
| return atan2(cast(real) y, cast(real) x); |
| } |
| |
| /// ditto |
| float atan2(float y, float x) @safe pure nothrow @nogc |
| { |
| return atan2(cast(real) y, cast(real) x); |
| } |
| |
| @system unittest |
| { |
| assert(atan2(1.0, sqrt(3.0)).approxEqual(PI / 6)); |
| } |
| |
| /*********************************** |
| * Calculates the hyperbolic cosine of x. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH cosh(x)) $(TH invalid?)) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) ) |
| * ) |
| */ |
| real cosh(real x) @safe pure nothrow @nogc |
| { |
| // cosh = (exp(x)+exp(-x))/2. |
| // The naive implementation works correctly. |
| const real y = exp(x); |
| return (y + 1.0/y) * 0.5; |
| } |
| |
| /// ditto |
| double cosh(double x) @safe pure nothrow @nogc { return cosh(cast(real) x); } |
| |
| /// ditto |
| float cosh(float x) @safe pure nothrow @nogc { return cosh(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(equalsDigit(cosh(1.0), (E + 1.0 / E) / 2, useDigits)); |
| } |
| |
| /*********************************** |
| * Calculates the hyperbolic sine of x. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH sinh(x)) $(TH invalid?)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no)) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no)) |
| * ) |
| */ |
| real sinh(real x) @safe pure nothrow @nogc |
| { |
| // sinh(x) = (exp(x)-exp(-x))/2; |
| // Very large arguments could cause an overflow, but |
| // the maximum value of x for which exp(x) + exp(-x)) != exp(x) |
| // is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80. |
| if (fabs(x) > real.mant_dig * LN2) |
| { |
| return copysign(0.5 * exp(fabs(x)), x); |
| } |
| |
| const real y = expm1(x); |
| return 0.5 * y / (y+1) * (y+2); |
| } |
| |
| /// ditto |
| double sinh(double x) @safe pure nothrow @nogc { return sinh(cast(real) x); } |
| |
| /// ditto |
| float sinh(float x) @safe pure nothrow @nogc { return sinh(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(sinh(1.0).approxEqual((E - 1.0 / E) / 2)); |
| } |
| |
| /*********************************** |
| * Calculates the hyperbolic tangent of x. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH tanh(x)) $(TH invalid?)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) ) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no)) |
| * ) |
| */ |
| real tanh(real x) @safe pure nothrow @nogc |
| { |
| // tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x)) |
| if (fabs(x) > real.mant_dig * LN2) |
| { |
| return copysign(1, x); |
| } |
| |
| const real y = expm1(2*x); |
| return y / (y + 2); |
| } |
| |
| /// ditto |
| double tanh(double x) @safe pure nothrow @nogc { return tanh(cast(real) x); } |
| |
| /// ditto |
| float tanh(float x) @safe pure nothrow @nogc { return tanh(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(equalsDigit(tanh(1.0), sinh(1.0) / cosh(1.0), 15)); |
| } |
| |
| package: |
| |
| /* Returns cosh(x) + I * sinh(x) |
| * Only one call to exp() is performed. |
| */ |
| creal coshisinh(real x) @safe pure nothrow @nogc |
| { |
| // See comments for cosh, sinh. |
| if (fabs(x) > real.mant_dig * LN2) |
| { |
| const real y = exp(fabs(x)); |
| return y * 0.5 + 0.5i * copysign(y, x); |
| } |
| else |
| { |
| const real y = expm1(x); |
| return (y + 1.0 + 1.0/(y + 1.0)) * 0.5 + 0.5i * y / (y+1) * (y+2); |
| } |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| creal c = coshisinh(3.0L); |
| assert(c.re == cosh(3.0L)); |
| assert(c.im == sinh(3.0L)); |
| } |
| |
| public: |
| |
| /*********************************** |
| * Calculates the inverse hyperbolic cosine of x. |
| * |
| * Mathematically, acosh(x) = log(x + sqrt( x*x - 1)) |
| * |
| * $(TABLE_DOMRG |
| * $(DOMAIN 1..$(INFIN)), |
| * $(RANGE 0..$(INFIN)) |
| * ) |
| * |
| * $(TABLE_SV |
| * $(SVH x, acosh(x) ) |
| * $(SV $(NAN), $(NAN) ) |
| * $(SV $(LT)1, $(NAN) ) |
| * $(SV 1, 0 ) |
| * $(SV +$(INFIN),+$(INFIN)) |
| * ) |
| */ |
| real acosh(real x) @safe pure nothrow @nogc |
| { |
| if (x > 1/real.epsilon) |
| return LN2 + log(x); |
| else |
| return log(x + sqrt(x*x - 1)); |
| } |
| |
| /// ditto |
| double acosh(double x) @safe pure nothrow @nogc { return acosh(cast(real) x); } |
| |
| /// ditto |
| float acosh(float x) @safe pure nothrow @nogc { return acosh(cast(real) x); } |
| |
| |
| @system unittest |
| { |
| assert(isNaN(acosh(0.9))); |
| assert(isNaN(acosh(real.nan))); |
| assert(acosh(1.0)==0.0); |
| assert(acosh(real.infinity) == real.infinity); |
| assert(isNaN(acosh(0.5))); |
| assert(equalsDigit(acosh(cosh(3.0)), 3, useDigits)); |
| } |
| |
| /*********************************** |
| * Calculates the inverse hyperbolic sine of x. |
| * |
| * Mathematically, |
| * --------------- |
| * asinh(x) = log( x + sqrt( x*x + 1 )) // if x >= +0 |
| * asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0 |
| * ------------- |
| * |
| * $(TABLE_SV |
| * $(SVH x, asinh(x) ) |
| * $(SV $(NAN), $(NAN) ) |
| * $(SV $(PLUSMN)0, $(PLUSMN)0 ) |
| * $(SV $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN)) |
| * ) |
| */ |
| real asinh(real x) @safe pure nothrow @nogc |
| { |
| return (fabs(x) > 1 / real.epsilon) |
| // beyond this point, x*x + 1 == x*x |
| ? copysign(LN2 + log(fabs(x)), x) |
| // sqrt(x*x + 1) == 1 + x * x / ( 1 + sqrt(x*x + 1) ) |
| : copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x); |
| } |
| |
| /// ditto |
| double asinh(double x) @safe pure nothrow @nogc { return asinh(cast(real) x); } |
| |
| /// ditto |
| float asinh(float x) @safe pure nothrow @nogc { return asinh(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(isIdentical(asinh(0.0), 0.0)); |
| assert(isIdentical(asinh(-0.0), -0.0)); |
| assert(asinh(real.infinity) == real.infinity); |
| assert(asinh(-real.infinity) == -real.infinity); |
| assert(isNaN(asinh(real.nan))); |
| assert(equalsDigit(asinh(sinh(3.0)), 3, useDigits)); |
| } |
| |
| /*********************************** |
| * Calculates the inverse hyperbolic tangent of x, |
| * returning a value from ranging from -1 to 1. |
| * |
| * Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2 |
| * |
| * $(TABLE_DOMRG |
| * $(DOMAIN -$(INFIN)..$(INFIN)), |
| * $(RANGE -1 .. 1) |
| * ) |
| * $(BR) |
| * $(TABLE_SV |
| * $(SVH x, acosh(x) ) |
| * $(SV $(NAN), $(NAN) ) |
| * $(SV $(PLUSMN)0, $(PLUSMN)0) |
| * $(SV -$(INFIN), -0) |
| * ) |
| */ |
| real atanh(real x) @safe pure nothrow @nogc |
| { |
| // log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) ) |
| return 0.5 * log1p( 2 * x / (1 - x) ); |
| } |
| |
| /// ditto |
| double atanh(double x) @safe pure nothrow @nogc { return atanh(cast(real) x); } |
| |
| /// ditto |
| float atanh(float x) @safe pure nothrow @nogc { return atanh(cast(real) x); } |
| |
| |
| @system unittest |
| { |
| assert(isIdentical(atanh(0.0), 0.0)); |
| assert(isIdentical(atanh(-0.0),-0.0)); |
| assert(isNaN(atanh(real.nan))); |
| assert(isNaN(atanh(-real.infinity))); |
| assert(atanh(0.0) == 0); |
| assert(equalsDigit(atanh(tanh(0.5L)), 0.5, useDigits)); |
| } |
| |
| /***************************************** |
| * Returns x rounded to a long value using the current rounding mode. |
| * If the integer value of x is |
| * greater than long.max, the result is |
| * indeterminate. |
| */ |
| long rndtol(real x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.rndtol(x); } |
| //FIXME |
| ///ditto |
| long rndtol(double x) @safe pure nothrow @nogc { return rndtol(cast(real) x); } |
| //FIXME |
| ///ditto |
| long rndtol(float x) @safe pure nothrow @nogc { return rndtol(cast(real) x); } |
| |
| @safe unittest |
| { |
| long function(real) prndtol = &rndtol; |
| assert(prndtol != null); |
| } |
| |
| /***************************************** |
| * Returns x rounded to a long value using the FE_TONEAREST rounding mode. |
| * If the integer value of x is |
| * greater than long.max, the result is |
| * indeterminate. |
| */ |
| extern (C) real rndtonl(real x); |
| |
| /*************************************** |
| * Compute square root of x. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?)) |
| * $(TR $(TD -0.0) $(TD -0.0) $(TD no)) |
| * $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes)) |
| * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no)) |
| * ) |
| */ |
| float sqrt(float x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); } |
| |
| /// ditto |
| double sqrt(double x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); } |
| |
| /// ditto |
| real sqrt(real x) @nogc @safe pure nothrow { pragma(inline, true); return core.math.sqrt(x); } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| //ctfe |
| enum ZX80 = sqrt(7.0f); |
| enum ZX81 = sqrt(7.0); |
| enum ZX82 = sqrt(7.0L); |
| |
| assert(isNaN(sqrt(-1.0f))); |
| assert(isNaN(sqrt(-1.0))); |
| assert(isNaN(sqrt(-1.0L))); |
| } |
| |
| @safe unittest |
| { |
| float function(float) psqrtf = &sqrt; |
| assert(psqrtf != null); |
| double function(double) psqrtd = &sqrt; |
| assert(psqrtd != null); |
| real function(real) psqrtr = &sqrt; |
| assert(psqrtr != null); |
| } |
| |
| creal sqrt(creal z) @nogc @safe pure nothrow |
| { |
| creal c; |
| real x,y,w,r; |
| |
| if (z == 0) |
| { |
| c = 0 + 0i; |
| } |
| else |
| { |
| const real z_re = z.re; |
| const real z_im = z.im; |
| |
| x = fabs(z_re); |
| y = fabs(z_im); |
| if (x >= y) |
| { |
| r = y / x; |
| w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r))); |
| } |
| else |
| { |
| r = x / y; |
| w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r))); |
| } |
| |
| if (z_re >= 0) |
| { |
| c = w + (z_im / (w + w)) * 1.0i; |
| } |
| else |
| { |
| if (z_im < 0) |
| w = -w; |
| c = z_im / (w + w) + w * 1.0i; |
| } |
| } |
| return c; |
| } |
| |
| /** |
| * Calculates e$(SUPERSCRIPT x). |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)) ) |
| * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) |
| * $(TR $(TD -$(INFIN)) $(TD +0.0) ) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) ) |
| * ) |
| */ |
| real exp(real x) @trusted pure nothrow @nogc |
| { |
| version (D_InlineAsm_X86) |
| { |
| // e^^x = 2^^(LOG2E*x) |
| // (This is valid because the overflow & underflow limits for exp |
| // and exp2 are so similar). |
| return exp2(LOG2E*x); |
| } |
| else version (D_InlineAsm_X86_64) |
| { |
| // e^^x = 2^^(LOG2E*x) |
| // (This is valid because the overflow & underflow limits for exp |
| // and exp2 are so similar). |
| return exp2(LOG2E*x); |
| } |
| else |
| { |
| alias F = floatTraits!real; |
| static if (F.realFormat == RealFormat.ieeeDouble) |
| { |
| // Coefficients for exp(x) |
| static immutable real[3] P = [ |
| 9.99999999999999999910E-1L, |
| 3.02994407707441961300E-2L, |
| 1.26177193074810590878E-4L, |
| ]; |
| static immutable real[4] Q = [ |
| 2.00000000000000000009E0L, |
| 2.27265548208155028766E-1L, |
| 2.52448340349684104192E-3L, |
| 3.00198505138664455042E-6L, |
| ]; |
| |
| // C1 + C2 = LN2. |
| enum real C1 = 6.93145751953125E-1; |
| enum real C2 = 1.42860682030941723212E-6; |
| |
| // Overflow and Underflow limits. |
| enum real OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024) |
| enum real UF = -7.451332191019412076235E2; // ln(2^-1075) |
| } |
| else static if (F.realFormat == RealFormat.ieeeExtended || |
| F.realFormat == RealFormat.ieeeExtended53) |
| { |
| // Coefficients for exp(x) |
| static immutable real[3] P = [ |
| 9.9999999999999999991025E-1L, |
| 3.0299440770744196129956E-2L, |
| 1.2617719307481059087798E-4L, |
| ]; |
| static immutable real[4] Q = [ |
| 2.0000000000000000000897E0L, |
| 2.2726554820815502876593E-1L, |
| 2.5244834034968410419224E-3L, |
| 3.0019850513866445504159E-6L, |
| ]; |
| |
| // C1 + C2 = LN2. |
| enum real C1 = 6.9314575195312500000000E-1L; |
| enum real C2 = 1.4286068203094172321215E-6L; |
| |
| // Overflow and Underflow limits. |
| enum real OF = 1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384) |
| enum real UF = -1.13994985314888605586758E4L; // ln(2^-16446) |
| } |
| else static if (F.realFormat == RealFormat.ieeeQuadruple) |
| { |
| // Coefficients for exp(x) - 1 |
| static immutable real[5] P = [ |
| 9.999999999999999999999999999999999998502E-1L, |
| 3.508710990737834361215404761139478627390E-2L, |
| 2.708775201978218837374512615596512792224E-4L, |
| 6.141506007208645008909088812338454698548E-7L, |
| 3.279723985560247033712687707263393506266E-10L |
| ]; |
| static immutable real[6] Q = [ |
| 2.000000000000000000000000000000000000150E0, |
| 2.368408864814233538909747618894558968880E-1L, |
| 3.611828913847589925056132680618007270344E-3L, |
| 1.504792651814944826817779302637284053660E-5L, |
| 1.771372078166251484503904874657985291164E-8L, |
| 2.980756652081995192255342779918052538681E-12L |
| ]; |
| |
| // C1 + C2 = LN2. |
| enum real C1 = 6.93145751953125E-1L; |
| enum real C2 = 1.428606820309417232121458176568075500134E-6L; |
| |
| // Overflow and Underflow limits. |
| enum real OF = 1.135583025911358400418251384584930671458833e4L; |
| enum real UF = -1.143276959615573793352782661133116431383730e4L; |
| } |
| else |
| static assert(0, "Not implemented for this architecture"); |
| |
| // Special cases. Raises an overflow or underflow flag accordingly, |
| // except in the case for CTFE, where there are no hardware controls. |
| if (isNaN(x)) |
| return x; |
| if (x > OF) |
| return real.infinity; |
| if (x < UF) |
| return 0.0; |
| |
| // Express: e^^x = e^^g * 2^^n |
| // = e^^g * e^^(n * LOG2E) |
| // = e^^(g + n * LOG2E) |
| int n = cast(int) floor(LOG2E * x + 0.5); |
| x -= n * C1; |
| x -= n * C2; |
| |
| // Rational approximation for exponential of the fractional part: |
| // e^^x = 1 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) |
| const real xx = x * x; |
| const real px = x * poly(xx, P); |
| x = px / (poly(xx, Q) - px); |
| x = 1.0 + ldexp(x, 1); |
| |
| // Scale by power of 2. |
| x = ldexp(x, n); |
| |
| return x; |
| } |
| } |
| |
| /// ditto |
| double exp(double x) @safe pure nothrow @nogc { return exp(cast(real) x); } |
| |
| /// ditto |
| float exp(float x) @safe pure nothrow @nogc { return exp(cast(real) x); } |
| |
| @system unittest |
| { |
| assert(exp(3.0).feqrel(E * E * E) > 16); |
| } |
| |
| /** |
| * Calculates the value of the natural logarithm base (e) |
| * raised to the power of x, minus 1. |
| * |
| * For very small x, expm1(x) is more accurate |
| * than exp(x)-1. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH e$(SUPERSCRIPT x)-1) ) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) ) |
| * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) |
| * $(TR $(TD -$(INFIN)) $(TD -1.0) ) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) ) |
| * ) |
| */ |
| real expm1(real x) @trusted pure nothrow @nogc |
| { |
| version (D_InlineAsm_X86) |
| { |
| enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 |
| asm pure nothrow @nogc |
| { |
| /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. |
| * Author: Don Clugston. |
| * |
| * expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x. |
| * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y)) |
| * and 2ym1 = (2^^(y-rndint(y))-1). |
| * If 2rndy < 0.5*real.epsilon, result is -1. |
| * Implementation is otherwise the same as for exp2() |
| */ |
| naked; |
| fld real ptr [ESP+4] ; // x |
| mov AX, [ESP+4+8]; // AX = exponent and sign |
| sub ESP, 12+8; // Create scratch space on the stack |
| // [ESP,ESP+2] = scratchint |
| // [ESP+4..+6, +8..+10, +10] = scratchreal |
| // set scratchreal mantissa = 1.0 |
| mov dword ptr [ESP+8], 0; |
| mov dword ptr [ESP+8+4], 0x80000000; |
| and AX, 0x7FFF; // drop sign bit |
| cmp AX, 0x401D; // avoid InvalidException in fist |
| jae L_extreme; |
| fldl2e; |
| fmulp ST(1), ST; // y = x*log2(e) |
| fist dword ptr [ESP]; // scratchint = rndint(y) |
| fisub dword ptr [ESP]; // y - rndint(y) |
| // and now set scratchreal exponent |
| mov EAX, [ESP]; |
| add EAX, 0x3fff; |
| jle short L_largenegative; |
| cmp EAX,0x8000; |
| jge short L_largepositive; |
| mov [ESP+8+8],AX; |
| f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1 |
| fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y) |
| fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1 |
| fld1; |
| fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1 |
| faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1 |
| add ESP,12+8; |
| ret PARAMSIZE; |
| |
| L_extreme: // Extreme exponent. X is very large positive, very |
| // large negative, infinity, or NaN. |
| fxam; |
| fstsw AX; |
| test AX, 0x0400; // NaN_or_zero, but we already know x != 0 |
| jz L_was_nan; // if x is NaN, returns x |
| test AX, 0x0200; |
| jnz L_largenegative; |
| L_largepositive: |
| // Set scratchreal = real.max. |
| // squaring it will create infinity, and set overflow flag. |
| mov word ptr [ESP+8+8], 0x7FFE; |
| fstp ST(0); |
| fld real ptr [ESP+8]; // load scratchreal |
| fmul ST(0), ST; // square it, to create havoc! |
| L_was_nan: |
| add ESP,12+8; |
| ret PARAMSIZE; |
| L_largenegative: |
| fstp ST(0); |
| fld1; |
| fchs; // return -1. Underflow flag is not set. |
| add ESP,12+8; |
| ret PARAMSIZE; |
| } |
| } |
| else version (D_InlineAsm_X86_64) |
| { |
| asm pure nothrow @nogc |
| { |
| naked; |
| } |
| version (Win64) |
| { |
| asm pure nothrow @nogc |
| { |
| fld real ptr [RCX]; // x |
| mov AX,[RCX+8]; // AX = exponent and sign |
| } |
| } |
| else |
| { |
| asm pure nothrow @nogc |
| { |
| fld real ptr [RSP+8]; // x |
| mov AX,[RSP+8+8]; // AX = exponent and sign |
| } |
| } |
| asm pure nothrow @nogc |
| { |
| /* expm1() for x87 80-bit reals, IEEE754-2008 conformant. |
| * Author: Don Clugston. |
| * |
| * expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x. |
| * = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y)) |
| * and 2ym1 = (2^(y-rndint(y))-1). |
| * If 2rndy < 0.5*real.epsilon, result is -1. |
| * Implementation is otherwise the same as for exp2() |
| */ |
| sub RSP, 24; // Create scratch space on the stack |
| // [RSP,RSP+2] = scratchint |
| // [RSP+4..+6, +8..+10, +10] = scratchreal |
| // set scratchreal mantissa = 1.0 |
| mov dword ptr [RSP+8], 0; |
| mov dword ptr [RSP+8+4], 0x80000000; |
| and AX, 0x7FFF; // drop sign bit |
| cmp AX, 0x401D; // avoid InvalidException in fist |
| jae L_extreme; |
| fldl2e; |
| fmul ; // y = x*log2(e) |
| fist dword ptr [RSP]; // scratchint = rndint(y) |
| fisub dword ptr [RSP]; // y - rndint(y) |
| // and now set scratchreal exponent |
| mov EAX, [RSP]; |
| add EAX, 0x3fff; |
| jle short L_largenegative; |
| cmp EAX,0x8000; |
| jge short L_largepositive; |
| mov [RSP+8+8],AX; |
| f2xm1; // 2^(y-rndint(y)) -1 |
| fld real ptr [RSP+8] ; // 2^rndint(y) |
| fmul ST(1), ST; |
| fld1; |
| fsubp ST(1), ST; |
| fadd; |
| add RSP,24; |
| ret; |
| |
| L_extreme: // Extreme exponent. X is very large positive, very |
| // large negative, infinity, or NaN. |
| fxam; |
| fstsw AX; |
| test AX, 0x0400; // NaN_or_zero, but we already know x != 0 |
| jz L_was_nan; // if x is NaN, returns x |
| test AX, 0x0200; |
| jnz L_largenegative; |
| L_largepositive: |
| // Set scratchreal = real.max. |
| // squaring it will create infinity, and set overflow flag. |
| mov word ptr [RSP+8+8], 0x7FFE; |
| fstp ST(0); |
| fld real ptr [RSP+8]; // load scratchreal |
| fmul ST(0), ST; // square it, to create havoc! |
| L_was_nan: |
| add RSP,24; |
| ret; |
| |
| L_largenegative: |
| fstp ST(0); |
| fld1; |
| fchs; // return -1. Underflow flag is not set. |
| add RSP,24; |
| ret; |
| } |
| } |
| else |
| { |
| // Coefficients for exp(x) - 1 and overflow/underflow limits. |
| static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple) |
| { |
| static immutable real[8] P = [ |
| 2.943520915569954073888921213330863757240E8L, |
| -5.722847283900608941516165725053359168840E7L, |
| 8.944630806357575461578107295909719817253E6L, |
| -7.212432713558031519943281748462837065308E5L, |
| 4.578962475841642634225390068461943438441E4L, |
| -1.716772506388927649032068540558788106762E3L, |
| 4.401308817383362136048032038528753151144E1L, |
| -4.888737542888633647784737721812546636240E-1L |
| ]; |
| |
| static immutable real[9] Q = [ |
| 1.766112549341972444333352727998584753865E9L, |
| -7.848989743695296475743081255027098295771E8L, |
| 1.615869009634292424463780387327037251069E8L, |
| -2.019684072836541751428967854947019415698E7L, |
| 1.682912729190313538934190635536631941751E6L, |
| -9.615511549171441430850103489315371768998E4L, |
| 3.697714952261803935521187272204485251835E3L, |
| -8.802340681794263968892934703309274564037E1L, |
| 1.0 |
| ]; |
| |
| enum real OF = 1.1356523406294143949491931077970764891253E4L; |
| enum real UF = -1.143276959615573793352782661133116431383730e4L; |
| } |
| else |
| { |
| static immutable real[5] P = [ |
| -1.586135578666346600772998894928250240826E4L, |
| 2.642771505685952966904660652518429479531E3L, |
| -3.423199068835684263987132888286791620673E2L, |
| 1.800826371455042224581246202420972737840E1L, |
| -5.238523121205561042771939008061958820811E-1L, |
| ]; |
| static immutable real[6] Q = [ |
| -9.516813471998079611319047060563358064497E4L, |
| 3.964866271411091674556850458227710004570E4L, |
| -7.207678383830091850230366618190187434796E3L, |
| 7.206038318724600171970199625081491823079E2L, |
| -4.002027679107076077238836622982900945173E1L, |
| 1.0 |
| ]; |
| |
| enum real OF = 1.1356523406294143949492E4L; |
| enum real UF = -4.5054566736396445112120088E1L; |
| } |
| |
| |
| // C1 + C2 = LN2. |
| enum real C1 = 6.9314575195312500000000E-1L; |
| enum real C2 = 1.428606820309417232121458176568075500134E-6L; |
| |
| // Special cases. Raises an overflow flag, except in the case |
| // for CTFE, where there are no hardware controls. |
| if (x > OF) |
| return real.infinity; |
| if (x == 0.0) |
| return x; |
| if (x < UF) |
| return -1.0; |
| |
| // Express x = LN2 (n + remainder), remainder not exceeding 1/2. |
| int n = cast(int) floor(0.5 + x / LN2); |
| x -= n * C1; |
| x -= n * C2; |
| |
| // Rational approximation: |
| // exp(x) - 1 = x + 0.5 x^^2 + x^^3 P(x) / Q(x) |
| real px = x * poly(x, P); |
| real qx = poly(x, Q); |
| const real xx = x * x; |
| qx = x + (0.5 * xx + xx * px / qx); |
| |
| // We have qx = exp(remainder LN2) - 1, so: |
| // exp(x) - 1 = 2^^n (qx + 1) - 1 = 2^^n qx + 2^^n - 1. |
| px = ldexp(1.0, n); |
| x = px * qx + (px - 1.0); |
| |
| return x; |
| } |
| } |
| |
| |
| |
| /** |
| * Calculates 2$(SUPERSCRIPT x). |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH exp2(x)) ) |
| * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) ) |
| * $(TR $(TD -$(INFIN)) $(TD +0.0) ) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) ) |
| * ) |
| */ |
| pragma(inline, true) |
| real exp2(real x) @nogc @trusted pure nothrow |
| { |
| version (InlineAsm_X86_Any) |
| { |
| if (!__ctfe) |
| return exp2Asm(x); |
| else |
| return exp2Impl(x); |
| } |
| else |
| { |
| return exp2Impl(x); |
| } |
| } |
| |
| version (InlineAsm_X86_Any) |
| private real exp2Asm(real x) @nogc @trusted pure nothrow |
| { |
| version (D_InlineAsm_X86) |
| { |
| enum PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC); // always a multiple of 4 |
| |
| asm pure nothrow @nogc |
| { |
| /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. |
| * Author: Don Clugston. |
| * |
| * exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x)) |
| * The trick for high performance is to avoid the fscale(28cycles on core2), |
| * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. |
| * |
| * We can do frndint by using fist. BUT we can't use it for huge numbers, |
| * because it will set the Invalid Operation flag if overflow or NaN occurs. |
| * Fortunately, whenever this happens the result would be zero or infinity. |
| * |
| * We can perform fscale by directly poking into the exponent. BUT this doesn't |
| * work for the (very rare) cases where the result is subnormal. So we fall back |
| * to the slow method in that case. |
| */ |
| naked; |
| fld real ptr [ESP+4] ; // x |
| mov AX, [ESP+4+8]; // AX = exponent and sign |
| sub ESP, 12+8; // Create scratch space on the stack |
| // [ESP,ESP+2] = scratchint |
| // [ESP+4..+6, +8..+10, +10] = scratchreal |
| // set scratchreal mantissa = 1.0 |
| mov dword ptr [ESP+8], 0; |
| mov dword ptr [ESP+8+4], 0x80000000; |
| and AX, 0x7FFF; // drop sign bit |
| cmp AX, 0x401D; // avoid InvalidException in fist |
| jae L_extreme; |
| fist dword ptr [ESP]; // scratchint = rndint(x) |
| fisub dword ptr [ESP]; // x - rndint(x) |
| // and now set scratchreal exponent |
| mov EAX, [ESP]; |
| add EAX, 0x3fff; |
| jle short L_subnormal; |
| cmp EAX,0x8000; |
| jge short L_overflow; |
| mov [ESP+8+8],AX; |
| L_normal: |
| f2xm1; |
| fld1; |
| faddp ST(1), ST; // 2^^(x-rndint(x)) |
| fld real ptr [ESP+8] ; // 2^^rndint(x) |
| add ESP,12+8; |
| fmulp ST(1), ST; |
| ret PARAMSIZE; |
| |
| L_subnormal: |
| // Result will be subnormal. |
| // In this rare case, the simple poking method doesn't work. |
| // The speed doesn't matter, so use the slow fscale method. |
| fild dword ptr [ESP]; // scratchint |
| fld1; |
| fscale; |
| fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint |
| fstp ST(0); // drop scratchint |
| jmp L_normal; |
| |
| L_extreme: // Extreme exponent. X is very large positive, very |
| // large negative, infinity, or NaN. |
| fxam; |
| fstsw AX; |
| test AX, 0x0400; // NaN_or_zero, but we already know x != 0 |
| jz L_was_nan; // if x is NaN, returns x |
| // set scratchreal = real.min_normal |
| // squaring it will return 0, setting underflow flag |
| mov word ptr [ESP+8+8], 1; |
| test AX, 0x0200; |
| jnz L_waslargenegative; |
| L_overflow: |
| // Set scratchreal = real.max. |
| // squaring it will create infinity, and set overflow flag. |
| mov word ptr [ESP+8+8], 0x7FFE; |
| L_waslargenegative: |
| fstp ST(0); |
| fld real ptr [ESP+8]; // load scratchreal |
| fmul ST(0), ST; // square it, to create havoc! |
| L_was_nan: |
| add ESP,12+8; |
| ret PARAMSIZE; |
| } |
| } |
| else version (D_InlineAsm_X86_64) |
| { |
| asm pure nothrow @nogc |
| { |
| naked; |
| } |
| version (Win64) |
| { |
| asm pure nothrow @nogc |
| { |
| fld real ptr [RCX]; // x |
| mov AX,[RCX+8]; // AX = exponent and sign |
| } |
| } |
| else |
| { |
| asm pure nothrow @nogc |
| { |
| fld real ptr [RSP+8]; // x |
| mov AX,[RSP+8+8]; // AX = exponent and sign |
| } |
| } |
| asm pure nothrow @nogc |
| { |
| /* exp2() for x87 80-bit reals, IEEE754-2008 conformant. |
| * Author: Don Clugston. |
| * |
| * exp2(x) = 2^(rndint(x))* 2^(y-rndint(x)) |
| * The trick for high performance is to avoid the fscale(28cycles on core2), |
| * frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction. |
| * |
| * We can do frndint by using fist. BUT we can't use it for huge numbers, |
| * because it will set the Invalid Operation flag is overflow or NaN occurs. |
| * Fortunately, whenever this happens the result would be zero or infinity. |
| * |
| * We can perform fscale by directly poking into the exponent. BUT this doesn't |
| * work for the (very rare) cases where the result is subnormal. So we fall back |
| * to the slow method in that case. |
| */ |
| sub RSP, 24; // Create scratch space on the stack |
| // [RSP,RSP+2] = scratchint |
| // [RSP+4..+6, +8..+10, +10] = scratchreal |
| // set scratchreal mantissa = 1.0 |
| mov dword ptr [RSP+8], 0; |
| mov dword ptr [RSP+8+4], 0x80000000; |
| and AX, 0x7FFF; // drop sign bit |
| cmp AX, 0x401D; // avoid InvalidException in fist |
| jae L_extreme; |
| fist dword ptr [RSP]; // scratchint = rndint(x) |
| fisub dword ptr [RSP]; // x - rndint(x) |
| // and now set scratchreal exponent |
| mov EAX, [RSP]; |
| add EAX, 0x3fff; |
| jle short L_subnormal; |
| cmp EAX,0x8000; |
| jge short L_overflow; |
| mov [RSP+8+8],AX; |
| L_normal: |
| f2xm1; |
| fld1; |
| fadd; // 2^(x-rndint(x)) |
| fld real ptr [RSP+8] ; // 2^rndint(x) |
| add RSP,24; |
| fmulp ST(1), ST; |
| ret; |
| |
| L_subnormal: |
| // Result will be subnormal. |
| // In this rare case, the simple poking method doesn't work. |
| // The speed doesn't matter, so use the slow fscale method. |
| fild dword ptr [RSP]; // scratchint |
| fld1; |
| fscale; |
| fstp real ptr [RSP+8]; // scratchreal = 2^scratchint |
| fstp ST(0); // drop scratchint |
| jmp L_normal; |
| |
| L_extreme: // Extreme exponent. X is very large positive, very |
| // large negative, infinity, or NaN. |
| fxam; |
| fstsw AX; |
| test AX, 0x0400; // NaN_or_zero, but we already know x != 0 |
| jz L_was_nan; // if x is NaN, returns x |
| // set scratchreal = real.min |
| // squaring it will return 0, setting underflow flag |
| mov word ptr [RSP+8+8], 1; |
| test AX, 0x0200; |
| jnz L_waslargenegative; |
| L_overflow: |
| // Set scratchreal = real.max. |
| // squaring it will create infinity, and set overflow flag. |
| mov word ptr [RSP+8+8], 0x7FFE; |
| L_waslargenegative: |
| fstp ST(0); |
| fld real ptr [RSP+8]; // load scratchreal |
| fmul ST(0), ST; // square it, to create havoc! |
| L_was_nan: |
| add RSP,24; |
| ret; |
| } |
| } |
| else |
| static assert(0); |
| } |
| |
| private real exp2Impl(real x) @nogc @trusted pure nothrow |
| { |
| // Coefficients for exp2(x) |
| static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple) |
| { |
| static immutable real[5] P = [ |
| 9.079594442980146270952372234833529694788E12L, |
| 1.530625323728429161131811299626419117557E11L, |
| 5.677513871931844661829755443994214173883E8L, |
| 6.185032670011643762127954396427045467506E5L, |
| 1.587171580015525194694938306936721666031E2L |
| ]; |
| |
| static immutable real[6] Q = [ |
| 2.619817175234089411411070339065679229869E13L, |
| 1.490560994263653042761789432690793026977E12L, |
| 1.092141473886177435056423606755843616331E10L, |
| 2.186249607051644894762167991800811827835E7L, |
| 1.236602014442099053716561665053645270207E4L, |
| 1.0 |
| ]; |
| } |
| else |
| { |
| static immutable real[3] P = [ |
| 2.0803843631901852422887E6L, |
| 3.0286971917562792508623E4L, |
| 6.0614853552242266094567E1L, |
| ]; |
| static immutable real[4] Q = [ |
| 6.0027204078348487957118E6L, |
| 3.2772515434906797273099E5L, |
| 1.7492876999891839021063E3L, |
| 1.0000000000000000000000E0L, |
| ]; |
| } |
| |
| // Overflow and Underflow limits. |
| enum real OF = 16_384.0L; |
| enum real UF = -16_382.0L; |
| |
| // Special cases. Raises an overflow or underflow flag accordingly, |
| // except in the case for CTFE, where there are no hardware controls. |
| if (isNaN(x)) |
| return x; |
| if (x > OF) |
| return real.infinity; |
| if (x < UF) |
| return 0.0; |
| |
| // Separate into integer and fractional parts. |
| int n = cast(int) floor(x + 0.5); |
| x -= n; |
| |
| // Rational approximation: |
| // exp2(x) = 1.0 + 2x P(x^^2) / (Q(x^^2) - P(x^^2)) |
| const real xx = x * x; |
| const real px = x * poly(xx, P); |
| x = px / (poly(xx, Q) - px); |
| x = 1.0 + ldexp(x, 1); |
| |
| // Scale by power of 2. |
| x = ldexp(x, n); |
| |
| return x; |
| } |
| |
| /// |
| @safe unittest |
| { |
| assert(feqrel(exp2(0.5L), SQRT2) >= real.mant_dig -1); |
| assert(exp2(8.0L) == 256.0); |
| assert(exp2(-9.0L)== 1.0L/512.0); |
| } |
| |
| @safe unittest |
| { |
| version (CRuntime_Microsoft) {} else // aexp2/exp2f/exp2l not implemented |
| { |
| assert( core.stdc.math.exp2f(0.0f) == 1 ); |
| assert( core.stdc.math.exp2 (0.0) == 1 ); |
| assert( core.stdc.math.exp2l(0.0L) == 1 ); |
| } |
| } |
| |
| @system unittest |
| { |
| version (FloatingPointControlSupport) |
| { |
| FloatingPointControl ctrl; |
| if (FloatingPointControl.hasExceptionTraps) |
| ctrl.disableExceptions(FloatingPointControl.allExceptions); |
| ctrl.rounding = FloatingPointControl.roundToNearest; |
| } |
| |
| enum realFormat = floatTraits!real.realFormat; |
| static if (realFormat == RealFormat.ieeeQuadruple) |
| { |
| static immutable real[2][] exptestpoints = |
| [ // x exp(x) |
| [ 1.0L, E ], |
| [ 0.5L, 0x1.a61298e1e069bc972dfefab6df34p+0L ], |
| [ 3.0L, E*E*E ], |
| [ 0x1.6p+13L, 0x1.6e509d45728655cdb4840542acb5p+16250L ], // near overflow |
| [ 0x1.7p+13L, real.infinity ], // close overflow |
| [ 0x1p+80L, real.infinity ], // far overflow |
| [ real.infinity, real.infinity ], |
| [-0x1.18p+13L, 0x1.5e4bf54b4807034ea97fef0059a6p-12927L ], // near underflow |
| [-0x1.625p+13L, 0x1.a6bd68a39d11fec3a250cd97f524p-16358L ], // ditto |
| [-0x1.62dafp+13L, 0x0.cb629e9813b80ed4d639e875be6cp-16382L ], // near underflow - subnormal |
| [-0x1.6549p+13L, 0x0.0000000000000000000000000001p-16382L ], // ditto |
| [-0x1.655p+13L, 0 ], // close underflow |
| [-0x1p+30L, 0 ], // far underflow |
| ]; |
| } |
| else static if (realFormat == RealFormat.ieeeExtended || |
| realFormat == RealFormat.ieeeExtended53) |
| { |
| static immutable real[2][] exptestpoints = |
| [ // x exp(x) |
| [ 1.0L, E ], |
| [ 0.5L, 0x1.a61298e1e069bc97p+0L ], |
| [ 3.0L, E*E*E ], |
| [ 0x1.1p+13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow |
| [ 0x1.7p+13L, real.infinity ], // close overflow |
| [ 0x1p+80L, real.infinity ], // far overflow |
| [ real.infinity, real.infinity ], |
| [-0x1.18p+13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow |
| [-0x1.625p+13L, 0x1.a6bd68a39d11f35cp-16358L ], // ditto |
| [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L ], // near underflow - subnormal |
| [-0x1.643p+13L, 0x1p-16444L ], // ditto |
| [-0x1.645p+13L, 0 ], // close underflow |
| [-0x1p+30L, 0 ], // far underflow |
| ]; |
| } |
| else static if (realFormat == RealFormat.ieeeDouble) |
| { |
| static immutable real[2][] exptestpoints = |
| [ // x, exp(x) |
| [ 1.0L, E ], |
| [ 0.5L, 0x1.a61298e1e069cp+0L ], |
| [ 3.0L, E*E*E ], |
| [ 0x1.6p+9L, 0x1.93bf4ec282efbp+1015L ], // near overflow |
| [ 0x1.7p+9L, real.infinity ], // close overflow |
| [ 0x1p+80L, real.infinity ], // far overflow |
| [ real.infinity, real.infinity ], |
| [-0x1.6p+9L, 0x1.44a3824e5285fp-1016L ], // near underflow |
| [-0x1.64p+9L, 0x0.06f84920bb2d3p-1022L ], // near underflow - subnormal |
| [-0x1.743p+9L, 0x0.0000000000001p-1022L ], // ditto |
| [-0x1.8p+9L, 0 ], // close underflow |
| [-0x1p30L, 0 ], // far underflow |
| ]; |
| } |
| else |
| static assert(0, "No exp() tests for real type!"); |
| |
| const minEqualMantissaBits = real.mant_dig - 13; |
| real x; |
| version (IeeeFlagsSupport) IeeeFlags f; |
| foreach (ref pair; exptestpoints) |
| { |
| version (IeeeFlagsSupport) resetIeeeFlags(); |
| x = exp(pair[0]); |
| assert(feqrel(x, pair[1]) >= minEqualMantissaBits); |
| } |
| |
| // Ideally, exp(0) would not set the inexact flag. |
| // Unfortunately, fldl2e sets it! |
| // So it's not realistic to avoid setting it. |
| assert(exp(0.0L) == 1.0); |
| |
| // NaN propagation. Doesn't set flags, bcos was already NaN. |
| version (IeeeFlagsSupport) |
| { |
| resetIeeeFlags(); |
| x = exp(real.nan); |
| f = ieeeFlags; |
| assert(isIdentical(abs(x), real.nan)); |
| assert(f.flags == 0); |
| |
| resetIeeeFlags(); |
| x = exp(-real.nan); |
| f = ieeeFlags; |
| assert(isIdentical(abs(x), real.nan)); |
| assert(f.flags == 0); |
| } |
| else |
| { |
| x = exp(real.nan); |
| assert(isIdentical(abs(x), real.nan)); |
| |
| x = exp(-real.nan); |
| assert(isIdentical(abs(x), real.nan)); |
| } |
| |
| x = exp(NaN(0x123)); |
| assert(isIdentical(x, NaN(0x123))); |
| |
| // High resolution test (verified against GNU MPFR/Mathematica). |
| assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6_DF34p+0L); |
| } |
| |
| |
| /** |
| * Calculate cos(y) + i sin(y). |
| * |
| * On many CPUs (such as x86), this is a very efficient operation; |
| * almost twice as fast as calculating sin(y) and cos(y) separately, |
| * and is the preferred method when both are required. |
| */ |
| creal expi(real y) @trusted pure nothrow @nogc |
| { |
| version (InlineAsm_X86_Any) |
| { |
| version (Win64) |
| { |
| asm pure nothrow @nogc |
| { |
| naked; |
| fld real ptr [ECX]; |
| fsincos; |
| fxch ST(1), ST(0); |
| ret; |
| } |
| } |
| else |
| { |
| asm pure nothrow @nogc |
| { |
| fld y; |
| fsincos; |
| fxch ST(1), ST(0); |
| } |
| } |
| } |
| else |
| { |
| return cos(y) + sin(y)*1i; |
| } |
| } |
| |
| /// |
| @safe pure nothrow @nogc unittest |
| { |
| assert(expi(1.3e5L) == cos(1.3e5L) + sin(1.3e5L) * 1i); |
| assert(expi(0.0L) == 1L + 0.0Li); |
| } |
| |
| /********************************************************************* |
| * Separate floating point value into significand and exponent. |
| * |
| * Returns: |
| * Calculate and return $(I x) and $(I exp) such that |
| * value =$(I x)*2$(SUPERSCRIPT exp) and |
| * .5 $(LT)= |$(I x)| $(LT) 1.0 |
| * |
| * $(I x) has same sign as value. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH value) $(TH returns) $(TH exp)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0)) |
| * $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max)) |
| * $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min)) |
| * $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min)) |
| * ) |
| */ |
| T frexp(T)(const T value, out int exp) @trusted pure nothrow @nogc |
| if (isFloatingPoint!T) |
| { |
| Unqual!T vf = value; |
| ushort* vu = cast(ushort*)&vf; |
| static if (is(Unqual!T == float)) |
| int* vi = cast(int*)&vf; |
| else |
| long* vl = cast(long*)&vf; |
| int ex; |
| alias F = floatTraits!T; |
| |
| ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; |
| static if (F.realFormat == RealFormat.ieeeExtended || |
| F.realFormat == RealFormat.ieeeExtended53) |
| { |
| if (ex) |
| { // If exponent is non-zero |
| if (ex == F.EXPMASK) // infinity or NaN |
| { |
| if (*vl & 0x7FFF_FFFF_FFFF_FFFF) // NaN |
| { |
| *vl |= 0xC000_0000_0000_0000; // convert NaNS to NaNQ |
| exp = int.min; |
| } |
| else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity |
| exp = int.min; |
| else // positive infinity |
| exp = int.max; |
| |
| } |
| else |
| { |
| exp = ex - F.EXPBIAS; |
| vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) | 0x3FFE; |
| } |
| } |
| else if (!*vl) |
| { |
| // vf is +-0.0 |
| exp = 0; |
| } |
| else |
| { |
| // subnormal |
| |
| vf *= F.RECIP_EPSILON; |
| ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; |
| exp = ex - F.EXPBIAS - T.mant_dig + 1; |
| vu[F.EXPPOS_SHORT] = ((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FFE; |
| } |
| return vf; |
| } |
| else static if (F.realFormat == RealFormat.ieeeQuadruple) |
| { |
| if (ex) // If exponent is non-zero |
| { |
| if (ex == F.EXPMASK) |
| { |
| // infinity or NaN |
| if (vl[MANTISSA_LSB] | |
| (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN |
| { |
| // convert NaNS to NaNQ |
| vl[MANTISSA_MSB] |= 0x0000_8000_0000_0000; |
| exp = int.min; |
| } |
| else if (vu[F.EXPPOS_SHORT] & 0x8000) // negative infinity |
| exp = int.min; |
| else // positive infinity |
| exp = int.max; |
| } |
| else |
| { |
| exp = ex - F.EXPBIAS; |
| vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); |
| } |
| } |
| else if ((vl[MANTISSA_LSB] | |
| (vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) |
| { |
| // vf is +-0.0 |
| exp = 0; |
| } |
| else |
| { |
| // subnormal |
| vf *= F.RECIP_EPSILON; |
| ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; |
| exp = ex - F.EXPBIAS - T.mant_dig + 1; |
| vu[F.EXPPOS_SHORT] = F.EXPBIAS | (0x8000 & vu[F.EXPPOS_SHORT]); |
| } |
| return vf; |
| } |
| else static if (F.realFormat == RealFormat.ieeeDouble) |
| { |
| if (ex) // If exponent is non-zero |
| { |
| if (ex == F.EXPMASK) // infinity or NaN |
| { |
| if (*vl == 0x7FF0_0000_0000_0000) // positive infinity |
| { |
| exp = int.max; |
| } |
| else if (*vl == 0xFFF0_0000_0000_0000) // negative infinity |
| exp = int.min; |
| else |
| { // NaN |
| *vl |= 0x0008_0000_0000_0000; // convert NaNS to NaNQ |
| exp = int.min; |
| } |
| } |
| else |
| { |
| exp = (ex - F.EXPBIAS) >> 4; |
| vu[F.EXPPOS_SHORT] = cast(ushort)((0x800F & vu[F.EXPPOS_SHORT]) | 0x3FE0); |
| } |
| } |
| else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) |
| { |
| // vf is +-0.0 |
| exp = 0; |
| } |
| else |
| { |
| // subnormal |
| vf *= F.RECIP_EPSILON; |
| ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; |
| exp = ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1; |
| vu[F.EXPPOS_SHORT] = |
| cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3FE0); |
| } |
| return vf; |
| } |
| else static if (F.realFormat == RealFormat.ieeeSingle) |
| { |
| if (ex) // If exponent is non-zero |
| { |
| if (ex == F.EXPMASK) // infinity or NaN |
| { |
| if (*vi == 0x7F80_0000) // positive infinity |
| { |
| exp = int.max; |
| } |
| else if (*vi == 0xFF80_0000) // negative infinity |
| exp = int.min; |
| else |
| { // NaN |
| *vi |= 0x0040_0000; // convert NaNS to NaNQ |
| exp = int.min; |
| } |
| } |
| else |
| { |
| exp = (ex - F.EXPBIAS) >> 7; |
| vu[F.EXPPOS_SHORT] = cast(ushort)((0x807F & vu[F.EXPPOS_SHORT]) | 0x3F00); |
| } |
| } |
| else if (!(*vi & 0x7FFF_FFFF)) |
| { |
| // vf is +-0.0 |
| exp = 0; |
| } |
| else |
| { |
| // subnormal |
| vf *= F.RECIP_EPSILON; |
| ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; |
| exp = ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1; |
| vu[F.EXPPOS_SHORT] = |
| cast(ushort)(((-1 - F.EXPMASK) & vu[F.EXPPOS_SHORT]) | 0x3F00); |
| } |
| return vf; |
| } |
| else // static if (F.realFormat == RealFormat.ibmExtended) |
| { |
| assert(0, "frexp not implemented"); |
| } |
| } |
| |
| /// |
| @system unittest |
| { |
| int exp; |
| real mantissa = frexp(123.456L, exp); |
| |
| // check if values are equal to 19 decimal digits of precision |
| assert(equalsDigit(mantissa * pow(2.0L, cast(real) exp), 123.456L, 19)); |
| |
| assert(frexp(-real.nan, exp) && exp == int.min); |
| assert(frexp(real.nan, exp) && exp == int.min); |
| assert(frexp(-real.infinity, exp) == -real.infinity && exp == int.min); |
| assert(frexp(real.infinity, exp) == real.infinity && exp == int.max); |
| assert(frexp(-0.0, exp) == -0.0 && exp == 0); |
| assert(frexp(0.0, exp) == 0.0 && exp == 0); |
| } |
| |
| @safe unittest |
| { |
| import std.meta : AliasSeq; |
| import std.typecons : tuple, Tuple; |
| |
| foreach (T; AliasSeq!(real, double, float)) |
| { |
| Tuple!(T, T, int)[] vals = // x,frexp,exp |
| [ |
| tuple(T(0.0), T( 0.0 ), 0), |
| tuple(T(-0.0), T( -0.0), 0), |
| tuple(T(1.0), T( .5 ), 1), |
| tuple(T(-1.0), T( -.5 ), 1), |
| tuple(T(2.0), T( .5 ), 2), |
| tuple(T(float.min_normal/2.0f), T(.5), -126), |
| tuple(T.infinity, T.infinity, int.max), |
| tuple(-T.infinity, -T.infinity, int.min), |
| tuple(T.nan, T.nan, int.min), |
| tuple(-T.nan, -T.nan, int.min), |
| |
| // Phobos issue #16026: |
| tuple(3 * (T.min_normal * T.epsilon), T( .75), (T.min_exp - T.mant_dig) + 2) |
| ]; |
| |
| foreach (elem; vals) |
| { |
| T x = elem[0]; |
| T e = elem[1]; |
| int exp = elem[2]; |
| int eptr; |
| T v = frexp(x, eptr); |
| assert(isIdentical(e, v)); |
| assert(exp == eptr); |
| |
| } |
| |
| static if (floatTraits!(T).realFormat == RealFormat.ieeeExtended) |
| { |
| static T[3][] extendedvals = [ // x,frexp,exp |
| [0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal |
| [0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063], |
| [T.min_normal, .5, -16381], |
| [T.min_normal/2.0L, .5, -16382] // subnormal |
| ]; |
| foreach (elem; extendedvals) |
| { |
| T x = elem[0]; |
| T e = elem[1]; |
| int exp = cast(int) elem[2]; |
| int eptr; |
| T v = frexp(x, eptr); |
| assert(isIdentical(e, v)); |
| assert(exp == eptr); |
| |
| } |
| } |
| } |
| } |
| |
| @safe unittest |
| { |
| import std.meta : AliasSeq; |
| void foo() { |
| foreach (T; AliasSeq!(real, double, float)) |
| { |
| int exp; |
| const T a = 1; |
| immutable T b = 2; |
| auto c = frexp(a, exp); |
| auto d = frexp(b, exp); |
| } |
| } |
| } |
| |
| /****************************************** |
| * Extracts the exponent of x as a signed integral value. |
| * |
| * If x is not a special value, the result is the same as |
| * $(D cast(int) logb(x)). |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?)) |
| * $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes)) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no)) |
| * $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no)) |
| * ) |
| */ |
| int ilogb(T)(const T x) @trusted pure nothrow @nogc |
| if (isFloatingPoint!T) |
| { |
| import core.bitop : bsr; |
| alias F = floatTraits!T; |
| |
| union floatBits |
| { |
| T rv; |
| ushort[T.sizeof/2] vu; |
| uint[T.sizeof/4] vui; |
| static if (T.sizeof >= 8) |
| ulong[T.sizeof/8] vul; |
| } |
| floatBits y = void; |
| y.rv = x; |
| |
| int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK; |
| static if (F.realFormat == RealFormat.ieeeExtended || |
| F.realFormat == RealFormat.ieeeExtended53) |
| { |
| if (ex) |
| { |
| // If exponent is non-zero |
| if (ex == F.EXPMASK) // infinity or NaN |
| { |
| if (y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) // NaN |
| return FP_ILOGBNAN; |
| else // +-infinity |
| return int.max; |
| } |
| else |
| { |
| return ex - F.EXPBIAS - 1; |
| } |
| } |
| else if (!y.vul[0]) |
| { |
| // vf is +-0.0 |
| return FP_ILOGB0; |
| } |
| else |
| { |
| // subnormal |
| return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(y.vul[0]); |
| } |
| } |
| else static if (F.realFormat == RealFormat.ieeeQuadruple) |
| { |
| if (ex) // If exponent is non-zero |
| { |
| if (ex == F.EXPMASK) |
| { |
| // infinity or NaN |
| if (y.vul[MANTISSA_LSB] | ( y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) // NaN |
| return FP_ILOGBNAN; |
| else // +- infinity |
| return int.max; |
| } |
| else |
| { |
| return ex - F.EXPBIAS - 1; |
| } |
| } |
| else if ((y.vul[MANTISSA_LSB] | (y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) |
| { |
| // vf is +-0.0 |
| return FP_ILOGB0; |
| } |
| else |
| { |
| // subnormal |
| const ulong msb = y.vul[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF; |
| const ulong lsb = y.vul[MANTISSA_LSB]; |
| if (msb) |
| return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(msb) + 64; |
| else |
| return ex - F.EXPBIAS - T.mant_dig + 1 + bsr(lsb); |
| } |
| } |
| else static if (F.realFormat == RealFormat.ieeeDouble) |
| { |
| if (ex) // If exponent is non-zero |
| { |
| if (ex == F.EXPMASK) // infinity or NaN |
| { |
| if ((y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FF0_0000_0000_0000) // +- infinity |
| return int.max; |
| else // NaN |
| return FP_ILOGBNAN; |
| } |
| else |
| { |
| return ((ex - F.EXPBIAS) >> 4) - 1; |
| } |
| } |
| else if (!(y.vul[0] & 0x7FFF_FFFF_FFFF_FFFF)) |
| { |
| // vf is +-0.0 |
| return FP_ILOGB0; |
| } |
| else |
| { |
| // subnormal |
| enum MANTISSAMASK_64 = ((cast(ulong) F.MANTISSAMASK_INT) << 32) | 0xFFFF_FFFF; |
| return ((ex - F.EXPBIAS) >> 4) - T.mant_dig + 1 + bsr(y.vul[0] & MANTISSAMASK_64); |
| } |
| } |
| else static if (F.realFormat == RealFormat.ieeeSingle) |
| { |
| if (ex) // If exponent is non-zero |
| { |
| if (ex == F.EXPMASK) // infinity or NaN |
| { |
| if ((y.vui[0] & 0x7FFF_FFFF) == 0x7F80_0000) // +- infinity |
| return int.max; |
| else // NaN |
| return FP_ILOGBNAN; |
| } |
| else |
| { |
| return ((ex - F.EXPBIAS) >> 7) - 1; |
| } |
| } |
| else if (!(y.vui[0] & 0x7FFF_FFFF)) |
| { |
| // vf is +-0.0 |
| return FP_ILOGB0; |
| } |
| else |
| { |
| // subnormal |
| const uint mantissa = y.vui[0] & F.MANTISSAMASK_INT; |
| return ((ex - F.EXPBIAS) >> 7) - T.mant_dig + 1 + bsr(mantissa); |
| } |
| } |
| else // static if (F.realFormat == RealFormat.ibmExtended) |
| { |
| core.stdc.math.ilogbl(x); |
| } |
| } |
| /// ditto |
| int ilogb(T)(const T x) @safe pure not
|
