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

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

assert([1e-20, 1].sum!"decimal" == 1);

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

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

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

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

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

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

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

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

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