raku:34 python:19 extreme math

Coming off the excellent raku weekly news, my curiosity was piqued by a tweet about big-endian smells that referenced a blog about “extreme math”. After getting my fill of COBOL mainframe nostalgia, the example of Muller’s Recurrence got me thinking.

The simple claim made in the tweet thread was:

Near the end it [the blog] states that no modern language has fixed point, but Raku (formerly Perl6) has a built in rational type which is quite an interesting comparison. It keeps two integers for the numerator and the denominator and no loss of precision occurs.

I have also covered some of the benefits of the raku approach to math in a previous blog Machine Math and Raku, often the example given is 0.1 + 0.2 =>0.3 which trips up a lot of languages. I like this example, but I am not entirely convinced by it – sure it can be odd when a programming newbie sees a slightly different result caused by floating point conversions – but it is too mickey mouse to be a serious concern.

The Muller Extreme Challenge

This challenge starts with seemingly innocuous equations and quickly descends into very substantial errors. To quote from the Technical Archaelogist blog:

Jean-Michel Muller is a French computer scientist with perhaps the best computer science job in the world. He finds ways to break computers using math. I’m sure he would say he studies reliability and accuracy problems, but no no no: He designs math problems that break computers. One such problem is his recurrence formula. Which looks something like this:

That doesn’t look so scary does it? The recurrence problem is useful for our purposes because:

  • It is straight forward math, no complicated formulas or concepts
  • We start off with two decimal places, so it’s easy to imagine this happening with a currency calculation.
  • The error produced is not a slight rounding error but orders of magnitude off.

And here’s a quick python script that produces floating point and fixed point versions of Muller’s Recurrence side by side:

from decimal import Decimal

def rec(y, z):
 return 108 - ((815-1500/z)/y)
 
def floatpt(N):
 x = [4, 4.25]
 for i in range(2, N+1):
  x.append(rec(x[i-1], x[i-2]))
 return x
 
def fixedpt(N):
 x = [Decimal(4), Decimal(17)/Decimal(4)]
 for i in range(2, N+1):
  x.append(rec(x[i-1], x[i-2]))
 return x

N = 30
flt = floatpt(N)
fxd = fixedpt(N)

for i in range(N):
 print( str(i) + ' | '+str(flt[i])+' | '+str(fxd[I]) )

Which gives us the following output:

i  | floating pt    | fixed pt
-- | -------------- | ---------------------------
0  | 4              | 4
1  | 4.25           | 4.25
2  | 4.47058823529  | 4.4705882352941176470588235
3  | 4.64473684211  | 4.6447368421052631578947362
4  | 4.77053824363  | 4.7705382436260623229461618
5  | 4.85570071257  | 4.8557007125890736342039857
6  | 4.91084749866  | 4.9108474990827932004342938
7  | 4.94553739553  | 4.9455374041239167246519529
8  | 4.96696240804  | 4.9669625817627005962571288
9  | 4.98004220429  | 4.9800457013556311118526582
10 | 4.9879092328   | 4.9879794484783912679439415
11 | 4.99136264131  | 4.9927702880620482067468253
12 | 4.96745509555  | 4.9956558915062356478184985
13 | 4.42969049831  | 4.9973912683733697540253088
14 | -7.81723657846 | 4.9984339437852482376781601
15 | 168.939167671  | 4.9990600687785413938424188
16 | 102.039963152  | 4.9994358732880376990501184
17 | 100.099947516  | 4.9996602467866575821700634
18 | 100.004992041  | 4.9997713526716167817979714
19 | 100.000249579  | 4.9993671517118171375788238
20 | 100.00001247862016 | 4.9897059157620938291040004
21 | 100.00000062392161 | 4.7951151851630947311130380
22 | 100.0000000311958  | 0.7281074924258006736651754
23 | 100.00000000155978 | -581.7081261405031229400219627
24 | 100.00000000007799 | 105.8595186892360167901632650
25 | 100.0000000000039  | 100.2767586430669099906187869
26 | 100.0000000000002  | 100.0137997241561168045699158
27 | 100.00000000000001 | 100.0006898905241097140861868
28 | 100.0 | 100.0000344942738135445216746
29 | 100.0 | 100.0000017247126631766583580
30 | 100.0 | 100.0000000862356186943169827

Up until about the 12th iteration the rounding error seems more or less negligible but things quickly go off the rails. Floating point math converges around a number twenty times the value of what the same calculation with fixed point math produces.

