std.numeric

This module is a port of a growing fragment of the numeric header in Alexander Stepanov's Standard Template Library, with a few additions.

var isIEEEQuadruple MakeLocalFft

tmpl CustomFloat CustomFloat CustomFloatParams CustomFloatParams FPTemporary isComplexLike secantMethod

alias lookup_t

enum CustomFloatFlags

struct CustomFloat GapWeightedSimilarityIncremental Stride

class Fft

union ToBinary

Types 7

enumCustomFloatFlags

Format flags for CustomFloat.

signed = 1Adds a sign bit to allow for signed numbers.

storeNormalized = 2Store values in normalized form by default. The actual precision of the significand is extended by 1 bit by assuming an implicit leading bit of 1 instead of 0. i.e. `1.nnnn` instead of `0.nnnn`. Tr...

allowDenorm = 4Stores the significand in https://en.wikipedia.org/wiki/IEEE754-1985#Denormalizednumbers form when the exponent is 0. Required to express the value 0.

infinity = 8Allows the storage of https://en.wikipedia.org/wiki/IEEE754-1985#Positiveandnegativeinfinity values.

nan = 16Allows the storage of https://en.wikipedia.org/wiki/NaN values.

probability = 32If set, select an exponent bias such that max_exp = 1. i.e. so that the maximum value is >= 1.0 and < 2.0. Ignored if the exponent bias is manually specified.

negativeUnsigned = 64If set, unsigned custom floats are assumed to be negative.

allowDenormZeroOnly = 128 | allowDenorm | storeNormalizedIf set, 0 is the only allowed https://en.wikipedia.org/wiki/IEEE754-1985#Denormalizednumbers number. Requires allowDenorm and storeNormalized.

ieee = signed | storeNormalized | allowDenorm | infinity | nanInclude all of the https://en.wikipedia.org/wiki/IEEEfloating_point options.

none = 0Include none of the above options.

private unionToBinary(F)

Fields

F set

CustomFloat!(CustomFloatParams!(min(F.sizeof * 8, 80))) get

Methods

typeof(get) opCall(F value)

