Report a bug
If you spot a problem with this page, click here to create a Github issue.
Improve this page
Quickly fork, edit online, and submit a pull request for this page. Requires a signed-in GitHub account. This works well for small changes. If you'd like to make larger changes you may want to consider using a local clone.

mir.math.sum

This module contains summation algorithms.
Authors:
Ilya Yaroshenko
Examples:
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: map;
auto ar = [1, 1e100, 1, -1e100].sliced.map!"a * 10_000";
const r = 20_000;
assert(r == ar.sum!"kbn");
assert(r == ar.sum!"kb2");
assert(r == ar.sum!"precise");
Examples:
import mir.ndslice.slice: sliced, slicedField;
import mir.ndslice.topology: map, iota, retro;
import mir.ndslice.concatenation: concatenation;
import mir.math.common;
auto ar = 1000
    .iota
    .map!(n => 1.7L.pow(n+1) - 1.7L.pow(n))
    ;
real d = 1.7L.pow(1000);
assert(sum!"precise"(concatenation(ar, [-d].sliced).slicedField) == -1);
assert(sum!"precise"(ar.retro, -d) == -1);
Examples:
Naive, Pairwise and Kahan algorithms can be used for user defined types.
static struct Quaternion(F)
    if (isFloatingPoint!F)
{
    F[4] rijk;

    /// + and - operator overloading
    Quaternion opBinary(string op)(auto ref const Quaternion rhs) const
        if (op == "+" || op == "-")
    {
        Quaternion ret = void;
        foreach (i, ref e; ret.rijk)
            mixin("e = rijk[i] "~op~" rhs.rijk[i];");
        return ret;
    }

    /// += and -= operator overloading
    Quaternion opOpAssign(string op)(auto ref const Quaternion rhs)
        if (op == "+" || op == "-")
    {
        foreach (i, ref e; rijk)
            mixin("e "~op~"= rhs.rijk[i];");
        return this;
    }

    ///constructor with single FP argument
    this(F f)
    {
        rijk[] = f;
    }

    ///assigment with single FP argument
    void opAssign(F f)
    {
        rijk[] = f;
    }
}

Quaternion!double q, p, r;
q.rijk = [0, 1, 2, 4];
p.rijk = [3, 4, 5, 9];
r.rijk = [3, 5, 7, 13];

assert(r == [p, q].sum!"naive");
assert(r == [p, q].sum!"pairwise");
assert(r == [p, q].sum!"kahan");
Examples:
All summation algorithms available for complex numbers.
cdouble[] ar = [1.0 + 2i, 2 + 3i, 3 + 4i, 4 + 5i];
cdouble r = 10 + 14i;
assert(r == ar.sum!"fast");
assert(r == ar.sum!"naive");
assert(r == ar.sum!"pairwise");
assert(r == ar.sum!"kahan");
version(LDC) // DMD Internal error: backend/cgxmm.c 628
{
    assert(r == ar.sum!"kbn");
    assert(r == ar.sum!"kb2");
}
assert(r == ar.sum!"precise");
Examples:
import mir.ndslice.topology: repeat, iota;

//simple integral summation
assert(sum([ 1, 2, 3, 4]) == 10);

//with initial value
assert(sum([ 1, 2, 3, 4], 5) == 15);

//with integral promotion
assert(sum([false, true, true, false, true]) == 3);
assert(sum(ubyte.max.repeat(100)) == 25_500);

//The result may overflow
assert(uint.max.repeat(3).sum()           ==  4_294_967_293U );
//But a seed can be used to change the summation primitive
assert(uint.max.repeat(3).sum(ulong.init) == 12_884_901_885UL);

//Floating point summation
assert(sum([1.0, 2.0, 3.0, 4.0]) == 10);

//Type overriding
static assert(is(typeof(sum!double([1F, 2F, 3F, 4F])) == double));
static assert(is(typeof(sum!double([1F, 2F, 3F, 4F], 5F)) == double));
assert(sum([1F, 2, 3, 4]) == 10);
assert(sum([1F, 2, 3, 4], 5F) == 15);

//Force pair-wise floating point summation on large integers
import std.math : approxEqual;
assert(iota([4096], uint.max / 2).sum(0.0)
           .approxEqual((uint.max / 2) * 4096.0 + 4096^^2 / 2));
Examples:
Precise summation
import mir.ndslice.topology: iota, map;
import core.stdc.tgmath: pow;
assert(iota(1000).map!(n => 1.7L.pow(real(n)+1) - 1.7L.pow(real(n)))
                 .sum!"precise" == -1 + 1.7L.pow(1000.0L));
