Pathological floating point problems: Difference between revisions
Pathological floating point problems (view source)
Revision as of 15:17, 4 August 2019
, 4 years ago→{{header|Phix}}: added gmp version
m (→{{header|Phix}}: added gmp version) |
|||
Line 2,087:
=={{header|Phix}}==
{{libheader|bigatom}}▼
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>
Line 2,157 ⟶ 2,156:
f(77617, 33096) = -0.827396059946821
</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}}==
| |||