# Numerically stable sum of doubles

It's known that if you naively sum up a collection of arbitrary floating-point numbers using any floating-point standard (e.g. floats or doubles with the C standard), the end result may differ, sometimes significantly, from the actual sum. This is due to the loss of precision involved in approximating a real number in a finite number of bits (e.g. 64). In some cases the rounding is catastrophic -- particularly if you have any cases of (large positive number + large negative number) or (large number + small number) as this causes a phenomenon known as "cancellation" -- basically, because precision is finite, there will always be cases where not enough precision is available.

Sometimes this numerical instability is acceptable. Sometimes it isn't. For those cases, I decided to write a class (stealing the "stable summation" code from Python):

EDIT: WARNING THIS CODE HAS A BUG IN IT, DON'T USE IT

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Snowbody
{
class StableSummer
{
List<double> partials;
public double Sum()
{
return partials.Sum();
}

public StableSummer(StableSummer m)
{
partials = new List<double>();
}

public StableSummer()
{
partials = new List<double>();
}

public void AddRange(IList<double> x)
{
foreach (var w in x)
{
}
}

public void Add(double xp)
{
// translation of http://code.activestate.com/recipes/393090/
// In the event Licensee prepares a derivative work that is based on or incorporates Python 2.7.7 or any part thereof, and wants to make the derivative work available to others as provided herein, then Licensee hereby agrees to include in any such work a brief summary of the changes made to Python 2.7.7.
// brief summary: translated from python to C#
// "Full precision summation using multiple floats for intermediate values"
// Rounded x+y stored in hi with the round-off stored in lo.  Together
// hi+lo are exactly equal to x+y.  The inner loop applies hi/lo summation
// to each partial so that the list of partial sums remains exact.
// Depends on IEEE-754 arithmetic guarantees.  See proof of correctness at:
// www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
int i = 0;
double x = xp;
foreach (var yp in partials)
{
double y;
if (Math.Abs(xp) < Math.Abs(yp))
{
x = yp;
y = xp;
}
else
{
y = yp;
}
double hi = x + y;
double lo = y - (hi - x);
if (lo != 0.0)
{
if (i == partials.Count)
{
}
else
{
partials[i] = lo;
}
i++;
}
x = hi;
}
partials.RemoveRange(i, partials.Count - i);
}
}
}


Usage: After instantiation, repeatedly call .Add() and/or .AddRange() and then call .Sum() when you want the (closest possible double to the) true sum.

I believe the code is giving the right answer. I found a bug while I was refactoring; it's supposed to reassign to xp at the end of the loop body.

My questions are: - Any refactoring necessary? - Should I make the class implement IList<double> or like an IEnumerable<double>? - Would it make sense to include a cast to double that implements the sum?

This code has a lot of one letter parameters. I would try to expand them a bit for readability. Of course, that's easier said than done when we're talking about primarily mathematical functions.

• (1) Thanks, and as I was doing this I found I had made a mistake in converting the code! (2) which lines have hard-to-read indentation? Which lines are inconsistent? – Snowbody Jul 9 '14 at 14:19
• My bad. My mobile broswer must have interpretted it wrong. There might be a mix of tabs & spaces in there. I'll pull that out of my answer. – RubberDuck Jul 9 '14 at 14:21
• I didn't think I had any tabs in there (my editor is set to convert tabs to spaces), I'll check. – Snowbody Jul 9 '14 at 14:22

Two small notes:

1. You should use constructor chaining for maintenance reason.
2. You should integrate a null check inside your overloaded constructor

public StableSummer(StableSummer m):this()
{
if (m != null)
{
}
}

public StableSummer()
{
partials = new List<double>();
}


Otherwise it looks good for me, except i wouldnt use an anonymous var in the foreach loop and i would name the vars in the Add() method more meaningful, but this is just a matter of taste.

• 1. thanks, I always forget about that. 2. Why is this necessary? why is "do nothing" better than a throw? – Snowbody Jul 9 '14 at 14:18
• For 2. It would just work without matter if the passed StableSummer is null or not – Heslacher Jul 9 '14 at 14:20
• I don't understand. If the passed StableSummer is null, why shouldn't that be an error? Why do you consider "working" to be the same as "doing nothing"? – Snowbody Jul 9 '14 at 14:22
• First for me "working" == "No Exception" and my preference: An avoided exception is the best exception. – Heslacher Jul 9 '14 at 14:30

You had a statement in your question that goes like this

Usage: After instantiation, repeatedly call .Add() and/or .AddRange() and then call .Sum() when you want the (closest possible double to the) true sum.

You had to write this because Add and Sum are synonyms. You could perhaps be better fitted if you change Add to something like AddAddend. Or you could only support one method with variable parameters like:

public void AddAddends(params double[] values){
foreach (var value in values)
{
//this method would be private

List<double> list = new List<double>(){1.4, 1.9, 2.5};

• The original reason I used .Add(), .AddRange(), and .Sum() was because this code was replacing a List<double>` which implements precisely those methods with the same functionality. – Snowbody Jul 14 '14 at 15:11