structCustomFloat(uint precision, // fraction bits (23 for float) uint exponentWidth, // exponent bits (8 for float) Exponent width CustomFloatFlags flags, uint bias) if (isCorrectCustomFloat(precision, exponentWidth, flags))

ditto

Fields

precision + ((flags & Flags.storeNormalized) != 0) mant_digthe number of bits in mantissa

exponent_max - bias - ((flags & (Flags.infinity | Flags.nan)) != 0) + 1 max_expmaximum int value such that 2max_exp-1 is representable

cast(T_signed_exp) -(cast(long) bias) + 1 + ((flags & Flags.allowDenorm) != 0) min_expminimum int value such that 2min_exp-1 is representable as a normalized value

Methods

void roundedShift(T, U)(ref T sig, U shift)

void toNormalized(T, U)(ref T sig, ref U exp) const

void fromNormalized(T, U)(ref T sig, ref U exp)

size_t dig() @propertyReturns: number of decimal digits of precision

CustomFloat epsilon() @propertyReturns: smallest increment to the value 1

int max_10_exp(){ @propertyReturns: maximum int value such that 10max10exp is representable

int min_10_exp(){ @propertyReturns: minimum int value such that 10min10exp is representable

CustomFloat max() @propertyReturns: largest representable value that's not infinity

CustomFloat min_normal() @propertyReturns: smallest representable normalized value that's not 0

CustomFloat re() @property constReturns: real part

CustomFloat im() @propertyReturns: imaginary part

void opAssign(F: CustomFloat)(F input)Self assignment

void opAssign(F)(F input) if (__traits(compiles, cast(real) input))Assigns from any `real` compatible type.

F get(F)() if (staticIndexOf!(immutable F, immutable float, immutable double, immutable real) >= 0) @property constFetches the stored value either as a `float`, `double` or `real`.

real opUnary(string op)() if (__traits(compiles, mixin(op ~ `(get!real)`)) || op == "++" || op == "--")Convert the CustomFloat to a real and perform the relevant operator on the result

real opBinary(string op, T)(T b) if (__traits(compiles, mixin(`get!real` ~ op ~ `b.get!real`))) constditto

real opBinary(string op, T)(T b) if (  __traits(compiles,  mixin(`get!real` ~ op ~ `b`))  &&
        !__traits(compiles,  mixin(`get!real` ~ op ~ `b.get!real`)))

constditto

real opBinaryRight(string op, T)(T a) if (  __traits(compiles,  mixin(`a` ~ op ~ `get!real`))  &&
        !__traits(compiles,  mixin(`get!real` ~ op ~ `b`))  &&
        !__traits(compiles,  mixin(`get!real` ~ op ~ `b.get!real`)))

constditto

int opCmp(T)(auto ref T b) if (__traits(compiles, cast(real) b)) constditto

void opOpAssign(string op, T)(auto ref T b) if (__traits(compiles, mixin(`get!real` ~ op ~ `cast(real) b`)))ditto

Constructors

this(F input)Initialize from any `real` compatible type.

Nested Templates

uType(uint bits)

sType(uint bits)

toString()ditto

structGapWeightedSimilarityIncremental(Range, F = double) if (isRandomAccessRange!(Range) && hasLength!(Range))

Similar to gapWeightedSimilarity, just works in an incremental manner by first revealing the matches of length 1, then gapped matches of length 2, and so on. The memory requirement is s.length * t.length. The time complexity is s.length * t.length time for computing each step. Continuing on the previous example:

The implementation is based on the pseudocode in Fig. 4 of the paper

"Efﬁcient Computation of Gapped Substring Kernels on Large Alphabets"

by Rousu et al., with additional algorithmic and systems-level optimizations.

Fields

Range s

F currentValue

F * kl

size_t gram

F lambda

Methods

GapWeightedSimilarityIncremental opSlice() refReturns: `this`.

void popFront()Computes the match of the popFront length. Completes in s.length * t.length time.

F front() @propertyReturns: The gapped similarity at the current match length (initially 1, grows with each call to `popFront`).

bool empty() @propertyReturns: Whether there are more matches.

Constructors

this(Range s, Range t, F lambda)Constructs an object given two ranges `s` and `t` and a penalty `lambda`. Constructor completes in s.length * t.length time and computes all matches of length 1.

private aliaslookup_t = float

classFft

A class for performing fast Fourier transforms of power of two sizes. This class encapsulates a large amount of state that is reusable when performing multiple FFTs of sizes smaller than or equal to that specified in the constructor. This results in substantial speedups when performing multiple FFTs with a known maximum size. However, a free function API is provided for convenience if you need to perform a one-off FFT.

References:

en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm

Fields

lookup_t[][] negSinLookup

Methods

void enforceSize(R)(R range) const

void fftImpl(Ret, R)(Stride!R range, Ret buf) const

void fftImplPureReal(Ret, R)(R range, Ret buf) const

void butterfly(R)(R buf) const

size_t size() @property const

Complex!F[] fft(F = double, R)(R range) if (isFloatingPoint!F && isRandomAccessRange!R) constCompute the Fourier transform of range using the N log N Cooley-Tukey Algorithm. `range` must be a random-access range with slicing and a length equal to `size` as provided at the construction of ...

void fft(Ret,  R)(R range,  Ret buf) if (isRandomAccessRange!Ret &&  isComplexLike!(ElementType!Ret)  &&  hasSlicing!Ret)

constSame as the overload, but allows for the results to be stored in a user- provided buffer. The buffer must be of the same length as range, must be a random-access range, must have slicing, and must...

Complex!F[] inverseFft(F =  double,  R)(R range) if (isRandomAccessRange!R &&  isComplexLike!(ElementType!R)  &&  isFloatingPoint!F)

constComputes the inverse Fourier transform of a range. The range must be a random access range with slicing, have a length equal to the size provided at construction of this object, and contain elemen...

void inverseFft(Ret,  R)(R range,  Ret buf) if (isRandomAccessRange!Ret &&  isComplexLike!(ElementType!Ret)  &&  hasSlicing!Ret)

constInverse FFT that allows a user-supplied buffer to be provided. The buffer must be a random access range with slicing, and its elements must be some complex-like type.

Constructors

this(lookup_t[] memSpace)

this(size_t size)Create an `Fft` object for computing fast Fourier transforms of power of two sizes of `size` or smaller. `size` must be a power of two.

structStride(R)

Fields

Unqual!R range

size_t _nSteps

size_t _length

Methods

size_t length() const @property

typeof(this) save() @property

E opIndex(size_t index)

E front() @property

void popFront()

void popHalf()

bool empty() const @property

size_t nSteps() const @property

void doubleSteps()

size_t nSteps(size_t newVal) @property

Constructors

this(R range, size_t nStepsIn)

Functions 39

private fnbool isCorrectCustomFloat(uint precision, uint exponentWidth, CustomFloatFlags flags) @safe pure nothrow @nogc

private fnbool oppositeSigns(T1, T2)(T1 a, T2 b)Return true if a and b have opposite sign.

T findRoot(T,  DF,  DT)(scope  DF f,  const  T a,  const  T b,
     scope  DT tolerance) if (
     isFloatingPoint!T &&
     is(typeof(tolerance(T.init,  T.init)) :  bool)  &&
     is(typeof(f(T.init))  ==  R,  R)  &&  isFloatingPoint!R
    )

Find a real root of a real function f(x) via bracketing.

fnT findRoot(T, DF)(scope DF f, const T a, const T b)ditto

Tuple!(T,  T,  R,  R) findRoot(T,  R,  DF,  DT)(scope  DF f,
     const  T ax,  const  T bx,  const  R fax,  const  R fbx,
     scope  DT tolerance) if (
     isFloatingPoint!T &&
     is(typeof(tolerance(T.init,  T.init)) :  bool)  &&
     is(typeof(f(T.init))  ==  R)  &&  isFloatingPoint!R
    )

Find root of a real function f(x) by bracketing, allowing the termination condition to be specified.

Tuple!(T,  T,  R,  R) findRoot(T,  R,  DF)(scope  DF f,
     const  T ax,  const  T bx,  const  R fax,  const  R fbx)

ditto

T findRoot(T,  R)(scope  R delegate(T)  f,  const  T a,  const  T b,
     scope  bool  delegate(T lo,  T hi)  tolerance = (T a,  T b)  =>  false)

ditto

Tuple!(T,  "x",  Unqual!(ReturnType!DF),  "y",  T,  "error") findLocalMin(T,  DF)(
         scope  DF f,
         const  T ax,
         const  T bx,
         const  T relTolerance =  sqrt(T.epsilon),
         const  T absTolerance =  sqrt(T.epsilon),
        ) if (isFloatingPoint!T
    &&  __traits(compiles,  { T _ =  DF.init(T.init); }))

Find a real minimum of a real function `f(x)` via bracketing. Given a function `f` and a range `(ax .. bx)`, returns the value of `x` in the range which is closest to a minimum of `f(x)`. `f` is ne...

CommonType!(ElementType!(Range1),  ElementType!(Range2)) euclideanDistance(Range1,  Range2)(Range1 a,  Range2 b) if (isInputRange!(Range1)  &&  isInputRange!(Range2))

Computes https://en.wikipedia.org/wiki/Euclidean_distance between input ranges `a` and `b`. The two ranges must have the same length. The three-parameter version stops computation as soon as the di...

CommonType!(ElementType!(Range1),  ElementType!(Range2)) euclideanDistance(Range1,  Range2,  F)(Range1 a,  Range2 b,  F limit) if (isInputRange!(Range1)  &&  isInputRange!(Range2))

Ditto

CommonType!(ElementType!(Range1),  ElementType!(Range2)) dotProduct(Range1,  Range2)(Range1 a,  Range2 b) if (isInputRange!(Range1)  &&  isInputRange!(Range2)  &&
    !(isArray!(Range1)  &&  isArray!(Range2)))

Computes the https://en.wikipedia.org/wiki/Dot_product of input ranges `a` and b. The two ranges must have the same length. If both ranges define length, the check is done once; otherwise, it is do...

fnCommonType!(F1, F2) dotProduct(F1, F2)(in F1[] avector, in F2[] bvector)Ditto

fnF dotProduct(F, uint N)(const ref scope F[N] a, const ref scope F[N] b) if (N <= 16)ditto

CommonType!(ElementType!(Range1),  ElementType!(Range2)) cosineSimilarity(Range1,  Range2)(Range1 a,  Range2 b) if (isInputRange!(Range1)  &&  isInputRange!(Range2))

Computes the https://en.wikipedia.org/wiki/Cosine_similarity of input ranges `a` and b. The two ranges must have the same length. If both ranges define length, the check is done once; otherwise, it...

fnbool normalize(R)(R range, ElementType!(R) sum = 1) if (isForwardRange!(R))Normalizes values in `range` by multiplying each element with a number chosen such that values sum up to `sum`. If elements in range sum to zero, assigns sum / range.length to all. Normalization ma...

fnElementType!Range sumOfLog2s(Range)(Range r) if (isInputRange!Range && isFloatingPoint!(ElementType!Range))Compute the sum of binary logarithms of the input range `r`. The error of this method is much smaller than with a naive sum of log2.

fnElementType!Range entropy(Range)(Range r) if (isInputRange!Range)Computes https://en.wikipedia.org/wiki/Entropy(informationtheory, entropy) of input range `r` in bits. This function assumes (without checking) that the values in `r` are all in [0. For the entropy...

ElementType!Range entropy(Range,  F)(Range r,  F max) if (isInputRange!Range &&
    !is(CommonType!(ElementType!Range,  F)  ==  void))

Ditto

CommonType!(ElementType!Range1,  ElementType!Range2) kullbackLeiblerDivergence(Range1,  Range2)(Range1 a,  Range2 b) if (isInputRange!(Range1)  &&  isInputRange!(Range2))

Computes the https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence between input ranges `a` and `b`, which is the sum ai log(ai / bi). The base of logarithm is 2. The ranges are assume...

CommonType!(ElementType!Range1,  ElementType!Range2) jensenShannonDivergence(Range1,  Range2)(Range1 a,  Range2 b) if (isInputRange!Range1 &&  isInputRange!Range2 &&
     is(CommonType!(ElementType!Range1,  ElementType!Range2)))

Computes the https://en.wikipedia.org/wiki/Jensen%E2%80%93Shannon_divergence between `a` and b, which is the sum (ai log(2 ai / (ai + bi) + bi log(2 bi / (ai + bi))) / 2). The base of logarithm...

CommonType!(ElementType!Range1,  ElementType!Range2) jensenShannonDivergence(Range1,  Range2,  F)(Range1 a,  Range2 b,  F limit) if (isInputRange!Range1 &&  isInputRange!Range2 &&
     is(typeof(CommonType!(ElementType!Range1,  ElementType!Range2).init
    >=  F.init) :  bool))

Ditto

F gapWeightedSimilarity(alias  comp =  "a == b",  R1,  R2,  F)(R1 s,  R2 t,  F lambda) if (isRandomAccessRange!(R1)  &&  hasLength!(R1)  &&
     isRandomAccessRange!(R2)  &&  hasLength!(R2))

The so-called "all-lengths gap-weighted string kernel" computes a similarity measure between `s` and `t` based on all of their common subsequences of all lengths. Gapped subsequences are also inclu...

Select!(isFloatingPoint!(F),  F,  double) gapWeightedSimilarityNormalized(alias  comp =  "a == b",  R1,  R2,  F)(R1 s,  R2 t,  F lambda,  F sSelfSim =  F.init,  F tSelfSim =  F.init) if (isRandomAccessRange!(R1)  &&  hasLength!(R1)  &&
     isRandomAccessRange!(R2)  &&  hasLength!(R2))

The similarity per `gapWeightedSimilarity` has an issue in that it grows with the lengths of the two strings, even though the strings are not actually very similar. For example, the range ["Hello" ...

fnGapWeightedSimilarityIncremental!(R, F) gapWeightedSimilarityIncremental(R, F)(R r1, R r2, F penalty)Ditto

fntypeof(Unqual!(T).init % Unqual!(U).init) gcd(T, U)(T a, U b) if (isIntegral!T && isIntegral!U)Computes the greatest common divisor of `a` and `b` by using an efficient algorithm such as en.wikipedia.org/wiki/Euclideanalgorithm or en.wikipedia.org/wiki/BinaryGCD_algorithm algorithm.

private fntypeof(T.init % T.init) gcdImpl(T)(T a, T b) if (isIntegral!T)

auto gcd(T)(T a,  T b) if (!isIntegral!T &&
         is(typeof(T.init %  T.init))  &&
         is(typeof(T.init ==  0 ||  T.init >  0)))

ditto

private fnauto gcdImpl(T)(T a, T b) if (!isIntegral!T)

fntypeof(Unqual!(T).init % Unqual!(U).init) lcm(T, U)(T a, U b) if (isIntegral!T && isIntegral!U)Computes the least common multiple of `a` and `b`. Arguments are the same as gcd.

auto lcm(T)(T a,  T b) if (!isIntegral!T &&
         is(typeof(T.init %  T.init))  &&
         is(typeof(T.init ==  0 ||  T.init >  0)))

ditto

fnComplex!F[] fft(F = double, R)(R range)Convenience functions that create an `Fft` object, run the FFT or inverse FFT and return the result. Useful for one-off FFTs.

fnvoid fft(Ret, R)(R range, Ret buf)ditto

fnComplex!F[] inverseFft(F = double, R)(R range)ditto

fnvoid inverseFft(Ret, R)(R range, Ret buf)ditto

fnC swapRealImag(C)(C input)

fnsize_t decimalToFactorial(ulong decimal, ref ubyte[21] fac) @safe pure nothrow @nogcThis function transforms `decimal` value into a value in the factorial number system stored in `fac`.

fnvoid slowFourier2(Ret, R)(R range, Ret buf) if (isComplexLike!(ElementType!Ret))

fnvoid slowFourier4(Ret, R)(R range, Ret buf)

fnN roundDownToPowerOf2(N)(N num) if (isScalarType!N && !isFloatingPoint!N)

Variables 2

private enumvarisIEEEQuadruple = floatTraits!real.realFormat == RealFormat.ieeeQuadruple

private enumvar

MakeLocalFft = q{
    import core.stdc.stdlib;
    import core.exception : onOutOfMemoryError;

    auto lookupBuf = (cast(lookup_t*) malloc(range.length * 2 * lookup_t.sizeof))
                     [0 .. 2 * range.length];
    if (!lookupBuf.ptr)
        onOutOfMemoryError();

    scope(exit) free(cast(void*) lookupBuf.ptr);
    auto fftObj = scoped!Fft(lookupBuf);
}

Templates 7

tmplCustomFloatParams(uint bits)

tmplCustomFloatParams(uint precision, uint exponentWidth, CustomFloatFlags flags)

tmplCustomFloat(uint bits) if (bits == 8 || bits == 16 || bits == 32 || bits == 64 || bits == 80)

Allows user code to define custom floating-point formats. These formats are for storage only; all operations on them are performed by first implicitly extracting them to real first. After the operation is completed the result can be stored in a custom floating-point value via assignment.

tmplCustomFloat(uint precision, uint exponentWidth, CustomFloatFlags flags = CustomFloatFlags.ieee) if (((flags & flags.signed) + precision + exponentWidth) % 8 == 0 && precision + exponentWidth > 0)

ditto

tmplFPTemporary(F) if (isFloatingPoint!F)

Defines the fastest type to use when storing temporaries of a calculation intended to ultimately yield a result of type F (where F must be one of float, double, or real). When doing a multi-step computation, you may want to store intermediate results as FPTemporary!F.

The necessity of FPTemporary stems from the optimized floating-point operations and registers present in virtually all processors. When adding numbers in the example above, the addition may in fact be done in real precision internally. In that case, storing the intermediate result in double format is not only less precise, it is also (surprisingly) slower, because a conversion from real to double is performed every pass through the loop. This being a lose-lose situation, FPTemporary!F has been defined as the fastest type to use for calculations at precision F. There is no need to define a type for the most accurate calculations, as that is always real.

Finally, there is no guarantee that using FPTemporary!F will always be fastest, as the speed of floating-point calculations depends on very many factors.

tmplsecantMethod(alias fun)

Implements the secant method for finding a root of the function fun starting from points [xn_1, x_n] (ideally close to the root). Num may be float, double, or real.

Functions

Num secantMethod(Num)(Num xn_1, Num xn)

tmplisComplexLike(T)