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 std.algorithm.iteration: map;
auto ar = [1, 1e100, 1, -1e100].map!(a => a*10_000);
const r = 20_000;
assert(r == ar.sum!(Summation.kbn));
assert(r == ar.sum!(Summation.kb2));
assert(r == ar.sum!(Summation.precise));
Examples:
import std.algorithm.iteration: map;
import std.math, std.range;
auto ar = 1000
    .iota
    .map!(n => 1.7L.pow(n+1) - 1.7L.pow(n))
    .array
    ;
real d = 1.7L.pow(1000);
assert(sum!(Summation.precise)(ar.chain([-d])) == -1);
assert(sum!(Summation.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!(Summation.naive));
assert(r == [p, q].sum!(Summation.pairwise));
assert(r == [p, q].sum!(Summation.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!(Summation.fast));
assert(r == ar.sum!(Summation.naive));
assert(r == ar.sum!(Summation.pairwise));
assert(r == ar.sum!(Summation.kahan));
version(LDC) // DMD Internal error: backend/cgxmm.c 628
{
    assert(r == ar.sum!(Summation.kbn));
    assert(r == ar.sum!(Summation.kb2));
}
assert(r == ar.sum!(Summation.precise));
Examples:
import std.range: 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(uint.max / 2, uint.max / 2 + 4096).sum(0.0)
           .approxEqual((uint.max / 2) * 4096.0 + 4096^^2 / 2));
Examples:
Precise summation
import std.range: iota;
import std.algorithm.iteration: map;
import core.stdc.tgmath: pow;
assert(iota(1000).map!(n => 1.7L.pow(real(n)+1) - 1.7L.pow(real(n)))
                 .sum!(Summation.precise) == -1 + 1.7L.pow(1000.0L));
Examples:
Precise summation with output range
import std.range: iota;
import std.algorithm.iteration: map;
import std.math;
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;
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 std.range;
import mir.math.sum;

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.empty);
        circularBuffer = buffer;
        summator = 0;
        summator.put(buffer);
    }

    ///operation without rounding
    void put(double x)
    {
        import std.algorithm.mutation : 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(iota(0.0, 1000.0).array);
assert(ma.avg == (1000 * 999 / 2) / 1000.0);
/// move by 10 elements
put(ma, iota(1000.0, 1010.0));
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 + 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 (isInputRange!Range);
Adds n to the internal partial sums.
void unsafePut(F x);
Adds x to the internal partial sums. This operation doesn't re-establish special value semantics across iterations (i.e. handling -inf + inf).

Preconditions: isFinite(x).

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 std.math, std.range;
import std.algorithm.iteration: map;
auto r1 = iota(500).map!(a => 1.7L.pow(a+1) - 1.7L.pow(a));
auto r2 = iota(500, 1000).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)
Sums elements of r, which must be a finite . Although conceptually sum(r) is equivalent to reduce!((a, b) => a + b)(0, r), sum uses specialized algorithms to maximize accuracy.
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.