// 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:
* Does not work correctly for signed intergal types and value `Num`.min.
*/
auto abs(Num)(Num x) @nogc pure nothrow
if ((is(immutable Num == immutable short) || is(immutable Num == immutable byte)) ||
(is(typeof(Num.init >= 0)) && is(typeof(-Num.init))))
{
static if (isFloatingPoint!(Num))
return fabs(x);
else
{
static if (is(immutable Num == immutable short) || is(immutable Num == immutable byte))
return x >= 0 ? x : cast(Num) -int(x);
else
return x >= 0 ? x : -x;
}
}
/// ditto
@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);
}
@safe pure nothrow @nogc unittest
{
short s = -8;
byte b = -8;
assert(abs(s) == 8);
assert(abs(b) == 8);
immutable(byte) c = -8;
assert(abs(c) == 8);
}
@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 : 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(Unqual!T) (val < 0 ? -(T(1) << bsr(0 - val) + type) : T(1) << bsr(val) + type);
else
return cast(Unqual!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;
}