If you spot a problem with this page, click here to create a GitHub issue.

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.

License:

Apache-2.0

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");
assert(r == ar.sum!"decimal");

Examples:

Decimal precise summation

auto ar = [777.7, -777];
assert(ar.sum!"decimal" == 0.7);
assert(sum!"decimal"(777.7, -777) == 0.7);

// The exact binary reuslt is 0.7000000000000455
assert(ar[0] + ar[1] == 0.7000000000000455);
assert(ar.sum!"fast" == 0.7000000000000455);
assert(ar.sum!"kahan" == 0.7000000000000455);
assert(ar.sum!"kbn" == 0.7000000000000455);
assert(ar.sum!"kb2" == 0.7000000000000455);
assert(ar.sum!"precise" == 0.7000000000000455);

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.

import mir.internal.utility: isFloatingPoint;
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 ;
        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.

import mir.complex: Complex;

auto ar = [Complex!double(1.0, 2), Complex!double(2.0, 3), Complex!double(3.0, 4), Complex!double(4.0, 5)];
Complex!double r = Complex!double(10.0, 14);
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");
assert(r == ar.sum!"decimal");

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 mir.math : approxEqual;
assert(iota!long([4096], uint.max / 2).sum(0.0)
           .approxEqual((uint.max / 2) * 4096.0 + 4096.0 * 4096.0 / 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 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(e.sum < double.infinity);
assert(N+N == e.sum()); //finite number

Examples:

Moving mean

import mir.internal.utility: isFloatingPoint;
import mir.math.sum;
import mir.ndslice.topology: linspace;
import mir.rc.array: rcarray;

struct MovingAverage(T)
    if (isFloatingPoint!T)
{
    import mir.math.stat: MeanAccumulator;

    MeanAccumulator!(T, Summation.precise) meanAccumulator;
    double[] circularBuffer;
    size_t frontIndex;

    @disable this(this);

    auto avg() @property const
    {
        return meanAccumulator.mean;
    }

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

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

/// ma always keeps precise average of last 1000 elements
auto x = linspace!double([1000], [0.0, 999]).rcarray;
auto ma = MovingAverage!double(x[]);
assert(ma.avg == (1000 * 999 / 2) / 1000.0);
/// move by 10 elements
foreach(e; linspace!double([10], [1000.0, 1009.0]))
    ma.put(e);
assert(ma.avg == (1010 * 1009 / 2 - 10 * 9 / 2) / 1000.0);

Examples:

Arbitrary sum

import mir.complex;
alias C = Complex!double;
assert(sum(1, 2, 3, 4) == 10);
assert(sum!float(1, 2, 3, 4) == 10f);
assert(sum(1f, 2, 3, 4) == 10f);
assert(sum(C(1.0, 2), C(2, 3), C(3, 4), C(4, 5)) == C(10, 14));

enum Summation: int;

Summation algorithms.

appropriate: Performs pairwise summation for floating point based types and fast summation for integral based types.
pairwise: Pairwise summation algorithm.
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 32bit x86, nextDown(proir) <= result && result <= nextUp(proir). The current implementation re-establish special value semantics across iterations (i.e. handling ±inf).
References "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", Jonathan Richard Shewchuk, Mark Dickinson's post at bugs.python.org.
decimal: Precise decimal summation algorithm.
The elements of the sum are converted to a shortest decimal representation that being converted back would result the same floating-point number. The resulting decimal elements are summed without rounding. The decimal sum is converted back to a binary floating point representation using round-half-to-even rule.

See Also:
The Ryu algorithm
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 && !is(Range : __vector(V[N]), V, size_t N)); void put(Range : Slice!(Iterator, N, kind), Iterator, size_t N, SliceKind kind)(Range r);

Adds n to the internal partial sums.

const scope 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 (isMutable!F) template sum(Summation summation = Summation.appropriate) template sum(F, string summation) if (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 Summation.fast.

Returns:

The sum of all the elements in the range r.

F sum(Range)(Range r) if (isIterable!Range);
F sum(scope const F[] r...);

template fillCollapseSums(Summation summation, alias combineParts, combineElements...)

@property ref auto fillCollapseSums(Iterator, SliceKind kind)(Slice!(Iterator, 1, kind) data);

Library Reference

mir.math.sum