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.

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

LikeLiked by 1 person

hi Markus – I see you refer to the ‘e’ reflecting that the Num is now FP – also the sequence code version – super tight!

LikeLike

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

LikeLike

When we explicily specify FatRat, it even works up to 5000 (and maybe more)

my $N = 5000;

my FatRat @x = [];

@x[0] = FatRat(4);

@x[1] = FatRat(4.25);

LikeLike

Hi Jörn – that’s a good point … my article specifically focuses on Rat as an improvement over the default Float used in many languages and the benefits in terms of users getting the answers they expect 0.1+0.2=0.3 and a more natural number type that balances performance against this kind of algorithmic degradation. As you mention, the built in FatRat does much better … ! (IMO the language designers steered the right path by making Rat the default and letting us reach for the FatRat wrench when we need it.)

LikeLike

I don’t like the way is is done in the moment. We don’t use Raku because it is lightning fast, but because it works. So the logical way would be to automatically turn a Rat into a FatRat instead of a Float. Or better only use FatRats per default.

LikeLike

I guess we will have to agree to disagree on this one. Having spent large chunks of time writing test vectors for FPU silicon, IMO it would be a shame for to just to ignore all that performant silicon real estate and to know that whatever jnthn does raku can never match other VM language speeds. 😉

LikeLike