This module contains summation algorithms.

Authors

Ilya Yaroshenko

Copyright

Copyright © 2015-, Ilya Yaroshenko

Example

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");

Example

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);

Example

aive
, 
Pairwise
and 
Kahan` algorithms can be used for user defined types.

import std.traits : 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");

Example

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");

Example

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));

Example

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));

Example

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);

Example

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

Example

Moving mean

import mir.ndslice.topology: linspace;
import mir.math.sum;
import mir.array.allocation: 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);

Example

SIMD Vectors

Bugs

ICE 1662 (dmd only)
import core.simd;
import std.meta : AliasSeq;
double2 a = 1, b = 2, c = 3, d = 6;
with(Summation)
{
    foreach (algo; AliasSeq!(naive, fast, pairwise, kahan))
    {
        assert([a, b, c].sum!algo.array == d.array);
        assert([a, b].sum!algo(c).array == d.array);
    }
}

Enums

Summation

Functions

sumPrecise

Structs

Summator

Templates

sumReturns:

contains detailed documentation and examples about available summation algorithms.

sumReturns:

contains detailed documentation and examples about available summation algorithms.

sumReturns:

contains detailed documentation and examples about available summation algorithms.

sumReturns:

contains detailed documentation and examples about available summation algorithms.

fillCollapseSums