mir.math.numeric

struct ProdAccumulator(T) if (isFloatingPoint!T);

Examples:

import mir.ndslice.slice: sliced;

ProdAccumulator!float x;
x.put([1, 2, 3].sliced);
assert(x.prod == 6f);
x.put(4);
assert(x.prod == 24f);

long exp;

F x;

alias mantissa = x;

pure nothrow @nogc @safe this(F value);

pure nothrow @nogc @safe this(long exp, F x);

pure nothrow @nogc @safe void put(U)(U e)
if (is(U : T));

pure nothrow @nogc @safe void put(ProdAccumulator!T value);

void put(Range)(Range r)
if (isIterable!Range);

void put(Range : Slice!(Iterator, N, kind), Iterator, size_t N, SliceKind kind)(Range r);

const pure nothrow @nogc @property @safe T prod() scope;

const pure nothrow @nogc @safe ProdAccumulator!T ldexp(long exp);

const pure nothrow @nogc @safe ProdAccumulator!T opUnary(string op : "-")();

const pure nothrow @nogc @safe ProdAccumulator!T opUnary(string op : "+")();

F prod(F, Range)(Range r)
if (isFloatingPoint!F && isIterable!Range);

Calculates the product of the elements of the input.

This function uses a separate exponential accumulation algorithm to calculate the product. A consequence of this is that the result must be a floating point type. To calculate the product of a type that is not implicitly convertible to a floating point type, use mir.algorithm.iteration.reduce or mir.algorithm.iteration.fold.

/++

Parameters:

Range r

finite iterable range

Returns:

The prduct of all the elements in r /

See Also:

mir.algorithm.iteration.reduce mir.algorithm.iteration.fold

F prod(F, Range)(Range r, ref long exp)
if (isFloatingPoint!F && isIterable!Range);

Parameters:

Range r	finite iterable range
long exp	value of exponent

Returns:

The mantissa, such that the product equals the mantissa times 2^^exp

prodType!Range prod(Range)(Range r)
if (isIterable!Range);

Parameters:

Range r

finite iterable range

Returns:

The prduct of all the elements in r

prodType!Range prod(Range)(Range r, ref long exp)
if (isIterable!Range);

Parameters:

Range r	finite iterable range
long exp	value of exponent

Returns:

The mantissa, such that the product equals the mantissa times 2^^exp

prodType!T prod(T)(scope const T[] ar...);

Parameters:

T[] ar

values

Returns:

The prduct of all the elements in ar

Examples:

Product of arbitrary inputs

assert(prod(1.0, 3, 4) == 12.0);
assert(prod!float(1, 3, 4) == 12f);

Examples:

Product of arrays and ranges

import mir.math.common: approxEqual;

enum l = 2.0 ^^ (double.max_exp - 1);
enum s = 2.0 ^^ -(double.max_exp - 1);
auto r = [l, l, l, s, s, s, 0.8 * 2.0 ^^ 10];

assert(r.prod == 0.8 * 2.0 ^^ 10);

// Can get the mantissa and exponent
long e;
assert(r.prod(e).approxEqual(0.8));
assert(e == 10);

Examples:

Product of vector

import mir.ndslice.slice: sliced;
import mir.algorithm.iteration: reduce;
import mir.math.common: approxEqual;

enum l = 2.0 ^^ (double.max_exp - 1);
enum s = 2.0 ^^ -(double.max_exp - 1);
auto c = 0.8;
auto u = c * 2.0 ^^ 10;
auto r = [l, l, l, s, s, s, u, u, u].sliced;
          
assert(r.prod == reduce!"a * b"(1.0, [u, u, u]));

long e;
assert(r.prod(e).approxEqual(reduce!"a * b"(1.0, [c, c, c])));
assert(e == 30);

Examples:

Product of matrix

import mir.ndslice.fuse: fuse;
import mir.algorithm.iteration: reduce;

enum l = 2.0 ^^ (double.max_exp - 1);
enum s = 2.0 ^^ -(double.max_exp - 1);
auto c = 0.8;
auto u = c * 2.0 ^^ 10;
auto r = [
    [l, l, l],
    [s, s, s],
    [u, u, u]
].fuse;
          
assert(r.prod == reduce!"a * b"(1.0, [u, u, u]));

long e;
assert(r.prod(e) == reduce!"a * b"(1.0, [c, c, c]));
assert(e == 30);

Examples:

Column prod of matrix

