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).

2 Likes

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