Order of Operation Matters Thanks to Rounding Error

My previous post was supposed to be followed up much quicker than this. Oh Well.

If that was The Good Stuff Is Hidden, then this is The Hidden Good Stuff Can Hurt You.

Joe Landman of “thanks for the awesome meeting space for MDLUG at SGI” fame had an excellent post on floating point arithmetic.

His post is what gave me the idea to do the same thing with Enumerable methods in the previous post.

I’d like to show that Joe’s assertion that programmers should worry about rounding error is true for .NET developers especially in this world of forthcoming parallel LINQ.

Lets write our own Range method which takes a start, end and increment parameters. You could use Enumerable.Range(start,count) and Enumerable.Reverse(this source). But reversing an Enumerable so big doesn’t sound like a good idea to me.

My brain things Range start & end, not start & count. I blame Python. This means the Range implementation is a little long because I try to detect infinite ranges.

public IEnumerable<int> Range(int start, int end, int increment)

{

    int i = start;

    if (start.CompareTo(end) <= 0 && increment > 0)

        while (i.CompareTo(end) <= 0)//(i <= end)

        {

            yield return i;

            i = i + increment;

        }

    else if (start.CompareTo(end) >= 0 && increment < 0)

        while (i.CompareTo(end) >= 0)//(i <= end)

        {

            yield return i;

            i = i + increment;

        }

    else

        throw new ArgumentOutOfRangeException();

}

Now lets use this method to count and sum the same exact same values twice. Once from start to end and again in a different order from end to start.

float x = 0f, y = 0f;

int N = 100000000;

 

x = Range(1, N, 1).Aggregate(x, (a, b) => a + 1 / (float)b);

y = Range(N, 1, -1).Aggregate(y, (a, b) => a + 1 / (float)b);

 

Console.WriteLine(“[increasing] sum = {0}”, x);

Console.WriteLine(“[decreasing] sum = {0}”, y);

This outputs the following:

[increasing] sum = 15.40368
[decreasing] sum = 18.80792

Not what you would expect is it? But before we go there…

I find the F# version of this to be easier and more readable.

let si = {1..100000000}

let f a b = a+1.0f/b

let i = Seq.fold f 0.0f (Seq.map (fun a-> (float32 a)) si)

let sd= StandardRanges.int 100000000 (-1) 1 //{100000000..1}

let d= Seq.fold f 0.0f (Seq.map (fun a-> (float32 a)) sd)

printf “[increasing] sum = %f\n” i

printf “[decreasing] sum = %f\n” d

This outputs the same :

[increasing] sum = 15.40368
[decreasing] sum = 18.80792

Now the question… why aren’t these numbers the same. It is the result of adding all of the same numbers.

Well if you read Joe Landman’s post that I linked above, then you already know. Order of operations matters! Floating point rounding error accumulates and it adds up big here. See Joe’s post for more.

When I read this, the thought that PLINQ might be easy quickly left my mind. Results changing just because the order is different!

Now what is worse, lets say you try to do the same thing in F#, but rather than using integers and casting like we did in C#, we just use F#’s StandardRanges library to generate floating point values for us.

let n = 100000000.0f

printf “[increasing] sum = %f\n” (Seq.fold (fun a b -> a+1.0f/b) 0.0f (StandardRanges.generate 0.0f (+) 1.0f 1.0f n))

printf “[decreasing] sum = %f\n” (Seq.fold (fun a b -> a+1.0f/b) 0.0f (StandardRanges.float32 n -1.0f 1.0f))

The problem here is the way the floating point numbers are represented is slightly different than the above integer cast version. The result is  disappointingly different:

[increasing] sum = 15.40368
[decreasing] sum = Infinity

Infinity? Maybe I’ll have to look into this some more.

Some interesting things to note:

• StandardRanges had some issues in 1.9.2.9 Thanks to Dustin Cambell for trying the same code which ran fine on his 1.9.3.7 and making me upgrade.
• F#’s Seq is just an alias for IEnuemerable<> but it is augmented with more things.
• examining the source in prim-types.fs was awesome. thanks for the source Microsoft.
• Looks like prim-types.fs and StandardRanges has changed a bit from 1.9.2.9 to 1.9.3.7 and there is some float specialty stuff in there now. I need to investigate.
• float in F# is double in C#. float32 in F# is float in C#. or if you remember you SATs : C#float:F#float32::C#double:F#float

Wondering about using a 64bit type instead of a 32bit type? Scroll up and click my link Joe Landman’s post. It doesn’t make this issue go away. We need to actually think about what we are doing as programmers *gasp*.