Examples:
Precise summation with output range
import mir.ndslice.topology: iota, map;
import mir.math.common;
auto r = iota(1000).map!(n => 1.7L.pow(n+1) - 1.7L.pow(n));
Summator!(real, Summation.precise) s = 0.0;
s.put(r);
s -= 1.7L.pow(1000);
assert(s.sum() == -1);
Examples:
Precise summation with output range
import std.math: isFinite;
import mir.math.common;
float  M = 2.0f ^^ (float.max_exp-1);
double N = 2.0  ^^ (float.max_exp-1);
auto s = Summator!(float, Summation.precise)(0);
s += M;
s += M;
assert(float.infinity == s.sum()); //infinity
auto e = cast(Summator!(double, Summation.precise)) s;
assert(isFinite(e.sum()));
assert(N+N == e.sum()); //finite number
Examples:
Moving mean
import mir.ndslice.topology: linspace;
import mir.math.sum;
import std.array: array;

class MovingAverage
{
    Summator!(double, Summation.precise) summator;
    double[] circularBuffer;
    size_t frontIndex;

    double avg() @property const
    {
        return summator.sum() / circularBuffer.length;
    }

    this(double[] buffer)
    {
        assert(buffer.length);
        circularBuffer = buffer;
        summator = 0;
        summator.put(buffer);
    }

    ///operation without rounding
    void put(double x)
    {
        import mir.utility: swap;
        summator += x;
        swap(circularBuffer[frontIndex++], x);
        summator -= x;
        frontIndex %= circularBuffer.length;
    }
}

/// ma always keeps precise average of last 1000 elements
auto ma = new MovingAverage(linspace!double([1000], [0.0, 999]).array);
assert(ma.avg == (1000 * 999 / 2) / 1000.0);
/// move by 10 elements
foreach(x; linspace!double([10], [1000.0, 1009.0]))
    ma.put(x);
assert(ma.avg == (1010 * 1009 / 2 - 10 * 9 / 2) / 1000.0);
enum Summation: int;
Summation algorithms.
appropriate
Performs pairwise summation for floating point based types and fast summation for integral based types.
pairwise
precise
Precise summation algorithm. The value of the sum is rounded to the nearest representable floating-point number using the round-half-to-even rule. The result can differ from the exact value on X86, nextDown(proir) <= result && result <= nextUp(proir). The current implementation re-establish special value semantics across iterations (i.e. handling ±inf).
kahan
Kahan summation algorithm.
kbn
Kahan-Babuška-Neumaier summation algorithm. KBN gives more accurate results then Kahan.
kb2
Generalized Kahan-Babuška summation algorithm, order 2. KB2 gives more accurate results then Kahan and KBN.
naive
Naive algorithm (one by one).
fast
SIMD optimized summation algorithm.
struct Summator(T, Summation summation) if (isMutable!T);
Output range for summation.
this()(T n);
void put(N)(N n)
if (__traits(compiles, () { T a = n; a = n; a += n; } ));

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

void put(Range : Slice!(kind, packs, Iterator), SliceKind kind, size_t[] packs, Iterator)(Range r);
Adds n to the internal partial sums.
const T sum()();
Returns the value of the sum.
const C opCast(C : Summator!(P, _summation), P, Summation _summation)()
if (_summation == summation && isMutable!C && P.max_exp >= T.max_exp && P.mant_dig >= T.mant_dig);
Returns Summator with extended internal partial sums.
const C opCast(C)()
if (is(Unqual!C == T));
cast(C) operator overloading. Returns cast(C)sum(). See also: cast
void opAssign(T rhs);

void opOpAssign(string op : "+")(T rhs);

void opOpAssign(string op : "+")(ref const Summator rhs);

void opOpAssign(string op : "-")(T rhs);

void opOpAssign(string op : "-")(ref const Summator rhs);
Operator overloading.
Examples:
import mir.math.common;
import mir.ndslice.topology: iota, map;
auto r1 = iota(500).map!(a => 1.7L.pow(a+1) - 1.7L.pow(a));
auto r2 = iota([500], 500).map!(a => 1.7L.pow(a+1) - 1.7L.pow(a));
Summator!(real, Summation.precise) s1 = 0, s2 = 0.0;
foreach (e; r1) s1 += e;
foreach (e; r2) s2 -= e;
s1 -= s2;
s1 -= 1.7L.pow(1000);
assert(s1.sum() == -1);
const bool isNaN()();
Returns true if current sum is a NaN.
const bool isFinite()();
Returns true if current sum is finite (not infinite or NaN).
const bool isInfinity()();
Returns true if current sum is ±∞.
template sum(F, Summation summation = Summation.appropriate) if (isFloatingPoint!F && isMutable!F)

template sum(Summation summation = Summation.appropriate)

template sum(F, string summation) if (isFloatingPoint!F && isMutable!F)

template sum(string summation)
Sums elements of r, which must be a finite iterable.
A seed may be passed to sum. Not only will this seed be used as an initial value, but its type will be used if it is not specified.
Note that these specialized summing algorithms execute more primitive operations than vanilla summation. Therefore, if in certain cases maximum speed is required at expense of precision, one can use .
Returns:
The sum of all the elements in the range r.
See Also:
contains detailed documentation and examples about available summation algorithms.