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.
License:
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`
`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).
`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:
`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`);