diff --git a/posix.mak b/posix.mak
index c509da368af..13bcd79e250 100644
--- a/posix.mak
+++ b/posix.mak
@@ -91,7 +91,7 @@ DOCSRC = ../dlang.org
 WEBSITE_DIR = ../web
 DOC_OUTPUT_DIR = $(WEBSITE_DIR)/phobos-prerelease
 BIGDOC_OUTPUT_DIR = /tmp
-SRC_DOCUMENTABLES = index.d $(addsuffix .d,$(STD_MODULES) $(STD_NET_MODULES) $(STD_DIGEST_MODULES) $(STD_CONTAINER_MODULES) $(STD_RANGE_MODULES) $(STD_ALGO_MODULES) std/regex/package $(EXTRA_DOCUMENTABLES) $(STD_LOGGER_MODULES))
+SRC_DOCUMENTABLES = index.d $(addsuffix .d,$(STD_MODULES) $(STD_NET_MODULES) $(STD_DIGEST_MODULES) $(STD_CONTAINER_MODULES) $(STD_RANGE_MODULES) $(STD_ALGO_MODULES) std/regex/package $(EXTRA_DOCUMENTABLES) $(STD_LOGGER_MODULES) $(STD_NUMERIC_MODULES))
 STDDOC = $(DOCSRC)/html.ddoc $(DOCSRC)/dlang.org.ddoc $(DOCSRC)/std_navbar-prerelease.ddoc $(DOCSRC)/std.ddoc $(DOCSRC)/macros.ddoc
 BIGSTDDOC = $(DOCSRC)/std_consolidated.ddoc $(DOCSRC)/macros.ddoc
 # Set DDOC, the documentation generator
@@ -186,7 +186,7 @@ STD_MODULES = $(addprefix std/, array ascii base64 bigint \
 		bitmanip compiler complex concurrency conv		\
 		cstream csv datetime demangle encoding exception	\
 		file format functional getopt json math mathspecial	\
-		metastrings mmfile numeric outbuffer parallelism path	\
+		metastrings mmfile outbuffer parallelism path	\
 		process random signals socket socketstream	\
 		stdint stdio stdiobase stream string syserror system traits		\
 		typecons typetuple uni uri utf uuid variant xml zip zlib)
@@ -209,6 +209,8 @@ STD_DIGEST_MODULES = $(addprefix std/digest/, digest crc md ripemd sha)
 STD_CONTAINER_MODULES = $(addprefix std/container/, package array \
 		binaryheap dlist rbtree slist util)
 
+STD_NUMERIC_MODULES = $(addprefix std/numeric/, package summation)
+
 # OS-specific D modules
 EXTRA_MODULES_LINUX := $(addprefix std/c/linux/, linux socket)
 EXTRA_MODULES_OSX := $(addprefix std/c/osx/, socket)
@@ -236,7 +238,8 @@ EXTRA_MODULES += $(EXTRA_DOCUMENTABLES) $(addprefix			\
 # Aggregate all D modules relevant to this build
 D_MODULES = $(STD_MODULES) $(EXTRA_MODULES) $(STD_NET_MODULES) \
 	$(STD_DIGEST_MODULES) $(STD_CONTAINER_MODULES) $(STD_REGEX_MODULES) \
-	$(STD_RANGE_MODULES) $(STD_ALGO_MODULES) $(STD_LOGGER_MODULES)
+	$(STD_RANGE_MODULES) $(STD_ALGO_MODULES) $(STD_LOGGER_MODULES) \
+	$(STD_NUMERIC_MODULES)
 
 # Add the .d suffix to the module names
 D_FILES = $(addsuffix .d,$(D_MODULES))
@@ -442,6 +445,9 @@ $(DOC_OUTPUT_DIR)/std_algorithm_%.html : std/algorithm/%.d $(STDDOC)
 $(DOC_OUTPUT_DIR)/std_range_%.html : std/range/%.d $(STDDOC)
 	$(DDOC) project.ddoc $(STDDOC) -Df$@ $<
 
+$(DOC_OUTPUT_DIR)/std_numeric_%.html : std/numeric/%.d $(STDDOC)
+	$(DDOC) project.ddoc $(STDDOC) -Df$@ $<
+
 $(DOC_OUTPUT_DIR)/std_regex_%.html : std/regex/%.d $(STDDOC)
 	$(DDOC) project.ddoc $(STDDOC) -Df$@ $<
 
diff --git a/std/numeric.d b/std/numeric/package.d
similarity index 98%
rename from std/numeric.d
rename to std/numeric/package.d
index 4c059e1dc17..32b77b9f2b7 100644
--- a/std/numeric.d
+++ b/std/numeric/package.d
@@ -23,6 +23,8 @@ Distributed under the Boost Software License, Version 1.0.
 */
 module std.numeric;
 
+public import std.numeric.summation;
+
 import std.complex;
 import std.exception;
 import std.math;
@@ -812,6 +814,7 @@ T findRoot(T, DF, DT)(scope DF f, in T a, in T b,
         is(typeof(f(T.init)) == R, R) && isFloatingPoint!R
     )
 {
+    import std.math : fabs; // FIXME
     auto r = findRoot(f, a, b, f(a), f(b), tolerance);
     // Return the first value if it is smaller or NaN
     return !(fabs(r[2]) > fabs(r[3])) ? r[0] : r[1];
@@ -870,6 +873,8 @@ in
 }
 body
 {
+    import std.math : fabs; // FIXME
+
     // Author: Don Clugston. This code is (heavily) modified from TOMS748
     // (www.netlib.org).  The changes to improve the worst-cast performance are
     // entirely original.
@@ -1583,6 +1588,7 @@ unittest
     }
 }
 
