Pathological floating point problems: Difference between revisions
Content added Content deleted
m (→{{header|Phix}}: added gmp version) |
|||
Line 2,087: | Line 2,087: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
⚫ | |||
Standard maths with the IEEE-754 hardware floats fails predictably (that code left in as comments), so roll out the bigatoms.<br> |
Standard maths with the IEEE-754 hardware floats fails predictably (that code left in as comments), so roll out the bigatoms.<br> |
||
Task1: needs at least 120 decimal places to avoid serious divergance, and 196 to be accurate to 78 digits<br> |
Task1: needs at least 120 decimal places to avoid serious divergance, and 196 to be accurate to 78 digits<br> |
||
Line 2,157: | Line 2,156: | ||
f(77617, 33096) = -0.827396059946821 |
f(77617, 33096) = -0.827396059946821 |
||
</pre> |
</pre> |
||
=== gmp version === |
|||
⚫ | |||
Note that bigatom is formally deprecated in favour of gmp. However there are significant benefits to the |
|||
former on this task, being easier to read and debug.<br> |
|||
Again, I have left native/hopelessly incorrect versions in as comments for comparison.<br> |
|||
Task1: needs at least 120 decimal places to avoid serious divergance, and 196 to be accurate to 78 digits.<br> |
|||
Task2: needs at least 41 decimal places (and e-1 accurately specified to at least 42 decimal places).<br> |
|||
Task3: needs at least 36 decimal places (my suspicions above re bigatom in 15 now seem correct). |
|||
<lang Phix>include builtins\mpfr.e |
|||
puts(1,"Task 1\n") |
|||
constant {fns,fmts} = columnize({{3,1},{4,6},{5,6},{6,6},{7,6},{8,7},{20,16},{30,24},{50,40},{100,78}}) |
|||
mpfr_set_default_prec(-196) |
|||
sequence v = {mpfr_init(2),mpfr_init(-4)} |
|||
mpfr t1 = mpfr_init() |
|||
for n=3 to 100 do |
|||
-- v = append(v,111 - 1130/v[n-1] + 3000/(v[n-1]*v[n-2])) |
|||
mpfr_set_si(t1,1130) |
|||
mpfr_div(t1,t1,v[n-1]) |
|||
mpfr_si_sub(t1,111,t1) |
|||
mpfr t2 = mpfr_init() |
|||
mpfr_mul(t2,v[n-1],v[n-2]) |
|||
mpfr_si_div(t2,3000,t2) |
|||
mpfr_add(t2,t1,t2) |
|||
v = append(v,t2) |
|||
if n<9 or find(n,{20,30,50,100}) then |
|||
-- printf(1,"n = %-3d %20.16f\n", {n, v[n]}) |
|||
string fmt = sprintf("%%.%dRf",fmts[find(n,fns)]) |
|||
printf(1,"n = %-3d %s\n", {n, mpfr_sprintf(fmt,v[n])}) |
|||
end if |
|||
end for |
|||
puts(1,"\nTask 2\n") |
|||
--atom balance = exp(1)-1 |
|||
--for i=1 to 25 do balance = balance*i-1 end for |
|||
--printf(1,"\nTask 2\nBalance after 25 years: $%12.10f", balance) |
|||
mpfr_set_default_prec(-41) |
|||
mpfr balance = mpfr_init("1.71828182845904523536028747135266249775724709369995"& |
|||
"95749669676277240766303535475945713821785251664274") |
|||
for i=1 to 25 do |
|||
mpfr_mul_si(balance,balance,i) |
|||
mpfr_sub_si(balance,balance,1) |
|||
end for |
|||
mpfr_printf(1,"Balance after 25 years: $%.16Rf\n\n", balance) |
|||
puts(1,"Task 3\n") |
|||
mpfr_set_default_prec(-36) |
|||
integer a = 77617, |
|||
b = 33096 |
|||
--atom pa2 = power(a,2), |
|||
-- pb2a211 = 11*pa2*power(b,2), |
|||
-- pb4121 = 121*power(b,4), |
|||
-- pb6 = power(b,6), |
|||
-- pb855 = 5.5*power(b,8), |
|||
-- f_ab = 333.75 * pb6 + pa2 * (pb2a211 - pb6 - pb4121 - 2) + pb855 + a/(2*b) |
|||
--printf(1,"f(%d, %d) = %.15f\n\n", {a, b, f_ab}) |
|||
-- (translation of FreeBASIC) |
|||
mpfr {t2,t3,t4,t5,t6,t7} = mpfr_inits(6) |
|||
mpfr_set_d(t1, a) -- a |
|||
mpfr_set_d(t2, b) -- b |
|||
mpfr_set_d(t3, 333.75) -- 333.75 |
|||
mpfr_pow_si(t4, t2, 6) -- b ^ 6 |
|||
mpfr_mul(t3, t3, t4) -- 333.75 * b^6 |
|||
mpfr_pow_si(t5, t1, 2) -- a^2 |
|||
mpfr_pow_si(t6, t2, 2) -- b^2 |
|||
mpfr_mul_si(t7, t5, 11) -- 11 * a^2 |
|||
mpfr_mul(t7, t7, t6) -- 11 * a^2 * b^2 |
|||
mpfr_sub(t7, t7, t4) -- 11 * a^2 * b^2 - b^6 |
|||
mpfr_pow_si(t4, t2, 4) -- b^4 |
|||
mpfr_mul_si(t4, t4, 121) -- 121 * b^4 |
|||
mpfr_sub(t7, t7, t4) -- 11 * a^2 * b^2 - b^6 - 121 * b^4 |
|||
mpfr_sub_si(t7, t7, 2) -- 11 * a^2 * b^2 - b^6 - 121 * b^4 - 2 |
|||
mpfr_mul(t7, t7, t5) -- (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 |
|||
mpfr_add(t3, t3, t7) -- 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 |
|||
mpfr_set_d(t4, 5.5) -- 5.5 |
|||
mpfr_pow_si(t5, t2, 8) -- b^8 |
|||
mpfr_mul(t4, t4, t5) -- 5.5 * b^8 |
|||
mpfr_add(t3, t3, t4) -- 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 |
|||
mpfr_mul_si(t4, t2, 2) -- 2 * b |
|||
mpfr_div(t5, t1, t4) -- a / (2 * b) |
|||
mpfr_add(t3, t3, t5) -- 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 + a / (2 * b) |
|||
printf(1,"f(%d, %d) = %s\n", {a, b, mpfr_sprintf("%.15Rf",t3)}) |
|||
{t1,t2,t3,t4,t5,t6,t7} = mpfr_free({t1,t2,t3,t4,t5,t6,t7})</lang> |
|||
Identical output |
|||
=={{header|Python}}== |
=={{header|Python}}== |