A friend of mine recently proposed the following problem on twitter:
Given f y z =108 – (815 – 1500/z)/y, x0 = 4, x1 = 4.25 and xN+1 = f xN xN-1 then what is x80? Note: This might not be as easy as it looks
— Mårten Rånge (@marten_range) October 19, 2013
I didn’t pay much attention to his warning and fell for the trap. I though I could write a small program in two minutes to compute the series and find what was the value of x[80]. So here is (a slightly modified version of) the C++ code that I put together in a couple of minutes.
C++
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
#include <iostream> #include <iomanip> template <typename T, int Size = 80> T x(int n) { static T cache[Size + 1] = {0}; if(n == 0) cache[n] = 4.0; else if(n == 1) cache[n] = 4.25; else cache[n] = 108 - (815 - 1500 / cache[n - 2]) / cache[n - 1]; return cache[n]; } int main() { for(int i = 0; i <= 80; ++i) { std::cout << "x[" << i << "]=" << std::setprecision(15) << x<double>(i) << std::endl; } return 0; } |
When I ran it I was surprised to notice that the series was converging to 100 by x[26].
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
x[0]=4 x[1]=4.25 x[2]=4.47058823529412 x[3]=4.64473684210522 x[4]=4.77053824362508 x[5]=4.85570071256856 x[6]=4.91084749866063 x[7]=4.94553739553051 x[8]=4.966962408041 x[9]=4.98004220429301 x[10]=4.98790923279579 x[11]=4.99136264131455 x[12]=4.96745509555227 x[13]=4.42969049830883 x[14]=-7.81723657845932 x[15]=168.939167671065 x[16]=102.039963152059 x[17]=100.09994751625 x[18]=100.004992040972 x[19]=100.000249579237 x[20]=100.00001247862 x[21]=100.000000623922 x[22]=100.000000031196 x[23]=100.00000000156 x[24]=100.000000000078 x[25]=100.000000000004 x[26]=100 ... x[80]=100 |
Actually, the initial program didn’t call std::setprecision and the numbers you get without that are less precise, but that doesn’t change the convergence, since it is just a printing artifact.
|
1 2 3 4 5 6 7 |
x[0]=4 x[1]=4.25 x[2]=4.47059 x[3]=4.64474 x[4]=4.77054 x[5]=4.8557 ... |
Finding the series interesting I searched a bit and then I understood his warning. I found this was a well known problem proposed around 1980 by Jean-Michel Muller and discussed in several papers by Prof. W. Kahan.
Given the function
Œ(y, z) := 108 – ( 815 – 1500/z )/y
and initial values x0 := 4 and x1 := 4.25 , define xn+1 := Œ(xn, xn-1) for n = 1, 2, 3, … in turn.
Our task is to compute xN for some moderately big preassigned integer N, say N = 80 .
For details see How Futile are Mindless Assessments of Roundoff in Floating-Point Computation? and Three Problems for Math.
This exercise is intended to show the problem that arises in using floating-point numbers. The float and double (both an IEEE Standard for Floating-Point Arithmetic, IEEE 754) representations use inverse powers of 2, which means most numbers require a an infinite number of bits for a precise representation. Numbers such as 0.25 or 0.875 can be exactly encoded as 1/4 and 1/2+1/4+1/8, but numbers such as 0.10 cannot be encoded with a finite sum of such terms. The result is problems with accuracy of calculations. Rand-offs can propagate through calculations in unexpected ways, just like Muller’s recurrence shows.
The actual limit of Muller’s series is not 100, but 5.
I was curios then how the decimal type from .NET compares to double. decimal (that uses base 10 instead of 2) has more precision (but a smaller range) than float or double which makes is more suitable for some applications, such as financial ones. (For a discussion of when to use decimal and when to use double see decimal vs double! – Which one should I use and when?).
So here is my C# program that uses decimal.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
class MullersRecurrence { static decimal[] cache = new decimal[100]; public decimal x(int n) { if (n == 0) cache[n] = 4m; else if (n == 1) cache[n] = 4.25m; else cache[n] = 108 - (815 - 1500 / cache[n - 2]) / cache[n - 1]; return cache[n]; } } class Program { static void Main(string[] args) { var mr = new MullersRecurrence(); for(int i = 0; i <= 80; ++i) { Console.WriteLine("x[{0}]={1}", i, mr.x(i)); } } } |
The output of this program is:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
x[0]=4 x[1]=4.25 x[2]=4.47058823529411764705882353 x[3]=4.64473684210526315789473686 x[4]=4.77053824362606232294617603 x[5]=4.85570071258907363420428376 x[6]=4.91084749908279320044042857 x[7]=4.94553740412391672477683015 x[8]=4.96696258176270059878160878 x[9]=4.98004570135563116267108889 x[10]=4.98797944847839228829979003 x[11]=4.99277028806206866201151005 x[12]=4.99565589150664533306792637 x[13]=4.99739126838157043427422171 x[14]=4.99843394394934565979621707 x[15]=4.99906007206149646425952424 x[16]=4.99943593895922460992955065 x[17]=4.99966156035548033890851805 x[18]=4.99979762579572007199519838 x[19]=4.99989263769854913604459541 x[20]=5.00021692999623515255759378 x[21]=5.00575688343630115907717069 x[22]=5.11585535860978057839952266 x[23]=7.26513170553842597520695497 x[24]=36.178328937337879304087182981 x[25]=91.17958879988455033108590199 x[26]=99.51631713443793014723080822 x[27]=99.97569833055963020623148188 x[28]=99.99878462167868201734350518 x[29]=99.99993923036059445960870932 x[30]=99.99999696151664049461733529 x[31]=99.99999984807584112595945239 x[32]=99.99999999240379245628007687 x[33]=99.99999999962018963513083004 x[34]=99.99999999998100948212683970 x[35]=99.99999999999905047411745292 x[36]=99.99999999999995252370620598 x[37]=99.99999999999999762618532030 x[38]=99.99999999999999988130926632 x[39]=99.99999999999999999406546333 x[40]=99.99999999999999999970327317 x[41]=99.99999999999999999998516366 x[42]=99.99999999999999999999925818 x[43]=99.99999999999999999999996291 x[44]=99.99999999999999999999999815 x[45]=99.99999999999999999999999991 x[46]=100.00000000000000000000000000 x[47]=100 ... x[49]=100 |
This represents and improvement, but in the end suffers from the same accumulated round-offs problem. It takes more iterations, but eventually the series also converges to 100.
My friend then suggested trying a data type that doesn’t suffer from rounding issues. Such a type is BigRational for F# (it can be used with any .NET language). It is available in the F# PowerPack that is an open-source project available on CodePlex. Below is the F# equivalent of the previous program that uses BigRational.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
open Microsoft.FSharp.Math;; let cache = Array.create 100 BigRational.Zero let x n = match n with | 0 -> cache.[n] <- 4N | 1 -> cache.[n] <- 17N/4N | _ -> cache.[n] <- 108N - (815N - 1500N / cache.[n - 2]) / cache.[n - 1] cache.[n] [<EntryPoint>] let main argv = for i in 0 .. 80 do System.Console.WriteLine(double(x i)) 0 |
The output looks like this:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 |
x[0]=4 x[1]=4.25 x[2]=4.47058823529412 x[3]=4.64473684210526 x[4]=4.77053824362606 x[5]=4.85570071258907 x[6]=4.91084749908279 x[7]=4.94553740412392 x[8]=4.9669625817627 x[9]=4.98004570135563 x[10]=4.98797944847839 x[11]=4.99277028806207 x[12]=4.99565589150663 x[13]=4.99739126838134 x[14]=4.99843394394482 x[15]=4.99906007197089 x[16]=4.99943593714684 x[17]=4.99966152410377 x[18]=4.99979690071342 x[19]=4.99987813547793 x[20]=4.9999268795046 x[21]=4.99995612706116 x[22]=4.99997367600571 x[23]=4.99998420552027 x[24]=4.99999052328223 x[25]=4.99999431395856 x[26]=4.99999658837126 x[27]=4.99999795302136 x[28]=4.99999877181231 x[29]=4.99999926308721 x[30]=4.99999955785226 x[31]=4.99999973471133 x[32]=4.99999984082679 x[33]=4.99999990449607 x[34]=4.99999994269764 x[35]=4.99999996561859 x[36]=4.99999997937115 x[37]=4.99999998762269 x[38]=4.99999999257362 x[39]=4.99999999554417 x[40]=4.9999999973265 x[41]=4.9999999983959 x[42]=4.99999999903754 x[43]=4.99999999942252 x[44]=4.99999999965351 x[45]=4.99999999979211 x[46]=4.99999999987527 x[47]=4.99999999992516 x[48]=4.9999999999551 x[49]=4.99999999997306 x[50]=4.99999999998384 x[51]=4.9999999999903 x[52]=4.99999999999418 x[53]=4.99999999999651 x[54]=4.9999999999979 x[55]=4.99999999999874 x[56]=4.99999999999925 x[57]=4.99999999999955 x[58]=4.99999999999973 x[59]=4.99999999999984 x[60]=4.9999999999999 x[61]=4.99999999999994 x[62]=4.99999999999996 x[63]=4.99999999999998 x[64]=4.99999999999999 x[65]=4.99999999999999 x[66]=5 x[67]=5 x[68]=5 x[69]=5 x[70]=5 x[71]=5 x[72]=5 x[73]=5 x[74]=5 x[75]=5 x[76]=5 x[77]=5 x[78]=5 x[79]=5 x[80]=5 |
Now this is a totally different story. The values converge to the expected value of 5.
You probably noticed the casting to double for printing. It is necessary because otherwise the output would look like this:
|
1 2 3 4 5 6 7 |
x[0]=4 x[1]=17/4 x[2]=76/17 x[3]=353/76 ... x[79]=41359030627651383817474849310671104336332210648235594113/8271806125530276773348891823090615755005322810072671996 x[80]=206795153138256918939565417139009598365577843034794672964/41359030627651383817474849310671104336332210648235594113 |
That isn’t very helpful. I can’t even read these insane large numbers, not to mention dividing them. So to get the actual real number and be able to compare with the previous programs a conversion to double is necessary.
As I mentioned earlier, BigRational can also be used from C#.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
class MullersRecurrenceBigRational { BigRational [] cache = new BigRational[100]; public BigRational x(int n) { if (n == 0) cache[n] = BigRational.FromInt(4); else if (n == 1) cache[n] = BigRational.FromInt(17)/BigRational.FromInt(4); else cache[n] = BigRational.FromInt(108) - (BigRational.FromInt(815) - (BigRational.FromInt(1500) / cache[n - 2])) / cache[n - 1]; return cache[n]; } } class Program { static void Main(string[] args) { var mr = new MullersRecurrenceBigRational(); for(int i = 0; i <= 80; ++i) { Console.WriteLine("x[{0}] = {1}", i, (double)mr.x(i)); } } } |
It yields the very same output as the F# program so I will not list it again. However, below is a comparison table with the results for various number data types.
| index | C++ with float | C++/C# with double | C# with decimal | C#/F# with BigRational |
|---|---|---|---|---|
| 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 |
4 4.25 4.47058868408203 4.64474487304688 4.77070617675781 4.85921478271484 4.98312377929688 6.39543151855469 27.6326293945313 86.9937591552734 99.2555084228516 99.9625854492188 99.9981307983398 99.9999084472656 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 |
4 4.25 4.47058823529412 4.64473684210522 4.77053824362508 4.85570071256856 4.91084749866063 4.94553739553051 4.966962408041 4.98004220429301 4.98790923279579 4.99136264131455 4.96745509555227 4.42969049830883 -7.81723657845932 168.939167671065 102.039963152059 100.09994751625 100.004992040972 100.000249579237 100.00001247862 100.000000623922 100.000000031196 100.00000000156 100.000000000078 100.000000000004 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 |
4 4.25 4.47058823529411764705882353 4.64473684210526315789473686 4.77053824362606232294617603 4.85570071258907363420428376 4.91084749908279320044042857 4.94553740412391672477683015 4.96696258176270059878160878 4.98004570135563116267108889 4.98797944847839228829979003 4.99277028806206866201151005 4.99565589150664533306792637 4.99739126838157043427422171 4.99843394394934565979621707 4.99906007206149646425952424 4.99943593895922460992955065 4.99966156035548033890851805 4.99979762579572007199519838 4.99989263769854913604459541 5.00021692999623515255759378 5.00575688343630115907717069 5.11585535860978057839952266 7.26513170553842597520695497 36.178328937337879304087182981 91.17958879988455033108590199 99.51631713443793014723080822 99.97569833055963020623148188 99.99878462167868201734350518 99.99993923036059445960870932 99.99999696151664049461733529 99.99999984807584112595945239 99.99999999240379245628007687 99.99999999962018963513083004 99.99999999998100948212683970 99.99999999999905047411745292 99.99999999999995252370620598 99.99999999999999762618532030 99.99999999999999988130926632 99.99999999999999999406546333 99.99999999999999999970327317 99.99999999999999999998516366 99.99999999999999999999925818 99.99999999999999999999996291 99.99999999999999999999999815 99.99999999999999999999999991 100.00000000000000000000000000 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 |
4 4.25 4.47058823529412 4.64473684210526 4.77053824362606 4.85570071258907 4.91084749908279 4.94553740412392 4.9669625817627 4.98004570135563 4.98797944847839 4.99277028806207 4.99565589150663 4.99739126838134 4.99843394394482 4.99906007197089 4.99943593714684 4.99966152410377 4.99979690071342 4.99987813547793 4.9999268795046 4.99995612706116 4.99997367600571 4.99998420552027 4.99999052328223 4.99999431395856 4.99999658837126 4.99999795302136 4.99999877181231 4.99999926308721 4.99999955785226 4.99999973471133 4.99999984082679 4.99999990449607 4.99999994269764 4.99999996561859 4.99999997937115 4.99999998762269 4.99999999257362 4.99999999554417 4.9999999973265 4.9999999983959 4.99999999903754 4.99999999942252 4.99999999965351 4.99999999979211 4.99999999987527 4.99999999992516 4.9999999999551 4.99999999997306 4.99999999998384 4.9999999999903 4.99999999999418 4.99999999999651 4.9999999999979 4.99999999999874 4.99999999999925 4.99999999999955 4.99999999999973 4.99999999999984 4.9999999999999 4.99999999999994 4.99999999999996 4.99999999999998 4.99999999999999 4.99999999999999 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 |
The conclusion is that you should be aware that round-offs can accumulate and lead to unexpected results. Use the most appropriate data types suitable. Do not use double (not to mention float) for financial data.
Example in boost::multiprecision
#include “stdafx.h”
#include
#include
#include
using boost::multiprecision::cpp_dec_float_100;
template
T x(int n)
{
static T cache[Size + 1] = { 0 };
if (n == 0) cache[n] = 4.0;
else if (n == 1) cache[n] = 4.25;
else cache[n] = 108 – (815 – 1500 / cache[n – 2]) / cache[n – 1];
return cache[n];
}
int _tmain(int argc, _TCHAR* argv[])
{
for (int i = 0; i <= 80; ++i)
{
std::cout << "x[" << i << "]=" << std::setprecision(15) << x(i) << std::endl;
}
return 0;
}
Very nice problem, very nice solution, and very interesting article (also the original). Left me wondered how short/fast would it be in Python:
from fractions import Fraction as F
def fn(y,z):
return F(108) – (F(815) – F(1500)/z)/y
def run(N):
a, b = F(4), F(4.25)
for n in range(N):
print n, a
b, a = fn(b,a), b
return b
On my machine Python executes
run(1000)in 1.6 seconds (without the print) and `run(2000) in 10 seconds.