+
 /**
 Normalizes values in $(D range) by multiplying each element with a
 number chosen such that values sum up to $(D sum). If elements in $(D
@@ -1646,45 +1652,6 @@ unittest
     assert(a == [ 0.5, 0.5 ]);
 }
 
-/**
-Computes accurate sum of binary logarithms of input range $(D r).
- */
-ElementType!Range sumOfLog2s(Range)(Range r)
-    if (isInputRange!Range && isFloatingPoint!(ElementType!Range))
-{
-    long exp = 0;
-    Unqual!(typeof(return)) x = 1;
-    foreach (e; r)
-    {
-        if (e < 0)
-            return typeof(return).nan;
-        int lexp = void;
-        x *= frexp(e, lexp);
-        exp += lexp;
-        if (x < 0.5)
-        {
-            x *= 2;
-            exp--;
-        }
-    }
-    return exp + log2(x);
-}
-
-///
-unittest
-{
-    assert(sumOfLog2s(new double[0]) == 0);
-    assert(sumOfLog2s([0.0L]) == -real.infinity);
-    assert(sumOfLog2s([-0.0L]) == -real.infinity);
-    assert(sumOfLog2s([2.0L]) == 1);
-    assert(sumOfLog2s([-2.0L]).isNaN());
-    assert(sumOfLog2s([real.nan]).isNaN());
-    assert(sumOfLog2s([-real.nan]).isNaN());
-    assert(sumOfLog2s([real.infinity]) == real.infinity);
-    assert(sumOfLog2s([-real.infinity]).isNaN());
-    assert(sumOfLog2s([ 0.25, 0.25, 0.25, 0.125 ]) == -9);
-}
-
 /**
 Computes $(LUCKY _entropy) of input range $(D r) in bits. This
 function assumes (without checking) that the values in $(D r) are all
diff --git a/std/numeric/summation.d b/std/numeric/summation.d
new file mode 100644
index 00000000000..add9b33a503
--- /dev/null
+++ b/std/numeric/summation.d
@@ -0,0 +1,1252 @@
+/**
+This module contains basic summation algorithms.
+
+License: $(LINK2 http://boost.org/LICENSE_1_0.txt, Boost License 1.0).
+
+Authors: $(WEB 9il.github.io, Ilya Yaroshenko)
+
+Source: $(PHOBOSSRC std/numeric/_summation.d)
+*/
+module std.numeric.summation;
+
+import std.traits;
+import std.typecons;
+import std.range.primitives;
+import std.math;
+
+private template SummationType(F)
+    if (isFloatingPoint!F || isComplex!F)
+{
+    version(X86) //workaround for Issue 13474
+    {
+        static if(!is(Unqual!F == real) && (isComplex!F || !is(Unqual!(typeof(F.init.re)) == real)))
+            pragma(msg, "Warning: Summation algorithms on x86 use 80bit representation for single and double floating point numbers.");
+        static if(isComplex!F)
+        {
+            import std.complex : Complex;
+            alias SummationType = Complex!real;
+        }
+        else
+            alias SummationType = real;
+    }
+    else
+        alias SummationType = F;
+}
+
+/++
+Computes accurate sum of binary logarithms of input range $(D r).
++/
+ElementType!Range sumOfLog2s(Range)(Range r)
+    if (isInputRange!Range && isFloatingPoint!(ElementType!Range))
+{
+    long exp = 0;
+    Unqual!(typeof(return)) x = 1;
+    foreach (e; r)
+    {
+        if (e < 0)
+            return typeof(return).nan;
+        int lexp = void;
+        x *= frexp(e, lexp);
+        exp += lexp;
+        if (x < 0.5)
+        {
+            x *= 2;
+            exp--;
+        }
+    }
+    return exp + log2(x);
+}
+
+///
+unittest
+{
+    import std.math, std.numeric;
+    assert(sumOfLog2s(new double[0]) == 0);
+    assert(sumOfLog2s([0.0L]) == -real.infinity);
+    assert(sumOfLog2s([-0.0L]) == -real.infinity);
+    assert(sumOfLog2s([2.0L]) == 1);
+    assert(sumOfLog2s([-2.0L]).isNaN());
+    assert(sumOfLog2s([real.nan]).isNaN());
+    assert(sumOfLog2s([-real.nan]).isNaN());
+    assert(sumOfLog2s([real.infinity]) == real.infinity);
+    assert(sumOfLog2s([-real.infinity]).isNaN());
+    assert(sumOfLog2s([ 0.25, 0.25, 0.25, 0.125 ]) == -9);
+}
+
+
+/++
+Summation algorithms.
++/
+enum Summation
+{
+    /++
+    Fast summation algorithm.
+    +/
+    Fast,
+
+    /++
+    Naive algorithm (one by one).
+    +/
+    Naive,
+
+    /++
+    $(LUCKY Pairwise summation) algorithm. Range must be a finite sliceable range.
+    +/
+    Pairwise,
+
+    /++
+    $(LUCKY Kahan summation) algorithm.
+    +/
+    Kahan,
+
+    /++
+    $(LUCKY Kahan-Babuška-Neumaier summation algorithm). $(D КBN) gives more accurate results then $(D Kahan).
+    +/
+    KBN,
+
+    /++
+    $(LUCKY Generalized Kahan-Babuška summation algorithm), order 2. $(D КB2) gives more accurate results then $(D Kahan) and $(D КBN).
+    +/
+    KB2,
+
+    /++
+    Precise summation algorithm.
+    The value of the sum is rounded to the nearest representable
+    floating-point number using the $(LUCKY round-half-to-even rule).
+    Result can be differ from the exact value on $(D X86), $(D nextDown(proir) <= result &&  result <= nextUp(proir)).
+    +/
+    Precise,
+}
+
+
+/++
+Computes sum of range.
++/
+template fsum(F, Summation summation = Summation.Precise)
+    if (isFloatingPoint!F && isMutable!F)
+{
+    F fsum(Range)(Range r)
+        if (isSummable!(Range, F))
+    {
+        alias sum = Algo!summation;
+        return sum!(Range, F)(r);
+    }
+
+    F fsum(Range)(F seed, Range r)
+        if (isSummable!(Range, F))
+    {
+        alias sum = Algo!summation;
+        return sum!(Range, F)(r, seed);
+    }
+}
+
+///ditto
+template fsum(Summation summation = Summation.Precise)
+{
+    Unqual!(ForeachType!Range) fsum(Range)(Range r)
+        if (isSummable!(Range, Unqual!(ForeachType!Range)))
+    {
+        alias sum = Algo!summation;
+        return sum!(Range, typeof(return))(r);
+    }
+
+    F fsum(F, Range)(F seed, Range r)
+        if (isSummable!(Range, F))
+    {
+        alias sum = Algo!summation;
+        return sum!(Range, F)(r, seed);
+    }
+}
+
+///
+unittest {
+    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
+}
+
+// FIXME
+// Fails for 32bit systems.
+// See also https://issues.dlang.org/show_bug.cgi?id=13474#c7
+// and https://github.com/D-Programming-Language/phobos/pull/2513
+version(X86)
+{
+
+}
+else
+{
+    unittest {
+        import std.algorithm;
+        auto ar = [1, 1e100, 1, -1e100].map!(a => a*10000);
+        const r = 20000;
+        assert(r != ar.fsum!(Summation.Naive));
+        assert(r != ar.fsum!(Summation.Kahan));
+        assert(r != ar.fsum!(Summation.Pairwise));
+    }
+}
+
+///
+unittest
+{
+    import core.stdc.tgmath, std.algorithm, std.range;
+    auto ar = 1000.0
+        .iota
+        .map!(n => 1.7.pow(n+1.0) - 1.7.pow(n))
+        .array
+        ;
+    //Summation.Precise is default
+    double d = 1.7.pow(1000.0);
+    assert(fsum(ar.chain([-d]))  ==  -1.0);
+    assert(fsum!real(-d, ar.retro)  ==  -1.0);
+}
+
+/++
+$(D Naive), $(D Pairwise) and $(D Kahan) algorithms can be used for user defined types.
++/
+unittest {
+    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 == "-")
+        {
+            Quaternion ret = void;
+            foreach (i, ref e; array)
+                mixin("e "~op~"= rhs.array[i];");
+            return this;
+        }
+
+        ///constructor with single FP argument
+        this(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));
+}
+
+/++
+All summation algorithms available for complex numbers.
++/
+unittest
+{
+    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
+}
+
+
+/++
+Handler for full precise summation with $(D put) primitive.
+The current implementation re-establish special
+value semantics across iterations (i.e. handling -inf + inf).
+
+References: $(LINK2 http://www.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps,
+    "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", Jonathan Richard Shewchuk),
+    $(LINK2 http://bugs.python.org/file10357/msum4.py, Mark Dickinson's post at bugs.python.org).
++/
+/+
+Precise summation function as msum() by Raymond Hettinger in
+<http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
+enhanced with the exact partials sum and roundoff from Mark
+Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
+See those links for more details, proofs and other references.
+IEEE 754R floating point semantics are assumed.
++/
+struct Summator(T)
+    if (isFloatingPoint!T && isMutable!T)
+{
+    import std.internal.scopebuffer;
+
+private:
+    alias F = SummationType!T;
+    enum F M = (cast(F)(2)) ^^ (T.max_exp - 1);
+    F[32] scopeBufferArray = void;
+    ScopeBuffer!F partials;
+    //sum for NaN and infinity.
+    F s;
+    //Overflow Degree. Count of 2^^F.max_exp minus count of -(2^^F.max_exp)
+    sizediff_t o;
+
+
+    /++
+    Compute the sum of a list of nonoverlapping floats.
+    On input, partials is a list of nonzero, nonspecial,
+    nonoverlapping floats, strictly increasing in magnitude, but
+    possibly not all having the same sign.
+    On output, the sum of partials gives the error in the returned
+    result, which is correctly rounded (using the round-half-to-even
+    rule).
+    Two floating point values x and y are non-overlapping if the least significant nonzero
+    bit of x is more significant than the most significant nonzero bit of y, or vice-versa.
+    +/
+    static F partialsReduce(F s, in F[] partials)
+    in
+    {
+        debug(numeric) assert(!partials.length || s.isFinite);
+    }
+    body
+    {
+        bool _break = void;
+        foreach_reverse(i, y; partials)
+        {
+            s = partialsReducePred(s, y, i ? partials[i-1] : 0, _break);
+            if (_break)
+                break;
+            debug(numeric) assert(s.isFinite);
+        }
+        return s;
+    }
+
+    static F partialsReducePred(F s, F y, F z, out bool _break)
+    out(result)
+    {
+        debug(numeric) assert(result.isFinite);
+    }
+    body
+    {
+        F x = s;
+        s = x + y;
+        F d = s - x;
+        F l = y - d;
+        debug(numeric)
+        {
+            assert(x.isFinite);
+            assert(y.isFinite);
+            assert(s.isFinite);
+            assert(fabs(y) < fabs(x));
+        }
+        if (l)
+        {
+        //Make half-even rounding work across multiple partials.
+        //Needed so that sum([1e-16, 1, 1e16]) will round-up the last
+        //digit to two instead of down to zero (the 1e-16 makes the 1
+        //slightly closer to two). Can guarantee commutativity.
+            if (z && !signbit(l * z))
+            {
+                l *= 2;
+                x = s + l;
+                F t = x - s;
+                if (l == t)
+                    s = x;
+            }
+            _break = true;
+        }
+        return s;
+    }
+
+    //Returns corresponding infinity if is overflow and 0 otherwise.
+    F overflow() const
+    {
+        if (o == 0)
+            return 0;
+        if (partials.length && (o == -1 || o == 1)  && signbit(o * partials[$-1]))
+        {
+            // problem case: decide whether result is representable
+            F x = o * M;
+            F y = partials[$-1] / 2;
+            F h = x + y;
+            F d = h - x;
+            F l = (y - d) * 2;
+            y = h * 2;
+            d = h + l;
+            F t = d - h;
+            version(X86)
+            {
+                if (!.isInfinity(cast(T)y) || !.isInfinity(sum()))
+                    return 0;
+            }
+            else
+            {
+                if (!.isInfinity(cast(T)y) || partials.length > 1 && !signbit(l * partials[$-2]) && t == l)
+                    return 0;
+            }
+        }
+        return F.infinity * o;
+    }
+
+public:
+
+    ///
+    this(F x)
+    {
+        partials = scopeBuffer(scopeBufferArray);
+        s = 0;
+        o = 0;
+        if (x) put(x);
+    }
+
+    ///
+    @disable this();
+
+    // free ScopeBuffer
+    ~this()
+    {
+        partials.free;
+    }
+
+    // copy ScopeBuffer if necessary
+    this(this)
+    {
+        auto a = partials[];
+        if (scopeBufferArray.ptr !is a.ptr)
+        {
+            partials = scopeBuffer(scopeBufferArray);
+            partials.put(a);
+        }
+    }
+
+    ///Adds $(D x) to the internal partial sums.
+    void put(T _x)
+    {
+        F x = _x;
+        if (.isFinite(x))
+        {
+            size_t i;
+            foreach (y; partials[])
+            {
+                F h = x + y;
+                if (.isInfinity(cast(T)h))
+                {
+                    if (fabs(x) < fabs(y))
+                    {
+                        F t = x; x = y; y = t;
+                    }
+                    //h == -F.infinity
+                    if (signbit(h))
+                    {
+                        x += M;
+                        x += M;
+                        o--;
+                    }
+                    //h == +F.infinity
+                    else
+                    {
+                        x -= M;
+                        x -= M;
+                        o++;
+                    }
+                    debug(numeric) assert(x.isFinite);
+                    h = x + y;
+                }
+                debug(numeric) assert(h.isFinite);
+                F l;
+                if(fabs(x) < fabs(y))
+                {
+                    F t = h - y;
+                    l = x - t;
+                }
+                else
+                {
+                    F t = h - x;
+                    l = y - t;
+                }
+                debug(numeric) assert(l.isFinite);
+                if (l)
+                {
+                    partials[i++] = l;
+                }
+                x = h;
+            }
+            partials.length = i;
+            if (x)
+            {
+                partials.put(x);
+            }
+        }
+        else
+        {
+            s += x;
+        }
+    }
+
+    /++
+    Adds $(D x) to the internal partial sums.
+    This operation doesn't re-establish special
+    value semantics across iterations (i.e. handling -inf + inf).
+    Preconditions: $(D isFinite(x)).
+    +/
+    package void unsafePut(F x)
+    in {
+        assert(.isFinite(x));
+    }
+    body {
+        size_t i;
+        foreach (y; partials[])
+        {
+            F h = x + y;
+            debug(numeric) assert(.isFinite(h));
+            F l;
+            if(fabs(x) < fabs(y))
+            {
+                F t = h - y;
+                l = x - t;
+            }
+            else
+            {
+                F t = h - x;
+                l = y - t;
+            }
+            debug(numeric) assert(.isFinite(l));
+            if (l)
+            {
+                partials[i++] = l;
+            }
+            x = h;
+        }
+        partials.length = i;
+        if (x)
+        {
+            partials.put(x);
+        }
+    }
+
+    ///
+    unittest {
+        import core.stdc.tgmath, std.algorithm, std.range;
+        auto r = iota(1000.0).map!(a => 1.7.pow(a+1) - 1.7.pow(a));
+        Summator!double s = 0;
+        put(s, r);
+        s -= 1.7.pow(1000.0);
+        assert(s.sum() == -1);
+    }
+
+    /++
+    Returns the value of the sum, rounded to the nearest representable
+    floating-point number using the round-half-to-even rule.
+    Result can be differ from the exact value on $(D X86), $(D nextDown(proir) <= result &&  result <= nextUp(proir)).
+    +/
+    T sum() const
+    {
+        debug(numeric)
+        {
+            foreach (y; partials[])
+            {
+                assert(y);
+                assert(y.isFinite);
+            }
+            //TODO: Add Non-Overlapping check to std.math
+            import std.algorithm : isSorted, map;
+            assert(partials[].map!(a => fabs(a)).isSorted);
+        }
+
+        if (s)
+            return s;
+        auto parts = partials[];
+        F y = 0;
+        //pick last
+        if (parts.length)
+        {
+            y = parts[$-1];
+            parts = parts[0..$-1];
+        }
+        if (o)
+        {
+            immutable F of = o;
+            if (y && (o == -1 || o == 1)  && signbit(of * y))
+            {
+                // problem case: decide whether result is representable
+                y /= 2;
+                F x = of * M;
+                immutable F h = x + y;
+                F t = h - x;
+                F l = (y - t) * 2;
+                y = h * 2;
+                if (.isInfinity(cast(T)y))
+                {
+                    // overflow, except in edge case...
+                    x = h + l;
+                    t = x - h;
+                    y = parts.length && t == l && !signbit(l*parts[$-1]) ?
+                        x * 2 :
+                        F.infinity * of;
+                    parts = null;
+                }
+                else if (l)
+                {
+                    bool _break;
+                    y = partialsReducePred(y, l, parts.length ? parts[$-1] : 0, _break);
+                    if (_break)
+                        parts = null;
+                }
+            }
+            else
+            {
+                y = F.infinity * of;
+                parts = null;
+            }
+        }
+        return partialsReduce(y, parts);
+    }
+
+    version(none)
+    F partialsSum() const
+    {
+        debug(numeric) partialsDebug;
+        auto parts = partials[];
+        F y = 0;
+        //pick last
+        if (parts.length)
+        {
+            y = parts[$-1];
+            parts = parts[0..$-1];
+        }
+        return partialsReduce(y, parts);
+    }
+
+    ///Returns $(D Summator) with extended internal partial sums.
+    C opCast(C : Summator!P, P)()
+        if (
+            isMutable!C &&
+            P.max_exp >= T.max_exp &&
+            P.mant_dig >= T.mant_dig
+        )
+    {
+        static if (is(P == T))
+            return this;
+        else
+        {
+            typeof(return) ret = void;
+            ret.s = s;
+            ret.o = o;
+            ret.partials = scopeBuffer(ret.scopeBufferArray);
+            foreach (p; partials[])
+            {
+                ret.partials.put(p);
+            }
+            enum exp_diff = P.max_exp / T.max_exp;
+            static if (exp_diff)
+            {
+                if (ret.o)
+                {
+                    immutable f = ret.o / exp_diff;
+                    immutable t = cast(int)(ret.o % exp_diff);
+                    ret.o = f;
+                    ret.put((P(2) ^^ T.max_exp) * t);
+                }
+            }
+            return ret;
+        }
+    }
+
+    ///
+    unittest
+    {
+        import std.math;
+        float  M = 2.0f ^^ (float.max_exp-1);
+        double N = 2.0  ^^ (float.max_exp-1);
+        auto s = Summator!float(0); //float summator
+        s += M;
+        s += M;
+        assert(float.infinity == s.sum());
+        auto e = cast(Summator!double) s;
+        assert(isFinite(e.sum()));
+        assert(N+N == e.sum());
+    }
+
+    /++
+    $(D cast(F)) operator overloading. Returns $(D cast(T)sum()).
+    See also: $(D cast)
+    +/
+    F opCast(C)() if (is(Unqual!C == T))
+    {
+        return sum();
+    }
+
+    ///The assignment operator $(D =) overloading.
+    void opAssign(F rhs)
+    {
+        partials.length = 0;
+        s = 0;
+        o = 0;
+        if (rhs) put(rhs);
+    }
+
+    /// $(D +=) and $(D -=) operator overloading.
+    void opOpAssign(string op : "+")(F f)
+    {
+        put(f);
+    }
+
+    ///ditto
+    void opOpAssign(string op : "+")(ref const Summator rhs)
+    {
+        s += rhs.s;
+        o += rhs.o;
+        foreach (f; rhs.partials[])
+            put(f);
+    }
+
+    ///ditto
+    void opOpAssign(string op : "-")(F f)
+    {
+        put(-f);
+    }
+
+    ///ditto
+    void opOpAssign(string op : "-")(ref const Summator rhs)
+    {
+        s -= rhs.s;
+        o -= rhs.o;
+        foreach (f; rhs.partials[])
+            put(-f);
+    }
+
+    ///
+    unittest {
+        import core.stdc.tgmath, std.algorithm, std.range;
+        auto r1 = iota(500.0).map!(a => 1.7.pow(a+1) - 1.7.pow(a));
+        auto r2 = iota(500.0, 1000.0).map!(a => 1.7.pow(a+1) - 1.7.pow(a));
+        Summator!double s1 = 0.0, s2 = 0.0;
+        foreach (e; r1) s1 += e;
+        foreach (e; r2) s2 -= e;
+        s1 -= s2;
+        s1 -= 1.7.pow(1000.0);
+        assert(s1.sum() == -1);
+    }
+
+    ///Returns $(D true) if current sum is a NaN.
+    bool isNaN() const
+    {
+        return .isNaN(s);
+    }
+
+    ///Returns $(D true) if current sum is finite (not infinite or NaN).
+    bool isFinite() const
+    {
+        if (s)
+            return false;
+        return !overflow;
+    }
+
+    ///Returns $(D true) if current sum is ±∞.
+    bool isInfinity() const
+    {
+        return .isInfinity(s) || overflow();
+    }
+}
+
+unittest
+{
+    import std.range;
+    import std.algorithm;
+    import core.stdc.tgmath;
+
+    Summator!double summator = 0;
+
+    enum double M = (cast(double)2) ^^ (double.max_exp - 1);
+    Tuple!(double[], double)[] tests = [
+        tuple(new double[0], 0.0),
+        tuple([0.0], 0.0),
+        tuple([1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100),
+        tuple([1e308, 1e308, -1e308], 1e308),
+        tuple([-1e308, 1e308, 1e308], 1e308),
+        tuple([1e308, -1e308, 1e308], 1e308),
+        tuple([M, M, -2.0^^1000], 1.7976930277114552e+308),
+        tuple([M, M, M, M, -M, -M, -M], 8.9884656743115795e+307),
+        tuple([2.0^^53, -0.5, -2.0^^-54], 2.0^^53-1.0),
+        tuple([2.0^^53, 1.0, 2.0^^-100], 2.0^^53+2.0),
+        tuple([2.0^^53+10.0, 1.0, 2.0^^-100], 2.0^^53+12.0),
+        tuple([2.0^^53-4.0, 0.5, 2.0^^-54], 2.0^^53-3.0),
+        tuple([M-2.0^^970, -1.0, M], 1.7976931348623157e+308),
+        tuple([double.max, double.max*2.^^-54], double.max),
+        tuple([double.max, double.max*2.^^-53], double.infinity),
+        tuple(iota(1, 1001).map!(a => 1.0/a).array , 7.4854708605503451),
+        tuple(iota(1, 1001).map!(a => (-1.0)^^a/a).array, -0.69264743055982025), //0.693147180559945309417232121458176568075500134360255254120680...
+        tuple(iota(1000.0).map!(a => 1.7.pow(a+1) - 1.7.pow(a)).chain([-(1.7.pow(1000.0))]).array , -1.0),
+        tuple(iota(1, 1001).map!(a => 1.0/a).retro.array , 7.4854708605503451),
+        tuple(iota(1, 1001).map!(a => (-1.0)^^a/a).retro.array, -0.69264743055982025),
+        tuple(iota(1000.0).map!(a => 1.7.pow(a+1) - 1.7.pow(a)).chain([-(1.7.pow(1000.0))]).retro.array , -1.0),
+        tuple([double.infinity, -double.infinity, double.nan], double.nan),
+        tuple([double.nan, double.infinity, -double.infinity], double.nan),
+        tuple([double.infinity, double.nan, double.infinity], double.nan),
+        tuple([double.infinity, double.infinity], double.infinity),
+        tuple([double.infinity, -double.infinity], double.nan),
+        tuple([-double.infinity, 1e308, 1e308, -double.infinity], -double.infinity),
+        tuple([M-2.0^^970, 0.0, M], double.infinity),
+        tuple([M-2.0^^970, 1.0, M], double.infinity),
+        tuple([M, M], double.infinity),
+        tuple([M, M, -1.0], double.infinity),
+        tuple([M, M, M, M, -M, -M], double.infinity),
+        tuple([M, M, M, M, -M, M], double.infinity),
+        tuple([-M, -M, -M, -M], -double.infinity),
+        tuple([M, M, -2.^^971], double.max),
+        tuple([M, M, -2.^^970], double.infinity),
+        tuple([-2.^^970, M, M, -2.^^-1074], double.max),
+        tuple([M, M, -2.^^970, 2.^^-1074], double.infinity),
+        tuple([-M, 2.^^971, -M], -double.max),
+        tuple([-M, -M, 2.^^970], -double.infinity),
+        tuple([-M, -M, 2.^^970, 2.^^-1074], -double.max),
+        tuple([-2.^^-1074, -M, -M, 2.^^970], -double.infinity),
+        tuple([2.^^930, -2.^^980, M, M, M, -M], 1.7976931348622137e+308),
+        tuple([M, M, -1e307], 1.6976931348623159e+308),
+        tuple([1e16, 1., 1e-16], 10000000000000002.0),
+    ];
+    foreach (i, test; tests)
+    {
+        foreach (t; test[0]) summator.put(t);
+        auto r = test[1];
+        auto s = summator.sum;
+        version(X86)
+        {
+            assert(summator.isNaN() == r.isNaN());
+            assert(summator.isFinite() == r.isFinite() || r == -double.max && s == -double.infinity || r == double.max && s == double.infinity);
+            assert(summator.isInfinity() == r.isInfinity() || r == -double.max && s == -double.infinity || r == double.max && s == double.infinity);
+            assert(nextDown(s) <= r && nextUp(s) >= r || s.isNaN && r.isNaN);            
+        }
+        else
+        {
+            assert(summator.isNaN() == r.isNaN());
+            assert(summator.isFinite() == r.isFinite());
+            assert(summator.isInfinity() == r.isInfinity());
+            assert(s == r || s.isNaN && r.isNaN);                        
+        }
+        summator = 0;
+    }
+}
+
+
+private:
+
+template isComplex(C)
+{
+    import std.complex : Complex;
+    enum bool isComplex = is(C : Complex!F, F);
+}
+
+version(X86)
+{
+
+}
+else
+{
+    // FIXME (perfomance issue): fabs in std.math available only for for real.
+    F fabs(F)(F f) //+-0, +-NaN, +-inf doesn't matter
+    {
+        if (__ctfe)
+        {
+            return f < 0 ? -f : f;
+        }
+        else
+        {
+            version(LDC)
+            {
+                import ldc.intrinsics : llvm_fabs;
+                return llvm_fabs(f);
+            }
+            else
+            {
+                import core.stdc.tgmath : fabs;
+                return fabs(f);
+            }
+        }
+    }    
+}
+
+
+template isSummable(Range, F)
+{
+    enum bool isSummable =
+        isInputRange!Range &&
+        isImplicitlyConvertible!(Unqual!(ForeachType!Range), F) &&
+        !isInfinite!Range &&
+        __traits(compiles,
+        {
+            F a = 0.1, b, c;
+            c = a + b;
+            c = a - b;
+            a += b;
+            a -= b;
+        });
+}
+
+/++
+Naive summation algorithm.
++/
+F sumNaive(Range, F = Unqual!(ForeachType!Range))(Range r, F s = 0)
+{
+    foreach (x; r)
+    {
+        s += x;
+    }
+    return s;
+}
+
+///TODO
+alias sumFast = sumNaive;
+
+/++
+$(LUCKY Pairwise summation) algorithm. Range must be a finite sliceable range.
++/
+F sumPairwise(Range, F = Unqual!(ForeachType!Range))(Range r)
+    if (hasLength!Range && hasSlicing!Range)
+{
+    import std.range : hasLength, hasSlicing;
+    static assert(hasLength!Range && hasSlicing!Range);
+    switch (r.length)
+    {
+        case 0: return F(0);
+        case 1: return cast(F)r[0];
+        case 2: return cast(F)(r[0] + cast(F)r[1]);
+        default: auto n = r.length/2; return cast(F)(sumPairwise!(Range, F)(r[0..n]) + sumPairwise!(Range, F)(r[n..$]));
+    }
+}
+
+F sumPairwise(Range, F = Unqual!(ForeachType!Range))(Range r, F seed)
+{
+    F s = seed;
+    s += sumPairwise!Range(r);
+    return s;
+}
+
+
+/++
+$(LUCKY Kahan summation) algorithm.
++/
+/++
+---------------------
+s := x[1]
+c := 0
+FOR k := 2 TO n DO
+y := x[k] - c
+t := s + y
+c := (t - s) - y
+s := t
+END DO
+---------------------
++/
+F sumKahan(Range, F = Unqual!(ForeachType!Range))(Range r, F s = 0)
+{
+    F c = 0.0;
+    F y; // do not declare in the loop (algo can be used for matrixes and etc)
+    F t; // ditto
+    foreach (F x; r)
+    {
+        y = x - c;
+        t = s + y;
+        c = t - s;
+        c -= y;
+        s = t;
+    }
+    return s;
+}
+
+
+/++
+$(LUCKY Kahan-Babuška-Neumaier summation algorithm).
+$(D КBN) gives more accurate results then $(D Kahan).
++/
+/++
+---------------------
+s := x[1]
+c := 0
+FOR i := 2 TO n DO
+t := s + x[i]
+IF ABS(s) >= ABS(x[i]) THEN
+    c := c + ((s-t)+x[i])
+ELSE
+    c := c + ((x[i]-t)+s)
+END IF
+s := t
+END DO
+s := s + c
+---------------------
++/
+T sumKBN(Range, T = Unqual!(ForeachType!Range))(Range r, T _s = 0)
+    if (isFloatingPoint!T || isComplex!T)
+{
+    alias F = SummationType!T;
+    F s = _s;
+    F c = 0.0;
+    static if (isFloatingPoint!F)
+    {
+        foreach (F x; r)
+        {
+            F t = s + x;
+            if (fabs(s) >= fabs(x))
+            {
+                F d = s - t;
+                d += x;
+                c += d;
+            }
+            else
+            {
+                F d = x - t;
+                d += s;
+                c += d;
+            }
+            s = t;
+        }
+    }
+    else
+    {
+        foreach (F x; r)
+        {
+            F t = s + x;
+            if (fabs(s.re) < fabs(x.re))
+            {
+                auto t_re = s.re;
+                s.re = x.re;
+                x.re = t_re;
+            }
+            if (fabs(s.im) < fabs(x.im))
+            {
+                auto t_im = s.im;
+                s.im = x.im;
+                x.im = t_im;
+            }
+            F d = s - t;
+            d += x;
+            c += d;
+            s = t;
+        }
+    }
+    s += c;
+    return cast(T)s;
+}
+
+
+/++
+$(LUCKY Generalized Kahan-Babuška summation algorithm), order 2.
+$(D КB2) gives more accurate results then $(D Kahan) and $(D КBN).
++/
+/++
+---------------------
+s := 0 ; cs := 0 ; ccs := 0
+FOR j := 1 TO n DO
+    t := s + x[i]
+    IF ABS(s) >= ABS(x[i]) THEN
+        c := (s-t) + x[i]
+    ELSE
+        c := (x[i]-t) + s
+    END IF
+    s := t
+    t := cs + c
+    IF ABS(cs) >= ABS(c) THEN
+        cc := (cs-t) + c
+    ELSE
+        cc := (c-t) + cs
+    END IF
+    cs := t
+    ccs := ccs + cc
+END FOR
+RETURN s+cs+ccs
+---------------------
++/
+T sumKB2(Range, T = Unqual!(ForeachType!Range))(Range r, T _s = 0)
+    if (isFloatingPoint!T || isComplex!T)
+{
+    alias F = SummationType!T;
+    F s = _s;
+    F cs = 0.0;
+    F ccs = 0.0;
+    static if (isFloatingPoint!F)
+    {
+        foreach (F x; r)
+        {
+            F t = s + x;
+            F c = void;
+            if (fabs(s) >= fabs(x))
+            {
+                F d = s - t;
+                c = d + x;
+            }
+            else
+            {
+                F d = x - t;
+                c = d + s;
+            }
+            s = t;
+            t = cs + c;
+            if (fabs(cs) >= fabs(c))
+            {
+                F d = cs - t;
+                d += c;
+                ccs += d;
+            }
+            else
+            {
+                F d = c - t;
+                d += cs;
+                ccs += d;
+            }
+            cs = t;
+        }
+    }
+    else
+    {
+        foreach (F x; r)
+        {
+            F t = s + x;
+            if (fabs(s.re) < fabs(x.re))
+            {
+                auto t_re = s.re;
+                s.re = x.re;
+                x.re = t_re;
+            }
+            if (fabs(s.im) < fabs(x.im))
+            {
+                auto t_im = s.im;
+                s.im = x.im;
+                x.im = t_im;
+            }
+            F c = (s-t)+x;
+            s = t;
+            if (fabs(cs.re) < fabs(c.re))
+            {
+                auto t_re = cs.re;
+                cs.re = c.re;
+                c.re = t_re;
+            }
+            if (fabs(cs.im) < fabs(c.im))
+            {
+                auto t_im = cs.im;
+                cs.im = c.im;
+                c.im = t_im;
+            }
+            F d = cs - t;
+            d += c;
+            ccs += d;
+            cs = t;
+        }
+    }
+    cs += ccs;
+    s += cs;
+    return cast(T)s;
+}
+
+unittest
+{
+    import std.typetuple;
+    with(Summation)
+    foreach (F; TypeTuple!(float, double, real))
+    {
+        F[] ar = [1, 2, 3, 4];
+        F r = 10;
+        assert(r == ar.fsum!Fast());
+        assert(r == ar.fsum!Pairwise());
+        assert(r == ar.fsum!Kahan());
+        assert(r == ar.fsum!KBN());
+        assert(r == ar.fsum!KB2());
+    }
+}
+
+//@@@BUG@@@: DMD 2.066 Segmentation fault (core dumped)
+version(none)
+unittest
+{
+    import core.simd;
+    static if (__traits(compiles, double2.init + double2.init))
+    {
+        double2[] ar = [double2([1.0, 2]), double2([2, 3]), double2([3, 4]), double2([4, 6])];
+        assert(ar.sumFast().array == double2([10, 14]).array);
+        assert(ar.sumPairwise().array == double2([10, 14]).array);
+        assert(ar.sumKahan().array == double2([10, 14]).array);
+    }
+}
+
+
+/++
+Precise summation.
++/
+F sumPrecise(Range, F = Unqual!(ForeachType!Range))(Range r, F seed = 0)
+    if (isFloatingPoint!F || isComplex!F)
+{
+    static if (isFloatingPoint!F)
+    {
+        auto sum = Summator!F(seed);
+        foreach (e; r)
+        {
+            sum.put(e);
+        }
+        return sum.sum;
+    }
+    else
+    {
+        alias T = typeof(F.init.re);
+        static if (isForwardRange!Range)
+        {
+            auto s = r.save;
+            auto sum = Summator!T(seed.re);
+            foreach (e; r)
+            {
+                sum.put(e.re);
+            }
+            T sumRe = sum.sum;
+            sum = seed.im;
+            foreach (e; s)
+            {
+                sum.put(e.im);
+            }
+            return F(sumRe, sum.sum);
+        }
+        else
+        {
+            auto sumRe = Summator!T(seed.re);
+            auto sumIm = Summator!T(seed.im);
+            foreach (e; r)
+            {
+                sumRe.put(e.re);
+                sumIm.put(e.im);
+            }
+            return F(sumRe.sum, sumIm.sum);
+        }
+    }
+}
+
+template Algo(Summation summation)
+{
+    static if (summation == Summation.Fast)
+        alias Algo = sumFast;
+    else
+    static if (summation == Summation.Naive)
+        alias Algo = sumNaive;
+    else
+    static if (summation == Summation.Pairwise)
+        alias Algo = sumPairwise;
+    else
+    static if (summation == Summation.Kahan)
+        alias Algo = sumKahan;
+    else
+    static if (summation == Summation.KBN)
+        alias Algo = sumKBN;
+    else
+    static if (summation == Summation.KB2)
+        alias Algo = sumKB2;
+    else
+    static if (summation == Summation.Precise)
+        alias Algo = sumPrecise;
+    else
+    static assert(0);
+}
diff --git a/win32.mak b/win32.mak
index e2fefa38a01..7e9282a6660 100644
--- a/win32.mak
+++ b/win32.mak
@@ -110,7 +110,7 @@ SRC_STD_1_HEAVY= std\stdio.d std\stdiobase.d \
 
 SRC_STD_2a_HEAVY= std\array.d std\functional.d std\path.d std\outbuffer.d std\utf.d
 
-SRC_STD_3= std\csv.d std\math.d std\complex.d std\numeric.d std\bigint.d \
+SRC_STD_3= std\csv.d std\math.d std\complex.d std\bigint.d \
 	std\bitmanip.d std\typecons.d \
 	std\uni.d std\base64.d std\ascii.d \
 	std\demangle.d std\uri.d std\metastrings.d std\mmfile.d std\getopt.d
@@ -171,7 +171,7 @@ SRC_STD= std\zlib.d std\zip.d std\stdint.d std\conv.d std\utf.d std\uri.d \
 	std\cstream.d std\demangle.d \
 	std\signals.d std\typetuple.d std\traits.d \
 	std\getopt.d \
-	std\variant.d std\numeric.d std\bitmanip.d std\complex.d std\mathspecial.d \
+	std\variant.d std\bitmanip.d std\complex.d std\mathspecial.d \
 	std\functional.d std\array.d std\typecons.d \
 	std\json.d std\xml.d std\encoding.d std\bigint.d std\concurrency.d \
 	std\stdiobase.d std\parallelism.d \
@@ -185,6 +185,8 @@ SRC_STD_REGEX= std\regex\internal\ir.d std\regex\package.d std\regex\internal\pa
 SRC_STD_RANGE= std\range\package.d std\range\primitives.d \
 	std\range\interfaces.d
 
+SRC_STD_NUMERIC= std\numeric\package.d std\numeric\summation.d
+
 SRC_STD_NET= std\net\isemail.d std\net\curl.d
 
 SRC_STD_LOGGER= std\experimental\logger\core.d std\experimental\logger\filelogger.d \
@@ -242,6 +244,7 @@ SRC_TO_COMPILE_NOT_STD= \
 SRC_TO_COMPILE= $(SRC_STD_ALL) \
 	$(SRC_STD_ALGO) \
 	$(SRC_STD_RANGE) \
+	$(SRC_STD_NUMERIC) \
 	$(SRC_TO_COMPILE_NOT_STD)
 
 SRC_ZLIB= \
@@ -342,7 +345,8 @@ DOCS=	$(DOC)\object.html \
 	$(DOC)\std_math.html \
 	$(DOC)\std_mathspecial.html \
 	$(DOC)\std_mmfile.html \
-	$(DOC)\std_numeric.html \
+	$(DOC)\std_numeric_package.html \
+	$(DOC)\std_numeric_summation.html \
 	$(DOC)\std_outbuffer.html \
 	$(DOC)\std_parallelism.html \
 	$(DOC)\std_path.html \
@@ -417,6 +421,7 @@ UNITTEST_OBJS= \
 unittest : $(LIB)
 	$(DMD) $(UDFLAGS) -L/co -c -unittest -ofunittest1.obj $(SRC_STD_1_HEAVY)
 	$(DMD) $(UDFLAGS) -L/co -c -unittest -ofunittest2.obj $(SRC_STD_RANGE)
+	$(DMD) $(UDFLAGS) -L/co -c -unittest -ofunittest2.obj $(SRC_STD_NUMERIC)
 	$(DMD) $(UDFLAGS) -L/co -c -unittest -ofunittest2a.obj $(SRC_STD_2a_HEAVY)
 	$(DMD) $(UDFLAGS) -L/co -c -unittest -ofunittest3.obj $(SRC_STD_3)
 	$(DMD) $(UDFLAGS) -L/co -c -unittest -ofunittest3a.obj $(SRC_STD_3a)
@@ -455,7 +460,8 @@ cov : $(SRC_TO_COMPILE) $(LIB)
 	$(DMD) -conf= -cov=93 -unittest -main -run std\csv.d
 	$(DMD) -conf= -cov=91 -unittest -main -run std\math.d
 	$(DMD) -conf= -cov=95 -unittest -main -run std\complex.d
-	$(DMD) -conf= -cov=70 -unittest -main -run std\numeric.d
+	$(DMD) -conf= -cov=70 -unittest -main -run std\numeric\package.d
+	$(DMD) -conf= -cov=70 -unittest -main -run std\numeric\summation.d
 	$(DMD) -conf= -cov=94 -unittest -main -run std\bigint.d
 	$(DMD) -conf= -cov=95 -unittest -main -run std\bitmanip.d
 	$(DMD) -conf= -cov=82 -unittest -main -run std\typecons.d
@@ -520,6 +526,7 @@ cov : $(SRC_TO_COMPILE) $(LIB)
 	$(DMD) -conf= -cov=92 -unittest -main -run std\internal\math\errorfunction.d
 	$(DMD) -conf= -cov=31 -unittest -main -run std\internal\windows\advapi32.d
 	$(DMD) -conf= -cov=58 -unittest -main -run etc\c\zlib.d
+>>>>>>> upstream/master
 
 html : $(DOCS)
 
@@ -705,8 +712,11 @@ $(DOC)\std_mathspecial.html : $(STDDOC) std\mathspecial.d
 $(DOC)\std_mmfile.html : $(STDDOC) std\mmfile.d
 	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_mmfile.html $(STDDOC) std\mmfile.d
 
-$(DOC)\std_numeric.html : $(STDDOC) std\numeric.d
-	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_numeric.html $(STDDOC) std\numeric.d
+$(DOC)\std_numeric.html : $(STDDOC) std\numeric\package.d
+	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_numeric_package.html $(STDDOC) std\numeric\package.d
+
+$(DOC)\std_numeric_summation.html : $(STDDOC) std\numeric\summation.d
+	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_numeric_summation.html $(STDDOC) std\numeric\summation.d
 
 $(DOC)\std_outbuffer.html : $(STDDOC) std\outbuffer.d
 	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_outbuffer.html $(STDDOC) std\outbuffer.d
@@ -882,7 +892,7 @@ zip : win32.mak win64.mak posix.mak $(STDDOC) $(SRC) \
 	$(SRC_ETC) $(SRC_ETC_C) $(SRC_ZLIB) $(SRC_STD_NET) $(SRC_STD_DIGEST) $(SRC_STD_CONTAINER) \
 	$(SRC_STD_INTERNAL) $(SRC_STD_INTERNAL_DIGEST) $(SRC_STD_INTERNAL_MATH) \
 	$(SRC_STD_INTERNAL_WINDOWS) $(SRC_STD_REGEX) $(SRC_STD_RANGE) $(SRC_STD_ALGO) \
-	$(SRC_STD_LOGGER)
+	$(SRC_STD_LOGGER) $(SRC_STD_NUMERIC) 
 	del phobos.zip
 	zip32 -u phobos win32.mak win64.mak posix.mak $(STDDOC)
 	zip32 -u phobos $(SRC)
@@ -905,6 +915,7 @@ zip : win32.mak win64.mak posix.mak $(STDDOC) $(SRC) \
 	zip32 -u phobos $(SRC_STD_CONTAINER)
 	zip32 -u phobos $(SRC_STD_REGEX)
 	zip32 -u phobos $(SRC_STD_RANGE)
+	zip32 -u phobos $(SRC_STD_NUMERIC)
 	zip32 -u phobos $(SRC_STD_ALGO)
 
 phobos.zip : zip
diff --git a/win64.mak b/win64.mak
index a4fe3451a13..3769b90ddf0 100644
--- a/win64.mak
+++ b/win64.mak
@@ -111,7 +111,7 @@ SRC_STD_1_HEAVY= std\stdio.d std\stdiobase.d \
 SRC_STD_2a_HEAVY= std\array.d std\functional.d std\path.d std\outbuffer.d std\utf.d
 
 SRC_STD_math=std\math.d
-SRC_STD_3= std\csv.d std\complex.d std\numeric.d std\bigint.d
+SRC_STD_3= std\csv.d std\complex.d std\bigint.d
 SRC_STD_3c= std\datetime.d std\bitmanip.d std\typecons.d
 
 SRC_STD_3a= std\uni.d std\base64.d std\ascii.d \
@@ -188,7 +188,7 @@ SRC_STD= std\zlib.d std\zip.d std\stdint.d std\conv.d std\utf.d std\uri.d \
 	std\cstream.d std\demangle.d \
 	std\signals.d std\typetuple.d std\traits.d \
 	std\getopt.d \
-	std\variant.d std\numeric.d std\bitmanip.d std\complex.d std\mathspecial.d \
+	std\variant.d std\bitmanip.d std\complex.d std\mathspecial.d \
 	std\functional.d std\array.d std\typecons.d \
 	std\json.d std\xml.d std\encoding.d std\bigint.d std\concurrency.d \
 	std\stdiobase.d std\parallelism.d \
@@ -202,6 +202,8 @@ SRC_STD_REGEX= std\regex\internal\ir.d std\regex\package.d std\regex\internal\pa
 SRC_STD_RANGE= std\range\package.d std\range\primitives.d \
 	std\range\interfaces.d
 
+SRC_STD_NUMERIC= std\numeric\package.d std\numeric\summation.d
+
 SRC_STD_NET= std\net\isemail.d std\net\curl.d
 
 SRC_STD_LOGGER= std\experimental\logger\core.d std\experimental\logger\filelogger.d \
@@ -259,6 +261,7 @@ SRC_TO_COMPILE_NOT_STD= \
 SRC_TO_COMPILE= $(SRC_STD_ALL) \
 	$(SRC_STD_ALGO) \
 	$(SRC_STD_RANGE) \
+	$(SRC_STD_NUMERIC) \
 	$(SRC_TO_COMPILE_NOT_STD)
 
 SRC_ZLIB= \
@@ -359,7 +362,8 @@ DOCS=	$(DOC)\object.html \
 	$(DOC)\std_math.html \
 	$(DOC)\std_mathspecial.html \
 	$(DOC)\std_mmfile.html \
-	$(DOC)\std_numeric.html \
+	$(DOC)\std_numeric_package.html \
+	$(DOC)\std_numeric_summation.html \
 	$(DOC)\std_outbuffer.html \
 	$(DOC)\std_parallelism.html \
 	$(DOC)\std_path.html \
@@ -444,6 +448,7 @@ UNITTEST_OBJS= \
 unittest : $(LIB)
 	$(DMD) $(UDFLAGS) -c -unittest -ofunittest1.obj $(SRC_STD_1_HEAVY)
 	$(DMD) $(UDFLAGS) -c -unittest -ofunittest2.obj $(SRC_STD_RANGE)
+	$(DMD) $(UDFLAGS) -c -unittest -ofunittest2.obj $(SRC_STD_NUMERIC)
 	$(DMD) $(UDFLAGS) -c -unittest -ofunittest2a.obj $(SRC_STD_2a_HEAVY)
 	$(DMD) $(UDFLAGS) -c -unittest -ofunittestM.obj $(SRC_STD_math)
 	$(DMD) $(UDFLAGS) -c -unittest -ofunittest3.obj $(SRC_STD_3)
@@ -663,8 +668,11 @@ $(DOC)\std_mathspecial.html : $(STDDOC) std\mathspecial.d
 $(DOC)\std_mmfile.html : $(STDDOC) std\mmfile.d
 	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_mmfile.html $(STDDOC) std\mmfile.d
 
-$(DOC)\std_numeric.html : $(STDDOC) std\numeric.d
-	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_numeric.html $(STDDOC) std\numeric.d
+$(DOC)\std_numeric.html : $(STDDOC) std\numeric\package.d
+	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_numeric_pacakge.html $(STDDOC) std\numeric\package.d
+
+$(DOC)\std_numeric_summation.html : $(STDDOC) std\numeric\summation.d
+	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_numeric_summation.html $(STDDOC) std\numeric\summation.d
 
 $(DOC)\std_outbuffer.html : $(STDDOC) std\outbuffer.d
 	$(DMD) -c -o- $(DDOCFLAGS) -Df$(DOC)\std_outbuffer.html $(STDDOC) std\outbuffer.d
@@ -840,7 +848,7 @@ zip : win32.mak win64.mak posix.mak $(STDDOC) $(SRC) \
 	$(SRC_ETC) $(SRC_ETC_C) $(SRC_ZLIB) $(SRC_STD_NET) $(SRC_STD_DIGEST) $(SRC_STD_CONTAINER) \
 	$(SRC_STD_INTERNAL) $(SRC_STD_INTERNAL_DIGEST) $(SRC_STD_INTERNAL_MATH) \
 	$(SRC_STD_INTERNAL_WINDOWS) $(SRC_STD_REGEX) $(SRC_STD_RANGE) $(SRC_STD_ALGO) \
-	$(SRC_STD_LOGGER)
+	$(SRC_STD_LOGGER) $(SRC_STD_NUMERIC)
 	del phobos.zip
 	zip32 -u phobos win32.mak win64.mak posix.mak $(STDDOC)
 	zip32 -u phobos $(SRC)
@@ -863,6 +871,7 @@ zip : win32.mak win64.mak posix.mak $(STDDOC) $(SRC) \
 	zip32 -u phobos $(SRC_STD_CONTAINER)
 	zip32 -u phobos $(SRC_STD_REGEX)
 	zip32 -u phobos $(SRC_STD_RANGE)
+	zip32 -u phobos $(SRC_STD_NUMERIC)
 	zip32 -u phobos $(SRC_STD_ALGO)
 
 phobos.zip : zip