I’d like to state in advance that I have negligible familiarity with number representations, so please excuse any inaccuracies. I also enjoy getting rid of legacy systems.
I have recently read an article which pointed out that one of the reasons COBOL is still somewhat prevalent in the banking and governmental sectors is that it used fixed point number representations, which are less prone to accumulating inaccuracies during the typical calculations these sectors use. (There are other factors of course as many languages let you use fixed point numbers with varying levels of difficulty by now.)
This article describes Muller’s recursion formula, which is a good example of this erratic behavior. This series should converge to 5, but using doubles it converges to a 100, after some hectic jumps. I wanted to try this with Julia’s types, so I wrote this code.
module cobol using FixedPointNumbers mullers(y,z) = 108 - (815 - 1500/z)/y function recur(t,n) v = zeros(t, n) v[1:2] .= t[4, 4.25] for ii = 3:n v[ii] = mullers(v[ii-1], v[ii-2]) end v end end # module
Float64 it reproduces the behavior of other languages, as expected. With
BigFloat it properly converges to 5. With
Float64.(recur(Rational,20)) it also converges to 5, but overflows before the double representation would reach it’s fixed point at iteration 29.
What I am interested in (apart from your insights in general) is that by using a fixed point representation
N41f23 it also erroneously converges to 100 (albeit differently than doubles - at least no sign change) - is there a way to do it right with fixed point representation that I haven’t found?