jzakiya
November 16, 2021, 4:34am
1
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.
jzakiya
November 16, 2021, 7:10pm
4
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
jzakiya
November 16, 2021, 9:21pm
5
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
jzakiya
November 17, 2021, 8:38am
7
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