Least you think it is unlikely that anyone would do a recursive calculation so many times over. This is exactly what happened in 1991 when the Patriot Missile control system miscalculated the time and killed 28 people. And it turns out floating point math has blown lots of stuff up completely by accident. Mark Stadtherr gave an incredible talk about this called High Performance Computing: are we just getting wrong answers faster? You should read it if you want more examples and a more detailed history of the issue than I can offer here.

[endquote]

So, basically, python Float dies at iteration #12 and python Fixed/Decimal dies at iteration #19. According to the source text COBOL dies at iteration #18. Then the argument focuses on the need for the Decimal library.

How does raku Measure Up?

I do not buy the no loss of precision occurs claim made on twitter beyond the simpler examples, but I do think that Rats should fare well in the face of this kind of challenge. Here’s my code with raku default math:

my \N = 30;
my \x = []; 
x[0] = 4; 
x[1] = 4.25;

sub f(\y,\z) { 
    108 - ( (815 - 1500/z ) / y ) }

for 2..N -> \i { 
    x[i] = f(x[i-1],x[i-2])   }

for 0..N -> \i {
    say( i ~ ' | ' ~ x[i] )   }

Quick impression is that raku is a little more faithful to the mathematical description and a little less cramped than the python.

The raku output gives:

0 | 4
1 | 4.25
2 | 4.470588
3 | 4.644737
4 | 4.770538
5 | 4.855701
6 | 4.910847
7 | 4.945537
8 | 4.9669626
9 | 4.9800457
10 | 4.98797945
11 | 4.992770288
12 | 4.9956558915
13 | 4.9973912684
14 | 4.99843394394
15 | 4.999060071971
16 | 4.999435937147
17 | 4.9996615241038
18 | 4.99979690071342
19 | 4.99987813547793
20 | 4.9999268795046
21 | 4.9999561270611577
22 | 4.99997367600571244
23 | 4.99998420552027271
24 | 4.999990523282227659
25 | 4.9999943139585595936
26 | 4.9999965883712560237
27 | 4.99999795302135690799
28 | 4.999998771812315
29 | 4.99999926308729
30 | 4.999999557853926

So, 30 iterations with no loss of precision – and with the native raku math defaults. Nice!

Eventually raku breaks at 34 iterations, so raku:34, python:19.

~p6steve

PS. And to reflect the harsh reality of life, Victor Ejikhout’s comment can have the final word: so know your own limits!

This is not a problem of fixed point vs floating point. I think your examples favor Fix because you give it more digit of accuracy. What would happen if you used a Float format where the mantissa is equally long as the total Fix length? Objection #2: I think Cobol / Fix would converge away from 5 if you ran more iterations. The Muller equation has three fixed points: x_n==3, x_n==5, and x_n==100. If you start close enough to 5 it will converge there for a while, but (I’m guessing here; didn’t run all the tests) it will converge to the 100 solution. Since you give the float solution less precision it simply converges there faster.The only real lesson here is not to code unstable recursions.

7 thoughts on “raku:34 python:19 extreme math”

  1. You can actually see when and why we start to lose precision (after 28 iterations):

    > dd $_ for (4, 4.25, 108 – (815 – 1500 / * ) / * … *)[^30].kv
    0|4
    1|4.25
    2
    […]
    27
    28|4.999998771812315e0
    29|4.99999926308729e0

    Liked by 1 person

  2. Apparently the smaller and greater than characters in what I posted above are tripping up the blog software. Apologies.
    I removed them.

    dd $_ for (4, 4.25, 108 – (815 – 1500 / * ) / * … *)[^30].kv
    0
    4
    1
    4.25
    2
    76/17
    3
    353/76
    4
    1684/353
    5
    8177/1684
    6
    40156/8177
    7
    198593/40156
    8
    986404/198593
    9
    4912337/986404
    10
    24502636/4912337
    11
    122336033/24502636
    12
    611148724/122336033
    13
    3054149297/611148724
    14
    15265963516/3054149297
    15
    76315468673/15265963516
    16
    381534296644/76315468673
    17
    1907542343057/381534296644
    18
    9537324294796/1907542343057
    19
    47685459212513/9537324294796
    20
    238423809278164/47685459212513
    21
    1192108586037617/238423809278164
    22
    5960511549128476/1192108586037617
    23
    29802463602463553/5960511549128476
    24
    149012035582781284/29802463602463553
    25
    745059330625296977/149012035582781284
    26
    3725294111260656556/745059330625296977
    27
    18626462930705797793/3725294111260656556
    28
    4.999998771812315e0
    29
    4.99999926308729e0

    Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s