Pathological floating point problems: Difference between revisions

m (→‎{{header|C}}: Remove vanity tags)
Line 1,876:
 
3rd: Rump's example: f(77617.0, 33096.0) = -0.827396059946821
</pre>
 
=={{header|Phix}}==
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>
Task2: needs at least 41 decimal places (and 42 passed as the first argument to ba_euler)<br>
Task3: apparently only needs just 15 decimal places, then again ba_scale() defines the minimum, and it may (as in is permitted to) use more.
<lang Phix>include builtins\bigatom.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}})
{} = ba_scale(196)
sequence v = {2,-4}
for n=3 to 100 do
-- v = append(v,111 - 1130/v[n-1] + 3000/(v[n-1]*v[n-2]))
v = append(v,ba_add(ba_sub(111,ba_divide(1130,v[n-1])),ba_divide(3000,ba_multiply(v[n-1],v[n-2]))))
if n<9 or find(n,{20,30,50,100}) then
-- printf(1,"n = %-3d %20.16f\n", {n, v[n]})
string fmt = sprintf("%%.%dB",fmts[find(n,fns)])
printf(1,"n = %-3d %s\n", {n, ba_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)
{} = ba_scale(41)
bigatom balance = ba_sub(ba_euler(42,true),1)
for i=1 to 25 do balance = ba_sub(ba_multiply(balance,i),1) end for
ba_printf(1,"Balance after 25 years: $%.16B\n\n", balance)
 
puts(1,"Task 3\n")
{} = ba_scale(15) -- fine!
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})
bigatom pa2 = ba_power(a,2),
pb2a211 = ba_multiply(11,ba_multiply(pa2,ba_power(b,2))),
pb4121 = ba_multiply(121,ba_power(b,4)),
pb6 = ba_power(b,6),
pa2mid = ba_multiply(pa2,ba_sub(ba_sub(ba_sub(pb2a211,pb6),pb4121),2)),
pb633375 = ba_multiply(333.75,pb6),
pb855 = ba_multiply(5.5,ba_power(b,8)),
f_ab = ba_add(ba_add(ba_add(pb633375,pa2mid),pb855),ba_divide(a,ba_multiply(2,b)))
printf(1,"f(%d, %d) = %s\n", {a, b, ba_sprintf("%.15B",f_ab)})</lang>
Aside: obviously you don't ''have'' to use lots of temps like that, but I find that style often makes things much easier to debug.
{{out}}
<pre>
Task 1
n = 3 18.5
n = 4 9.378378
n = 5 7.801153
n = 6 7.154414
n = 7 6.806785
n = 8 6.5926328
n = 20 6.0435521101892689
n = 30 6.006786093031205758530554
n = 50 6.0001758466271871889456140207471954695237
n = 100 6.000000019319477929104086803403585715024350675436952458072592750856521767230266
 
Task 2
Balance after 25 years: $0.0399387296732302
 
Task 3
f(77617, 33096) = -0.827396059946821
</pre>
 
7,806

edits