import mir.ndslice.fuse: fuse;
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.topology: alongDim, byDim, map;

auto x = [
    [2.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 5.0]
].fuse;

auto result = [4.0, 7.5, 7.5, 2.0, 5.25, 21.25];

// Use byDim or alongDim with map to compute mean of row/column.
assert(x.byDim!1.map!prod.all!approxEqual(result));
assert(x.alongDim!0.map!prod.all!approxEqual(result));

// FIXME
// Without using map, computes the prod of the whole slice
// assert(x.byDim!1.prod.all!approxEqual(result));
// assert(x.alongDim!0.prod.all!approxEqual(result));

Examples:

Can also set output type

import mir.ndslice.slice: sliced;
import mir.math.common: approxEqual;
import mir.ndslice.topology: repeat;

auto x = [1, 2, 3].sliced;
assert(x.prod!float == 6f);

Examples:

Product of variables whose underlying types are implicitly convertible to double also have type double

static struct Foo
{
    int x;
    alias x this;
}

auto x = prod(1, 2, 3);
assert(x == 6.0);
static assert(is(typeof(x) == double));

auto y = prod([Foo(1), Foo(2), Foo(3)]);
assert(y == 6.0);
static assert(is(typeof(y) == double));

Unqual!(DeepElementType!Range) sumOfLog2s(Range)(Range r)
if (isFloatingPoint!(DeepElementType!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.

Examples:

alias isNaN = x => x != x;

assert(sumOfLog2s(new double[0]) == 0);
assert(sumOfLog2s([0.0L]) == -real.infinity);
assert(sumOfLog2s([-0.0L]) == -real.infinity);
assert(sumOfLog2s([2.0L]) == 1);
assert(isNaN(sumOfLog2s([-2.0L])));
assert(isNaN(sumOfLog2s([real.nan])));
assert(isNaN(sumOfLog2s([-real.nan])));
assert(sumOfLog2s([real.infinity]) == real.infinity);
assert(isNaN(sumOfLog2s([-real.infinity])));
assert(sumOfLog2s([ 0.25, 0.25, 0.25, 0.125 ]) == -9);

Unqual!(DeepElementType!Range) sumOfLogs(Range)(Range r)
if (isFloatingPoint!(DeepElementType!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 log.

Examples:

import mir.math: LN2, approxEqual;
assert([ 0.25, 0.25, 0.25, 0.125 ].sumOfLogs.approxEqual(-9 * double(LN2)));

auto factorial(uint coefficientSize = 128)(ulong count, ulong start = 1)
if (coefficientSize % (size_t.sizeof * 8) == 0 && (coefficientSize >= size_t.sizeof * 8));

Quickly computes factorial using extended precision floating point type mir.bignum.fp.

Parameters:

ulong count	number of product members
ulong start	initial member value (optional)

Returns:

(count + start - 1)! / (start - 1)!

Complexity O(count)

Examples:

import mir.bignum.fp: Fp;
import mir.math.common: approxEqual;
import mir.math.numeric: prod;
import mir.ndslice.topology: iota;

static assert(is(typeof(factorial(33)) == Fp!128));
static assert(is(typeof(factorial!256(33)) == Fp!256));
static immutable double f33 = 8.68331761881188649551819440128e+36;
static assert(approxEqual(cast(double) factorial(33), f33));

assert(cast(double) factorial(0) == 1);
assert(cast(double) factorial(0, 100) == 1);
assert(cast(double) factorial(1, 100) == 100);
assert(approxEqual(cast(double) factorial(100, 1000), iota([100], 1000).prod!double));

auto binomialCoefficient(uint coefficientSize = 128)(ulong n, uint k)
if (coefficientSize % (size_t.sizeof * 8) == 0 && (coefficientSize >= size_t.sizeof * 8));

Quickly computes binomial coefficient using extended precision floating point type mir.bignum.fp.

Parameters:

ulong n	number elements in the set
uint k	number elements in the subset

Returns:

n choose k

Complexity O(min(k, n - k))

Examples:

import mir.bignum.fp: Fp;
import mir.math.common: approxEqual;

static assert(is(typeof(binomialCoefficient(30, 18)) == Fp!128), typeof(binomialCoefficient(30, 18)).stringof);
static assert(cast(double) binomialCoefficient(30, 18) == 86493225);

assert(approxEqual(cast(double) binomialCoefficient(100, 40), 1.374623414580281150126736972e+28));