# Losing precision

Here is the Ruby code.

``````ar = [0, 2, -4]
100.times{ar << (111 - 1130.quo(ar[-1])+ 3000.quo(ar[-1]*ar[-2])) }

[3, 4, 5, 6, 7, 8, 20, 30, 50, 100].each do |n|
puts "%3d -> %0.16f" % [n, ar[n]]
end
``````

And its results.

``````➜ ruby mullersequence.rb
3 -> 18.5000000000000000
4 -> 9.3783783783783784
5 -> 7.8011527377521614
6 -> 7.1544144809752494
7 -> 6.8067847369236330
8 -> 6.5926327687044384
20 -> 6.0435521101892689
30 -> 6.0067860930312058
50 -> 6.0001758466271872
100 -> 6.0000000193194779
``````

Here’s my closest Crystal (1.2.2) translation.

``````require "big"

ar = [0.to_big_d, 2.to_big_d, -4.to_big_d]

100.times{ar << (111 - (1130 / (ar[-1])).to_big_d + (3000 / (ar[-1] * ar[-2]))).to_big_d }

[3, 4, 5, 6, 7, 8, 20, 30, 50, 100].each do |n|
puts "%3d -> %0.16f" % [n, ar[n]]
end
``````

And it’s results.

``````➜ crystal mullersequence.cr
3 -> 18.5000000000000000
4 -> 9.3783783783783790
5 -> 7.8011527377521617
6 -> 7.1544144809752490
7 -> 6.8067847369236327
8 -> 6.5926327687044388
20 -> 6.0435521101892693
30 -> 6.0067860930312058
50 -> 6.0001758466271875
100 -> 99.9999999999992468
``````

It’s losing precision on the last iteration.
How do I fix this?

Increasing the number of iterations for division improves matters.

``````require "big"

ar = [0.to_big_d, 2.to_big_d, -4.to_big_d]

100.times{
t1 = 1130.to_big_d.div(ar[-1], 200)
t2 = 3000.to_big_d.div(ar[-1] * ar[-2], 200)
ar << 111.to_big_d - t1 + t2
}

[3, 4, 5, 6, 7, 8, 20, 30, 50, 100].each do |n|
puts "%3d -> %0.16f" % [n, ar[n]]
end

``````
``````crystal run mullersequence.cr
3 -> 18.5000000000000000
4 -> 9.3783783783783790
5 -> 7.8011527377521617
6 -> 7.1544144809752490
7 -> 6.8067847369236327
8 -> 6.5926327687044388
20 -> 6.0435521101892693
30 -> 6.0067860930312058
50 -> 6.0001758466271875
100 -> 6.0000000193194776
``````

I checked with Optimize `BigDecimal#div` for inexact divisions by HertzDevil · Pull Request #10803 · crystal-lang/crystal · GitHub, it still produces the same output.

Just increasing the precision improves the results by a lot. I found the 100ths value starts becoming reasonable at around 120.

Maybe we should increase the default? I don’t know where the number (`100`) is coming from.

This is the first exercise that arises as part of a Rosetta Code problem.

https://rosettacode.org/wiki/Pathological_floating_point_problems

It’s a form of a Muller Recurrence problem, which will converge towards different results depending upon a systems/language floating point accuracy.

https://scipython.com/blog/mullers-recurrence/

Edit:
Currently: DEFAULT_MAX_DIV_ITERATIONS = `100_u64`

Changing `.div(x, 200)` to `.div(x, 132)` produced same results.
It’s interesting how different smaller values affect the results so widely.

1 Like

For the 2nd exercise:
Ruby sets the value of E to 50 decimal places: `BigMath.E(50)`

``````require 'bigdecimal/math'
balance = BigMath.E(50) - 1
1.upto(25){|y| balance = balance * y - 1}
puts "Bank balance after 25 years = #{balance.to_f}

=> Bank balance after 25 years = 0.03993872967323021
``````

I don’t know how to do that in Crystal.

``````require "big"
balance = Math::E.to_big_f - 1
1.upto(25){ |y| balance = (balance * y).to_big_f - 1 }
puts "Bank balance after 25 years = #{balance.to_f}"

=> Bank balance after 25 years = -2242373258.5701575
``````

`#quo` returns a `Rational` on Ruby, if you use `to_big_r` then the first example will converge:

``````ar = [0, 2, -4].map(&.to_big_r)
100.times { ar << (111 - 1130.to_big_r / ar[-1] + 3000.to_big_r / (ar[-1] * ar[-2])) }
"%.16f" % ar[100] # => 6.0000000193194776
``````

For the second example your best bet is to compute `E` yourself, in more or less the same way Ruby’s `BigMath` does it:

``````def e(precision)
y = BigDecimal.new(10.to_big_i ** precision, precision)
d = y
i = 1

while true
d = (d / i).scale_to(y)
y2 = y + d
return y if y2 == y
y = y2
i += 1
end
end

balance = e - 1
1.upto(25) { |y| balance = (balance * y) - 1 }
balance # => 0.03993872967323020890366866978532990188388352
``````

This means neither Ruby snippet exercises any kind of floating-point arithmetic (not sure if that’s the point of the original problem).

1 Like

Thanks @HertzDevil, that worked.

The last exercise was much more straightforward.

Ruby

``````def rump(a,b)
a, b = a.to_r, b.to_r
333.75r * b**6 + a**2 * (11 * a**2 * b**2 - b**6 - 121 * b**4 - 2)  + 5.5r * b**8 + a / (2 * b)
end

puts "rump(77617, 33096) = #{rump(77617, 33096).to_f}"
=> -0.8273960599468214
``````

Crystal

``````require "big"

def rump(a,b)
a, b = a.to_big_r, b.to_big_r
333.75.to_big_r * b**6 + a**2 * (11 * a**2 * b**2 - b**6 - 121 * b**4 - 2)  + 5.5.to_big_r * b**8 + a / (2 * b)
end

puts "rump(77617, 33096) = #{rump(77617, 33096).to_f}"
=> -0.8273960599468213
``````