fsum

import std.algorithm;
auto ar = [1, 1e100, 1, -1e100].map!(a => a*10000);
const r = 20000;
assert(r == ar.fsum!(Summation.KBN));
assert(r == ar.fsum!(Summation.KB2));
assert(r == ar.fsum); //Summation.Precise

import std.math, std.algorithm, std.range;
auto ar = 1000
    .iota
    .map!(n => 1.7L.pow(n+1) - 1.7L.pow(n))
    .array
    ;
//Summation.Precise is default
real d = 1.7L.pow(1000);
assert(fsum(ar.chain([-d])) == -1);
assert(fsum(-d, ar.retro) == -1);

static struct Quaternion(F)
    if (isFloatingPoint!F)
{
    F[3] array;

    /// + and - operator overloading
    Quaternion opBinary(string op)(auto ref Quaternion rhs) const
        if (op == "+" || op == "-")
    {
        Quaternion ret = void;
        foreach (i, ref e; ret.array)
            mixin("e = array[i] "~op~" rhs.array[i];");
        return ret;
    }

    /// += and -= operator overloading
    Quaternion opOpAssign(string op)(auto ref Quaternion rhs)
        if (op == "+" || op == "-")
    {
        foreach (i, ref e; array)
            mixin("e "~op~"= rhs.array[i];");
        return this;
    }

    ///constructor with single FP argument
    this(F f)
    {
        array[] = f;
    }

    ///assigment with single FP argument
    void opAssign(F f)
    {
        array[] = f;
    }
}

Quaternion!double q, p, r;
q.array = [0, 1, 2];
p.array = [3, 4, 5];
r.array = [3, 5, 7];

assert(r == [p, q].fsum!(Summation.Naive));
assert(r == [p, q].fsum!(Summation.Pairwise));
assert(r == [p, q].fsum!(Summation.Kahan));

import std.complex;
Complex!double[] ar = [complex(1.0, 2), complex(2, 3), complex(3, 4), complex(4, 5)];
Complex!double r = complex(10, 14);
assert(r == ar.fsum!(Summation.Fast));
assert(r == ar.fsum!(Summation.Naive));
assert(r == ar.fsum!(Summation.Pairwise));
assert(r == ar.fsum!(Summation.Kahan));
assert(r == ar.fsum!(Summation.KBN));
assert(r == ar.fsum!(Summation.KB2));
assert(r == ar.fsum); //Summation.Precise