| // Written in the D programming language. |
| |
| /** |
| This is a submodule of $(MREF std, math). |
| |
| It contains classical algebraic functions like `abs`, `sqrt`, and `poly`. |
| |
| Copyright: Copyright The D Language Foundation 2000 - 2011. |
| 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/algebraic.d) |
| |
| Macros: |
| TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> |
| <caption>Special Values</caption> |
| $0</table> |
| NAN = $(RED NAN) |
| POWER = $1<sup>$2</sup> |
| SUB = $1<sub>$2</sub> |
| PLUSMN = ± |
| INFIN = ∞ |
| PLUSMNINF = ±∞ |
| LT = < |
| |
| */ |
| |
| module std.math.algebraic; |
| |
| static import core.math; |
| static import core.stdc.math; |
| import std.traits : CommonType, isFloatingPoint, isIntegral, isSigned, Unqual; |
| |
| /*********************************** |
| * Calculates the absolute value of a number. |
| * |
| * Params: |
| * Num = (template parameter) type of number |
| * x = real number value |
| * |
| * Returns: |
| * The absolute value of the number. If floating-point or integral, |
| * the return type will be the same as the input. |
| * |
| * Limitations: |
| * When x is a signed integral equal to `Num.min` the value of x will be returned instead. |
| * Note for 2's complement; `-Num.min` (= `Num.max + 1`) is not representable due to overflow. |
| */ |
| auto abs(Num)(Num x) @nogc nothrow pure |
| if (isIntegral!Num || (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)))) |
| { |
| static if (isFloatingPoint!(Num)) |
| return fabs(x); |
| else |
| { |
| static if (isIntegral!Num) |
| return x >= 0 ? x : cast(Num) -x; |
| else |
| return x >= 0 ? x : -x; |
| } |
| } |
| |
| /// |
| @safe pure nothrow @nogc unittest |
| { |
| import std.math.traits : isIdentical, isNaN; |
| |
| assert(isIdentical(abs(-0.0L), 0.0L)); |
| assert(isNaN(abs(real.nan))); |
| assert(abs(-real.infinity) == real.infinity); |
| assert(abs(-56) == 56); |
| assert(abs(2321312L) == 2321312L); |
| assert(abs(23u) == 23u); |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| assert(abs(byte(-8)) == 8); |
| assert(abs(ubyte(8u)) == 8); |
| assert(abs(short(-8)) == 8); |
| assert(abs(ushort(8u)) == 8); |
| assert(abs(int(-8)) == 8); |
| assert(abs(uint(8u)) == 8); |
| assert(abs(long(-8)) == 8); |
| assert(abs(ulong(8u)) == 8); |
| assert(is(typeof(abs(byte(-8))) == byte)); |
| assert(is(typeof(abs(ubyte(8u))) == ubyte)); |
| assert(is(typeof(abs(short(-8))) == short)); |
| assert(is(typeof(abs(ushort(8u))) == ushort)); |
| assert(is(typeof(abs(int(-8))) == int)); |
| assert(is(typeof(abs(uint(8u))) == uint)); |
| assert(is(typeof(abs(long(-8))) == long)); |
| assert(is(typeof(abs(ulong(8u))) == ulong)); |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| import std.meta : AliasSeq; |
| static foreach (T; AliasSeq!(float, double, real)) |
| {{ |
| T f = 3; |
| assert(abs(f) == f); |
| assert(abs(-f) == f); |
| }} |
| } |
| |
| // see https://issues.dlang.org/show_bug.cgi?id=20205 |
| // to avoid falling into the trap again |
| @safe pure nothrow @nogc unittest |
| { |
| assert(50 - abs(-100) == -50); |
| } |
| |
| // https://issues.dlang.org/show_bug.cgi?id=19162 |
| @safe unittest |
| { |
| struct Vector(T, int size) |
| { |
| T x, y, z; |
| } |
| |
| static auto abs(T, int size)(auto ref const Vector!(T, size) v) |
| { |
| return v; |
| } |
| Vector!(int, 3) v; |
| assert(abs(v) == v); |
| } |
| |
| /******************************* |
| * Returns |x| |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH fabs(x))) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) ) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) ) |
| * ) |
| */ |
| pragma(inline, true) |
| real fabs(real x) @safe pure nothrow @nogc { return core.math.fabs(x); } |
| |
| ///ditto |
| pragma(inline, true) |
| double fabs(double x) @safe pure nothrow @nogc { return core.math.fabs(x); } |
| |
| ///ditto |
| pragma(inline, true) |
| float fabs(float x) @safe pure nothrow @nogc { return core.math.fabs(x); } |
| |
| /// |
| @safe unittest |
| { |
| import std.math.traits : isIdentical; |
| |
| assert(isIdentical(fabs(0.0f), 0.0f)); |
| assert(isIdentical(fabs(-0.0f), 0.0f)); |
| assert(fabs(-10.0f) == 10.0f); |
| |
| assert(isIdentical(fabs(0.0), 0.0)); |
| assert(isIdentical(fabs(-0.0), 0.0)); |
| assert(fabs(-10.0) == 10.0); |
| |
| assert(isIdentical(fabs(0.0L), 0.0L)); |
| assert(isIdentical(fabs(-0.0L), 0.0L)); |
| assert(fabs(-10.0L) == 10.0L); |
| } |
| |
| @safe unittest |
| { |
| real function(real) pfabs = &fabs; |
| assert(pfabs != null); |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| float f = fabs(-2.0f); |
| assert(f == 2); |
| |
| double d = fabs(-2.0); |
| assert(d == 2); |
| |
| real r = fabs(-2.0L); |
| assert(r == 2); |
| } |
| |
| /*************************************** |
| * 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)) |
| * ) |
| */ |
| pragma(inline, true) |
| float sqrt(float x) @nogc @safe pure nothrow { return core.math.sqrt(x); } |
| |
| /// ditto |
| pragma(inline, true) |
| double sqrt(double x) @nogc @safe pure nothrow { return core.math.sqrt(x); } |
| |
| /// ditto |
| pragma(inline, true) |
| real sqrt(real x) @nogc @safe pure nothrow { return core.math.sqrt(x); } |
| |
| /// |
| @safe pure nothrow @nogc unittest |
| { |
| import std.math.operations : feqrel; |
| import std.math.traits : isNaN; |
| |
| assert(sqrt(2.0).feqrel(1.4142) > 16); |
| assert(sqrt(9.0).feqrel(3.0) > 16); |
| |
| assert(isNaN(sqrt(-1.0f))); |
| assert(isNaN(sqrt(-1.0))); |
| assert(isNaN(sqrt(-1.0L))); |
| } |
| |
| @safe unittest |
| { |
| // https://issues.dlang.org/show_bug.cgi?id=5305 |
| float function(float) psqrtf = &sqrt; |
| assert(psqrtf != null); |
| double function(double) psqrtd = &sqrt; |
| assert(psqrtd != null); |
| real function(real) psqrtr = &sqrt; |
| assert(psqrtr != null); |
| |
| //ctfe |
| enum ZX80 = sqrt(7.0f); |
| enum ZX81 = sqrt(7.0); |
| enum ZX82 = sqrt(7.0L); |
| } |
| |
| @safe pure nothrow @nogc unittest |
| { |
| float f = sqrt(2.0f); |
| assert(fabs(f * f - 2.0f) < .00001); |
| |
| double d = sqrt(2.0); |
| assert(fabs(d * d - 2.0) < .00001); |
| |
| real r = sqrt(2.0L); |
| assert(fabs(r * r - 2.0) < .00001); |
| } |
| |
| /*************** |
| * Calculates the cube root of x. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH $(I x)) $(TH cbrt(x)) $(TH invalid?)) |
| * $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) ) |
| * $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) ) |
| * $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) ) |
| * ) |
| */ |
| real cbrt(real x) @trusted nothrow @nogc |
| { |
| version (CRuntime_Microsoft) |
| { |
| import std.math.traits : copysign; |
| import std.math.exponential : exp2; |
| |
| version (INLINE_YL2X) |
| return copysign(exp2(core.math.yl2x(fabs(x), 1.0L/3.0L)), x); |
| else |
| return core.stdc.math.cbrtl(x); |
| } |
| else |
| return core.stdc.math.cbrtl(x); |
| } |
| |
| /// |
| @safe unittest |
| { |
| import std.math.operations : feqrel; |
| |
| assert(cbrt(1.0).feqrel(1.0) > 16); |
| assert(cbrt(27.0).feqrel(3.0) > 16); |
| assert(cbrt(15.625).feqrel(2.5) > 16); |
| } |
| |
| /*********************************************************************** |
| * Calculates the length of the |
| * hypotenuse of a right-angled triangle with sides of length x and y. |
| * The hypotenuse is the value of the square root of |
| * the sums of the squares of x and y: |
| * |
| * sqrt($(POWER x, 2) + $(POWER y, 2)) |
| * |
| * Note that hypot(x, y), hypot(y, x) and |
| * hypot(x, -y) are equivalent. |
| * |
| * $(TABLE_SV |
| * $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?)) |
| * $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no)) |
| * $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no)) |
| * $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no)) |
| * ) |
| */ |
| T hypot(T)(const T x, const T y) @safe pure nothrow @nogc |
| if (isFloatingPoint!T) |
| { |
| // Scale x and y to avoid underflow and overflow. |
| // If one is huge and the other tiny, return the larger. |
| // If both are huge, avoid overflow by scaling by 2^^-N. |
| // If both are tiny, avoid underflow by scaling by 2^^N. |
| import core.math : fabs, sqrt; |
| import std.math.traits : floatTraits, RealFormat; |
| |
| alias F = floatTraits!T; |
| |
| T u = fabs(x); |
| T v = fabs(y); |
| if (!(u >= v)) // check for NaN as well. |
| { |
| v = u; |
| u = fabs(y); |
| if (u == T.infinity) return u; // hypot(inf, nan) == inf |
| if (v == T.infinity) return v; // hypot(nan, inf) == inf |
| } |
| |
| static if (F.realFormat == RealFormat.ieeeSingle) |
| { |
| enum SQRTMIN = 0x1p-60f; |
| enum SQRTMAX = 0x1p+60f; |
| enum SCALE_UNDERFLOW = 0x1p+90f; |
| enum SCALE_OVERFLOW = 0x1p-90f; |
| } |
| else static if (F.realFormat == RealFormat.ieeeDouble || |
| F.realFormat == RealFormat.ieeeExtended53 || |
| F.realFormat == RealFormat.ibmExtended) |
| { |
| enum SQRTMIN = 0x1p-450L; |
| enum SQRTMAX = 0x1p+500L; |
| enum SCALE_UNDERFLOW = 0x1p+600L; |
| enum SCALE_OVERFLOW = 0x1p-600L; |
| } |
| else static if (F.realFormat == RealFormat.ieeeExtended || |
| F.realFormat == RealFormat.ieeeQuadruple) |
| { |
| enum SQRTMIN = 0x1p-8000L; |
| enum SQRTMAX = 0x1p+8000L; |
| enum SCALE_UNDERFLOW = 0x1p+10000L; |
| enum SCALE_OVERFLOW = 0x1p-10000L; |
| } |
| else |
| assert(0, "hypot not implemented"); |
| |
| // Now u >= v, or else one is NaN. |
| T ratio = 1.0; |
| if (v >= SQRTMAX) |
| { |
| // hypot(huge, huge) -- avoid overflow |
| ratio = SCALE_UNDERFLOW; |
| u *= SCALE_OVERFLOW; |
| v *= SCALE_OVERFLOW; |
| } |
| else if (u <= SQRTMIN) |
| { |
| // hypot (tiny, tiny) -- avoid underflow |
| // This is only necessary to avoid setting the underflow |
| // flag. |
| ratio = SCALE_OVERFLOW; |
| u *= SCALE_UNDERFLOW; |
| v *= SCALE_UNDERFLOW; |
| } |
| |
| if (u * T.epsilon > v) |
| { |
| // hypot (huge, tiny) = huge |
| return u; |
| } |
| |
| // both are in the normal range |
| return ratio * sqrt(u*u + v*v); |
| } |
| |
| /// |
| @safe unittest |
| { |
| import std.math.operations : feqrel; |
| |
| assert(hypot(1.0, 1.0).feqrel(1.4142) > 16); |
| assert(hypot(3.0, 4.0).feqrel(5.0) > 16); |
| assert(hypot(real.infinity, 1.0L) == real.infinity); |
| assert(hypot(real.infinity, real.nan) == real.infinity); |
| } |
| |
| @safe unittest |
| { |
| import std.math.operations : feqrel; |
| |
| assert(hypot(1.0f, 1.0f).feqrel(1.4142f) > 16); |
| assert(hypot(3.0f, 4.0f).feqrel(5.0f) > 16); |
| assert(hypot(float.infinity, 1.0f) == float.infinity); |
| assert(hypot(float.infinity, float.nan) == float.infinity); |
| |
| assert(hypot(1.0L, 1.0L).feqrel(1.4142L) > 16); |
| assert(hypot(3.0L, 4.0L).feqrel(5.0L) > 16); |
| assert(hypot(double.infinity, 1.0) == double.infinity); |
| assert(hypot(double.infinity, double.nan) == double.infinity); |
| } |
| |
| @safe unittest |
| { |
| import std.math.operations : feqrel; |
| import std.math.traits : isIdentical; |
| import std.meta : AliasSeq; |
| |
| static foreach (T; AliasSeq!(float, double, real)) |
| {{ |
| static T[3][] vals = // x,y,hypot |
| [ |
| [ 0.0, 0.0, 0.0], |
| [ 0.0, -0.0, 0.0], |
| [ -0.0, -0.0, 0.0], |
| [ 3.0, 4.0, 5.0], |
| [ -300, -400, 500], |
| [0.0, 7.0, 7.0], |
| [9.0, 9*T.epsilon, 9.0], |
| [88/(64*sqrt(T.min_normal)), 105/(64*sqrt(T.min_normal)), 137/(64*sqrt(T.min_normal))], |
| [88/(128*sqrt(T.min_normal)), 105/(128*sqrt(T.min_normal)), 137/(128*sqrt(T.min_normal))], |
| [3*T.min_normal*T.epsilon, 4*T.min_normal*T.epsilon, 5*T.min_normal*T.epsilon], |
| [ T.min_normal, T.min_normal, sqrt(2.0L)*T.min_normal], |
| [ T.max/sqrt(2.0L), T.max/sqrt(2.0L), T.max], |
| [ T.infinity, T.nan, T.infinity], |
| [ T.nan, T.infinity, T.infinity], |
| [ T.nan, T.nan, T.nan], |
| [ T.nan, T.max, T.nan], |
| [ T.max, T.nan, T.nan], |
| ]; |
| for (int i = 0; i < vals.length; i++) |
| { |
| T x = vals[i][0]; |
| T y = vals[i][1]; |
| T z = vals[i][2]; |
| T h = hypot(x, y); |
| assert(isIdentical(z,h) || feqrel(z, h) >= T.mant_dig - 1); |
| } |
| }} |
| } |
| |
| /*********************************************************************** |
| * Calculates the distance of the point (x, y, z) from the origin (0, 0, 0) |
| * in three-dimensional space. |
| * The distance is the value of the square root of the sums of the squares |
| * of x, y, and z: |
| * |
| * sqrt($(POWER x, 2) + $(POWER y, 2) + $(POWER z, 2)) |
| * |
| * Note that the distance between two points (x1, y1, z1) and (x2, y2, z2) |
| * in three-dimensional space can be calculated as hypot(x2-x1, y2-y1, z2-z1). |
| * |
| * Params: |
| * x = floating point value |
| * y = floating point value |
| * z = floating point value |
| * |
| * Returns: |
| * The square root of the sum of the squares of the given arguments. |
| */ |
| T hypot(T)(const T x, const T y, const T z) @safe pure nothrow @nogc |
| if (isFloatingPoint!T) |
| { |
| import core.math : fabs, sqrt; |
| import std.math.operations : fmax; |
| const absx = fabs(x); |
| const absy = fabs(y); |
| const absz = fabs(z); |
| |
| // Scale all parameters to avoid overflow. |
| const ratio = fmax(absx, fmax(absy, absz)); |
| if (ratio == 0.0) |
| return ratio; |
| |
| return ratio * sqrt((absx / ratio) * (absx / ratio) |
| + (absy / ratio) * (absy / ratio) |
| + (absz / ratio) * (absz / ratio)); |
| } |
| |
| /// |
| @safe unittest |
| { |
| import std.math.operations : isClose; |
| |
| assert(isClose(hypot(1.0, 2.0, 2.0), 3.0)); |
| assert(isClose(hypot(2.0, 3.0, 6.0), 7.0)); |
| assert(isClose(hypot(1.0, 4.0, 8.0), 9.0)); |
| } |
| |
| @safe unittest |
| { |
| import std.meta : AliasSeq; |
| import std.math.traits : isIdentical; |
| import std.math.operations : isClose; |
| static foreach (T; AliasSeq!(float, double, real)) |
| {{ |
| static T[4][] vals = [ |
| [ 0.0L, 0.0L, 0.0L, 0.0L ], |
| [ 0.0L, 1.0L, 1.0L, sqrt(2.0L) ], |
| [ 1.0L, 1.0L, 1.0L, sqrt(3.0L) ], |
| [ 1.0L, 2.0L, 2.0L, 3.0L ], |
| [ 2.0L, 3.0L, 6.0L, 7.0L ], |
| [ 1.0L, 4.0L, 8.0L, 9.0L ], |
| [ 4.0L, 4.0L, 7.0L, 9.0L ], |
| [ 12.0L, 16.0L, 21.0L, 29.0L ], |
| [ 1e+8L, 1.0L, 1e-8L, 1e+8L+5e-9L ], |
| [ 1.0L, 1e+8L, 1e-8L, 1e+8L+5e-9L ], |
| [ 1e-8L, 1.0L, 1e+8L, 1e+8L+5e-9L ], |
| [ 1e-2L, 1e-4L, 1e-4L, 0.010000999950004999375L ], |
| [ 2147483647.0L, 2147483647.0L, 2147483647.0L, 3719550785.027307813987L ] |
| ]; |
| for (int i = 0; i < vals.length; i++) |
| { |
| T x = vals[i][0]; |
| T y = vals[i][1]; |
| T z = vals[i][2]; |
| T r = vals[i][3]; |
| T a = hypot(x, y, z); |
| assert(isIdentical(r, a) || isClose(r, a)); |
| } |
| }} |
| } |
| |
| /*********************************** |
| * Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2) + |
| * $(SUB a,3)$(POWER x,3); ... |
| * |
| * Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2) + |
| * x($(SUB a, 3) + ...))) |
| * Params: |
| * x = the value to evaluate. |
| * A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc. |
| */ |
| Unqual!(CommonType!(T1, T2)) poly(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc |
| if (isFloatingPoint!T1 && isFloatingPoint!T2) |
| in |
| { |
| assert(A.length > 0); |
| } |
| do |
| { |
| static if (is(immutable T2 == immutable real)) |
| { |
| return polyImpl(x, A); |
| } |
| else |
| { |
| return polyImplBase(x, A); |
| } |
| } |
| |
| /// ditto |
| Unqual!(CommonType!(T1, T2)) poly(T1, T2, int N)(T1 x, ref const T2[N] A) @safe pure nothrow @nogc |
| if (isFloatingPoint!T1 && isFloatingPoint!T2 && N > 0 && N <= 10) |
| { |
| // statically unrolled version for up to 10 coefficients |
| typeof(return) r = A[N - 1]; |
| static foreach (i; 1 .. N) |
| { |
| r *= x; |
| r += A[N - 1 - i]; |
| } |
| return r; |
| } |
| |
| /// |
| @safe nothrow @nogc unittest |
| { |
| real x = 3.1L; |
| static real[] pp = [56.1L, 32.7L, 6]; |
| |
| assert(poly(x, pp) == (56.1L + (32.7L + 6.0L * x) * x)); |
| } |
| |
| @safe nothrow @nogc unittest |
| { |
| double x = 3.1; |
| static double[] pp = [56.1, 32.7, 6]; |
| double y = x; |
| y *= 6.0; |
| y += 32.7; |
| y *= x; |
| y += 56.1; |
| assert(poly(x, pp) == y); |
| } |
| |
| @safe unittest |
| { |
| static assert(poly(3.0, [1.0, 2.0, 3.0]) == 34); |
| } |
| |
| private Unqual!(CommonType!(T1, T2)) polyImplBase(T1, T2)(T1 x, in T2[] A) @trusted pure nothrow @nogc |
| if (isFloatingPoint!T1 && isFloatingPoint!T2) |
| { |
| ptrdiff_t i = A.length - 1; |
| typeof(return) r = A[i]; |
| while (--i >= 0) |
| { |
| r *= x; |
| r += A[i]; |
| } |
| return r; |
| } |
| |
| version (linux) version = GenericPosixVersion; |
| else version (FreeBSD) version = GenericPosixVersion; |
| else version (OpenBSD) version = GenericPosixVersion; |
| else version (Solaris) version = GenericPosixVersion; |
| else version (DragonFlyBSD) version = GenericPosixVersion; |
| |
| private real polyImpl(real x, in real[] A) @trusted pure nothrow @nogc |
| { |
| version (D_InlineAsm_X86) |
| { |
| if (__ctfe) |
| { |
| return polyImplBase(x, A); |
| } |
| version (Windows) |
| { |
| // BUG: This code assumes a frame pointer in EBP. |
| asm pure nothrow @nogc // assembler by W. Bright |
| { |
| // EDX = (A.length - 1) * real.sizeof |
| mov ECX,A[EBP] ; // ECX = A.length |
| dec ECX ; |
| lea EDX,[ECX][ECX*8] ; |
| add EDX,ECX ; |
| add EDX,A+4[EBP] ; |
| fld real ptr [EDX] ; // ST0 = coeff[ECX] |
| jecxz return_ST ; |
| fld x[EBP] ; // ST0 = x |
| fxch ST(1) ; // ST1 = x, ST0 = r |
| align 4 ; |
| L2: fmul ST,ST(1) ; // r *= x |
| fld real ptr -10[EDX] ; |
| sub EDX,10 ; // deg-- |
| faddp ST(1),ST ; |
| dec ECX ; |
| jne L2 ; |
| fxch ST(1) ; // ST1 = r, ST0 = x |
| fstp ST(0) ; // dump x |
| align 4 ; |
| return_ST: ; |
| } |
| } |
| else version (GenericPosixVersion) |
| { |
| asm pure nothrow @nogc // assembler by W. Bright |
| { |
| // EDX = (A.length - 1) * real.sizeof |
| mov ECX,A[EBP] ; // ECX = A.length |
| dec ECX ; |
| lea EDX,[ECX*8] ; |
| lea EDX,[EDX][ECX*4] ; |
| add EDX,A+4[EBP] ; |
| fld real ptr [EDX] ; // ST0 = coeff[ECX] |
| jecxz return_ST ; |
| fld x[EBP] ; // ST0 = x |
| fxch ST(1) ; // ST1 = x, ST0 = r |
| align 4 ; |
| L2: fmul ST,ST(1) ; // r *= x |
| fld real ptr -12[EDX] ; |
| sub EDX,12 ; // deg-- |
| faddp ST(1),ST ; |
| dec ECX ; |
| jne L2 ; |
| fxch ST(1) ; // ST1 = r, ST0 = x |
| fstp ST(0) ; // dump x |
| align 4 ; |
| return_ST: ; |
| } |
| } |
| else version (OSX) |
| { |
| asm pure nothrow @nogc // assembler by W. Bright |
| { |
| // EDX = (A.length - 1) * real.sizeof |
| mov ECX,A[EBP] ; // ECX = A.length |
| dec ECX ; |
| lea EDX,[ECX*8] ; |
| add EDX,EDX ; |
| add EDX,A+4[EBP] ; |
| fld real ptr [EDX] ; // ST0 = coeff[ECX] |
| jecxz return_ST ; |
| fld x[EBP] ; // ST0 = x |
| fxch ST(1) ; // ST1 = x, ST0 = r |
| align 4 ; |
| L2: fmul ST,ST(1) ; // r *= x |
| fld real ptr -16[EDX] ; |
| sub EDX,16 ; // deg-- |
| faddp ST(1),ST ; |
| dec ECX ; |
| jne L2 ; |
| fxch ST(1) ; // ST1 = r, ST0 = x |
| fstp ST(0) ; // dump x |
| align 4 ; |
| return_ST: ; |
| } |
| } |
| else |
| { |
| static assert(0); |
| } |
| } |
| else |
| { |
| return polyImplBase(x, A); |
| } |
| } |
| |
| /** |
| * Gives the next power of two after `val`. `T` can be any built-in |
| * numerical type. |
| * |
| * If the operation would lead to an over/underflow, this function will |
| * return `0`. |
| * |
| * Params: |
| * val = any number |
| * |
| * Returns: |
| * the next power of two after `val` |
| */ |
| T nextPow2(T)(const T val) |
| if (isIntegral!T) |
| { |
| return powIntegralImpl!(PowType.ceil)(val); |
| } |
| |
| /// ditto |
| T nextPow2(T)(const T val) |
| if (isFloatingPoint!T) |
| { |
| return powFloatingPointImpl!(PowType.ceil)(val); |
| } |
| |
| /// |
| @safe @nogc pure nothrow unittest |
| { |
| assert(nextPow2(2) == 4); |
| assert(nextPow2(10) == 16); |
| assert(nextPow2(4000) == 4096); |
| |
| assert(nextPow2(-2) == -4); |
| assert(nextPow2(-10) == -16); |
| |
| assert(nextPow2(uint.max) == 0); |
| assert(nextPow2(uint.min) == 0); |
| assert(nextPow2(size_t.max) == 0); |
| assert(nextPow2(size_t.min) == 0); |
| |
| assert(nextPow2(int.max) == 0); |
| assert(nextPow2(int.min) == 0); |
| assert(nextPow2(long.max) == 0); |
| assert(nextPow2(long.min) == 0); |
| } |
| |
| /// |
| @safe @nogc pure nothrow unittest |
| { |
| assert(nextPow2(2.1) == 4.0); |
| assert(nextPow2(-2.0) == -4.0); |
| assert(nextPow2(0.25) == 0.5); |
| assert(nextPow2(-4.0) == -8.0); |
| |
| assert(nextPow2(double.max) == 0.0); |
| assert(nextPow2(double.infinity) == double.infinity); |
| } |
| |
| @safe @nogc pure nothrow unittest |
| { |
| assert(nextPow2(ubyte(2)) == 4); |
| assert(nextPow2(ubyte(10)) == 16); |
| |
| assert(nextPow2(byte(2)) == 4); |
| assert(nextPow2(byte(10)) == 16); |
| |
| assert(nextPow2(short(2)) == 4); |
| assert(nextPow2(short(10)) == 16); |
| assert(nextPow2(short(4000)) == 4096); |
| |
| assert(nextPow2(ushort(2)) == 4); |
| assert(nextPow2(ushort(10)) == 16); |
| assert(nextPow2(ushort(4000)) == 4096); |
| } |
| |
| @safe @nogc pure nothrow unittest |
| { |
| foreach (ulong i; 1 .. 62) |
| { |
| assert(nextPow2(1UL << i) == 2UL << i); |
| assert(nextPow2((1UL << i) - 1) == 1UL << i); |
| assert(nextPow2((1UL << i) + 1) == 2UL << i); |
| assert(nextPow2((1UL << i) + (1UL<<(i-1))) == 2UL << i); |
| } |
| } |
| |
| @safe @nogc pure nothrow unittest |
| { |
| import std.math.traits : isNaN; |
| import std.meta : AliasSeq; |
| |
| static foreach (T; AliasSeq!(float, double, real)) |
| {{ |
| enum T subNormal = T.min_normal / 2; |
| |
| static if (subNormal) assert(nextPow2(subNormal) == T.min_normal); |
| |
| assert(nextPow2(T(0.0)) == 0.0); |
| |
| assert(nextPow2(T(2.0)) == 4.0); |
| assert(nextPow2(T(2.1)) == 4.0); |
| assert(nextPow2(T(3.1)) == 4.0); |
| assert(nextPow2(T(4.0)) == 8.0); |
| assert(nextPow2(T(0.25)) == 0.5); |
| |
| assert(nextPow2(T(-2.0)) == -4.0); |
| assert(nextPow2(T(-2.1)) == -4.0); |
| assert(nextPow2(T(-3.1)) == -4.0); |
| assert(nextPow2(T(-4.0)) == -8.0); |
| assert(nextPow2(T(-0.25)) == -0.5); |
| |
| assert(nextPow2(T.max) == 0); |
| assert(nextPow2(-T.max) == 0); |
| |
| assert(nextPow2(T.infinity) == T.infinity); |
| assert(nextPow2(T.init).isNaN); |
| }} |
| } |
| |
| // https://issues.dlang.org/show_bug.cgi?id=15973 |
| @safe @nogc pure nothrow unittest |
| { |
| assert(nextPow2(uint.max / 2) == uint.max / 2 + 1); |
| assert(nextPow2(uint.max / 2 + 2) == 0); |
| assert(nextPow2(int.max / 2) == int.max / 2 + 1); |
| assert(nextPow2(int.max / 2 + 2) == 0); |
| assert(nextPow2(int.min + 1) == int.min); |
| } |
| |
| /** |
| * Gives the last power of two before `val`. $(T) can be any built-in |
| * numerical type. |
| * |
| * Params: |
| * val = any number |
| * |
| * Returns: |
| * the last power of two before `val` |
| */ |
| T truncPow2(T)(const T val) |
| if (isIntegral!T) |
| { |
| return powIntegralImpl!(PowType.floor)(val); |
| } |
| |
| /// ditto |
| T truncPow2(T)(const T val) |
| if (isFloatingPoint!T) |
| { |
| return powFloatingPointImpl!(PowType.floor)(val); |
| } |
| |
| /// |
| @safe @nogc pure nothrow unittest |
| { |
| assert(truncPow2(3) == 2); |
| assert(truncPow2(4) == 4); |
| assert(truncPow2(10) == 8); |
| assert(truncPow2(4000) == 2048); |
| |
| assert(truncPow2(-5) == -4); |
| assert(truncPow2(-20) == -16); |
| |
| assert(truncPow2(uint.max) == int.max + 1); |
| assert(truncPow2(uint.min) == 0); |
| assert(truncPow2(ulong.max) == long.max + 1); |
| assert(truncPow2(ulong.min) == 0); |
| |
| assert(truncPow2(int.max) == (int.max / 2) + 1); |
| assert(truncPow2(int.min) == int.min); |
| assert(truncPow2(long.max) == (long.max / 2) + 1); |
| assert(truncPow2(long.min) == long.min); |
| } |
| |
| /// |
| @safe @nogc pure nothrow unittest |
| { |
| assert(truncPow2(2.1) == 2.0); |
| assert(truncPow2(7.0) == 4.0); |
| assert(truncPow2(-1.9) == -1.0); |
| assert(truncPow2(0.24) == 0.125); |
| assert(truncPow2(-7.0) == -4.0); |
| |
| assert(truncPow2(double.infinity) == double.infinity); |
| } |
| |
| @safe @nogc pure nothrow unittest |
| { |
| assert(truncPow2(ubyte(3)) == 2); |
| assert(truncPow2(ubyte(4)) == 4); |
| assert(truncPow2(ubyte(10)) == 8); |
| |
| assert(truncPow2(byte(3)) == 2); |
| assert(truncPow2(byte(4)) == 4); |
| assert(truncPow2(byte(10)) == 8); |
| |
| assert(truncPow2(ushort(3)) == 2); |
| assert(truncPow2(ushort(4)) == 4); |
| assert(truncPow2(ushort(10)) == 8); |
| assert(truncPow2(ushort(4000)) == 2048); |
| |
| assert(truncPow2(short(3)) == 2); |
| assert(truncPow2(short(4)) == 4); |
| assert(truncPow2(short(10)) == 8); |
| assert(truncPow2(short(4000)) == 2048); |
| } |
| |
| @safe @nogc pure nothrow unittest |
| { |
| foreach (ulong i; 1 .. 62) |
| { |
| assert(truncPow2(2UL << i) == 2UL << i); |
| assert(truncPow2((2UL << i) + 1) == 2UL << i); |
| assert(truncPow2((2UL << i) - 1) == 1UL << i); |
| assert(truncPow2((2UL << i) - (2UL<<(i-1))) == 1UL << i); |
| } |
| } |
| |
| @safe @nogc pure nothrow unittest |
| { |
| import std.math.traits : isNaN; |
| import std.meta : AliasSeq; |
| |
| static foreach (T; AliasSeq!(float, double, real)) |
| { |
| assert(truncPow2(T(0.0)) == 0.0); |
| |
| assert(truncPow2(T(4.0)) == 4.0); |
| assert(truncPow2(T(2.1)) == 2.0); |
| assert(truncPow2(T(3.5)) == 2.0); |
| assert(truncPow2(T(7.0)) == 4.0); |
| assert(truncPow2(T(0.24)) == 0.125); |
| |
| assert(truncPow2(T(-2.0)) == -2.0); |
| assert(truncPow2(T(-2.1)) == -2.0); |
| assert(truncPow2(T(-3.1)) == -2.0); |
| assert(truncPow2(T(-7.0)) == -4.0); |
| assert(truncPow2(T(-0.24)) == -0.125); |
| |
| assert(truncPow2(T.infinity) == T.infinity); |
| assert(truncPow2(T.init).isNaN); |
| } |
| } |
| |
| private enum PowType |
| { |
| floor, |
| ceil |
| } |
| |
| pragma(inline, true) |
| private T powIntegralImpl(PowType type, T)(T val) |
| { |
| import core.bitop : bsr; |
| |
| if (val == 0 || (type == PowType.ceil && (val > T.max / 2 || val == T.min))) |
| return 0; |
| else |
| { |
| static if (isSigned!T) |
| return cast() cast(T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type); |
| else |
| return cast() cast(T) (T(1) << bsr(val) + type); |
| } |
| } |
| |
| private T powFloatingPointImpl(PowType type, T)(T x) |
| { |
| import std.math.traits : copysign, isFinite; |
| import std.math.exponential : frexp; |
| |
| if (!x.isFinite) |
| return x; |
| |
| if (!x) |
| return x; |
| |
| int exp; |
| auto y = frexp(x, exp); |
| |
| static if (type == PowType.ceil) |
| y = core.math.ldexp(cast(T) 0.5, exp + 1); |
| else |
| y = core.math.ldexp(cast(T) 0.5, exp); |
| |
| if (!y.isFinite) |
| return cast(T) 0.0; |
| |
| y = copysign(y, x); |
| |
| return y; |
| } |
