Pathological floating point problems: Difference between revisions
Content added Content deleted
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
|||
Line 1,303: | Line 1,303: | ||
77617 33096 f = -0.8273960599468214 |
77617 33096 f = -0.8273960599468214 |
||
</pre> |
</pre> |
||
=={{header|FreeBASIC}}== |
|||
<lang freebasic>' FB 1.05.0 Win64 |
|||
' As FB's native types have only 64 bit precision at most we need to use the |
|||
' C library, GMP v6.1.0, for arbitrary precision arithmetic |
|||
#Include Once "gmp.bi" |
|||
mpf_set_default_prec(640) '' 640 bit precision, enough for this exercise |
|||
Function v(n As UInteger, prev As __mpf_struct, prev2 As __mpf_struct) As __mpf_struct |
|||
Dim As __mpf_struct a, b, c |
|||
mpf_init(@a) : mpf_init(@b) : mpf_init(@c) |
|||
If n = 0 Then mpf_set_ui(@a, 0UL) |
|||
If n = 1 Then mpf_set_ui(@a, 2UL) |
|||
If n = 2 Then mpf_set_si(@a, -4L) |
|||
If n < 3 Then Return a |
|||
mpf_ui_div(@a, 1130UL, @prev) |
|||
mpf_mul(@b, @prev, @prev2) |
|||
mpf_ui_div(@c, 3000UL, @b) |
|||
mpf_ui_sub(@b, 111UL, @a) |
|||
mpf_add(@a, @b, @c) |
|||
mpf_clear(@b) |
|||
mpf_clear(@c) |
|||
Return a |
|||
End Function |
|||
Function f(a As Double, b As Double) As __mpf_Struct |
|||
Dim As __mpf_struct temp1, temp2, temp3, temp4, temp5, temp6, temp7, temp8 |
|||
mpf_init(@temp1) : mpf_init(@temp2) : mpf_init(@temp3) : mpf_init(@temp4) |
|||
mpf_init(@temp5) : mpf_init(@temp6) : mpf_init(@temp7) : mpf_init(@temp8) |
|||
mpf_set_d(@temp1, a) '' a |
|||
mpf_set_d(@temp2, b) '' b |
|||
mpf_set_d(@temp3, 333.75) '' 333.75 |
|||
mpf_pow_ui(@temp4, @temp2, 6UL) '' b ^ 6 |
|||
mpf_mul(@temp3, @temp3, @temp4) '' 333.75 * b^6 |
|||
mpf_pow_ui(@temp5, @temp1, 2UL) '' a^2 |
|||
mpf_pow_ui(@temp6, @temp2, 2UL) '' b^2 |
|||
mpf_mul_ui(@temp7, @temp5, 11UL) '' 11 * a^2 |
|||
mpf_mul(@temp7, @temp7, @temp6) '' 11 * a^2 * b^2 |
|||
mpf_sub(@temp7, @temp7, @temp4) '' 11 * a^2 * b^2 - b^6 |
|||
mpf_pow_ui(@temp4, @temp2, 4UL) '' b^4 |
|||
mpf_mul_ui(@temp4, @temp4, 121UL) '' 121 * b^4 |
|||
mpf_sub(@temp7, @temp7, @temp4) '' 11 * a^2 * b^2 - b^6 - 121 * b^4 |
|||
mpf_sub_ui(@temp7, @temp7, 2UL) '' 11 * a^2 * b^2 - b^6 - 121 * b^4 - 2 |
|||
mpf_mul(@temp7, @temp7, @temp5) '' (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 |
|||
mpf_add(@temp3, @temp3, @temp7) '' 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 |
|||
mpf_set_d(@temp4, 5.5) '' 5.5 |
|||
mpf_pow_ui(@temp5, @temp2, 8UL) '' b^8 |
|||
mpf_mul(@temp4, @temp4, @temp5) '' 5.5 * b^8 |
|||
mpf_add(@temp3, @temp3, @temp4) '' 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 |
|||
mpf_mul_ui(@temp4, @temp2, 2UL) '' 2 * b |
|||
mpf_div(@temp5, @temp1, @temp4) '' a / (2 * b) |
|||
mpf_add(@temp3, @temp3, @temp5) '' 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 + a / (2 * b) |
|||
mpf_clear(@temp1) : mpf_clear(@temp2) : mpf_clear(@temp4) : mpf_clear(@temp5) |
|||
mpf_clear(@temp6) : mpf_clear(@temp7) : mpf_clear(@temp8) |
|||
Return temp3 |
|||
End Function |
|||
Dim As Zstring * 60 z |
|||
Dim As __mpf_struct result, prev, prev2 |
|||
' We cache the two previous results to avoid recursive calls to v |
|||
For i As Integer = 1 To 100 |
|||
result = v(i, prev, prev2) |
|||
If (i >= 3 AndAlso i <= 8) OrElse i = 20 OrElse i = 30 OrElse i = 50 OrElse i = 100 Then |
|||
gmp_sprintf(@z,"%53.50Ff",@result) '' express result to 50 decimal places |
|||
Print "n ="; i , z |
|||
End If |
|||
prev2 = prev |
|||
prev = result |
|||
Next |
|||
mpf_clear(@prev) : mpf_clear(@prev2) '' note : prev = result |
|||
Dim As __mpf_struct e, balance, ii, temp |
|||
mpf_init(@e) : mpf_init(@balance) : mpf_init(@ii) : mpf_init(@temp) |
|||
mpf_set_str(@e, "2.71828182845904523536028747135266249775724709369995", 10) '' e to 50 decimal places |
|||
mpf_sub_ui(@balance, @e, 1UL) |
|||
For i As ULong = 1 To 25 |
|||
mpf_set_ui(@ii, i) |
|||
mpf_mul(@temp, @balance, @ii) |
|||
mpf_sub_ui(@balance, @temp, 1UL) |
|||
Next |
|||
Print |
|||
Print "Chaotic B/S balance after 25 years : "; |
|||
gmp_sprintf(@z,"%.16Ff",@balance) '' express balance to 16 decimal places |
|||
Print z |
|||
mpf_clear(@e) : mpf_clear(@balance) : mpf_clear(@ii) : mpf_clear(@temp) |
|||
Print |
|||
Dim rump As __mpf_struct |
|||
rump = f(77617.0, 33096.0) |
|||
gmp_sprintf(@z,"%.16Ff", @rump) '' express rump to 16 decimal places |
|||
Print "f(77617.0, 33096.0) = "; z |
|||
Print |
|||
Print "Press any key to quit" |
|||
Sleep</lang> |
|||
{{out}} |
|||
<pre> |
|||
n = 3 18.50000000000000000000000000000000000000000000000000 |
|||
n = 4 9.37837837837837837837837837837837837837837837837838 |
|||
n = 5 7.80115273775216138328530259365994236311239193083573 |
|||
n = 6 7.15441448097524935352789065386036202438123383819727 |
|||
n = 7 6.80678473692363298394175659627200908762327670780193 |
|||
n = 8 6.59263276870443839274200277636599482655298231773461 |
|||
n = 20 6.04355211018926886777747736409754013318771500000612 |
|||
n = 30 6.00678609303120575853055404795323970583307231443837 |
|||
n = 50 6.00017584662718718894561402074719546952373517709933 |
|||
n = 100 6.00000001931947792910408680340358571502435067543695 |
|||
Chaotic B/S balance after 25 years : 0.0399387296732302 |
|||
f(77617.0, 33096.0) = -0.8273960599468214 |
|||
</pre> |
|||
=={{header|Fōrmulæ}}== |
|||
In [http://wiki.formulae.org/Pathological_floating_point_problems this] page you can see the solution of this task. |
|||
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text ([http://wiki.formulae.org/Editing_F%C5%8Drmul%C3%A6_expressions more info]). Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for transportation effects more than visualization and edition. |
|||
The option to show Fōrmulæ programs and their results is showing images. Unfortunately images cannot be uploaded in Rosetta Code. |
|||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
Line 1,699: | Line 1,573: | ||
24**************** 48072289.9364179447293282 0.0415975491869292 0.16679020D-17 12 |
24**************** 48072289.9364179447293282 0.0415975491869292 0.16679020D-17 12 |
||
25**************** 1201807247.4104485511779785 0.0399387296732302 0.11269608D-17 12 |
25**************** 1201807247.4104485511779785 0.0399387296732302 0.11269608D-17 12 |
||
=={{header|FreeBASIC}}== |
|||
<lang freebasic>' FB 1.05.0 Win64 |
|||
' As FB's native types have only 64 bit precision at most we need to use the |
|||
' C library, GMP v6.1.0, for arbitrary precision arithmetic |
|||
#Include Once "gmp.bi" |
|||
mpf_set_default_prec(640) '' 640 bit precision, enough for this exercise |
|||
Function v(n As UInteger, prev As __mpf_struct, prev2 As __mpf_struct) As __mpf_struct |
|||
Dim As __mpf_struct a, b, c |
|||
mpf_init(@a) : mpf_init(@b) : mpf_init(@c) |
|||
If n = 0 Then mpf_set_ui(@a, 0UL) |
|||
If n = 1 Then mpf_set_ui(@a, 2UL) |
|||
If n = 2 Then mpf_set_si(@a, -4L) |
|||
If n < 3 Then Return a |
|||
mpf_ui_div(@a, 1130UL, @prev) |
|||
mpf_mul(@b, @prev, @prev2) |
|||
mpf_ui_div(@c, 3000UL, @b) |
|||
mpf_ui_sub(@b, 111UL, @a) |
|||
mpf_add(@a, @b, @c) |
|||
mpf_clear(@b) |
|||
mpf_clear(@c) |
|||
Return a |
|||
End Function |
|||
Function f(a As Double, b As Double) As __mpf_Struct |
|||
Dim As __mpf_struct temp1, temp2, temp3, temp4, temp5, temp6, temp7, temp8 |
|||
mpf_init(@temp1) : mpf_init(@temp2) : mpf_init(@temp3) : mpf_init(@temp4) |
|||
mpf_init(@temp5) : mpf_init(@temp6) : mpf_init(@temp7) : mpf_init(@temp8) |
|||
mpf_set_d(@temp1, a) '' a |
|||
mpf_set_d(@temp2, b) '' b |
|||
mpf_set_d(@temp3, 333.75) '' 333.75 |
|||
mpf_pow_ui(@temp4, @temp2, 6UL) '' b ^ 6 |
|||
mpf_mul(@temp3, @temp3, @temp4) '' 333.75 * b^6 |
|||
mpf_pow_ui(@temp5, @temp1, 2UL) '' a^2 |
|||
mpf_pow_ui(@temp6, @temp2, 2UL) '' b^2 |
|||
mpf_mul_ui(@temp7, @temp5, 11UL) '' 11 * a^2 |
|||
mpf_mul(@temp7, @temp7, @temp6) '' 11 * a^2 * b^2 |
|||
mpf_sub(@temp7, @temp7, @temp4) '' 11 * a^2 * b^2 - b^6 |
|||
mpf_pow_ui(@temp4, @temp2, 4UL) '' b^4 |
|||
mpf_mul_ui(@temp4, @temp4, 121UL) '' 121 * b^4 |
|||
mpf_sub(@temp7, @temp7, @temp4) '' 11 * a^2 * b^2 - b^6 - 121 * b^4 |
|||
mpf_sub_ui(@temp7, @temp7, 2UL) '' 11 * a^2 * b^2 - b^6 - 121 * b^4 - 2 |
|||
mpf_mul(@temp7, @temp7, @temp5) '' (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 |
|||
mpf_add(@temp3, @temp3, @temp7) '' 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 |
|||
mpf_set_d(@temp4, 5.5) '' 5.5 |
|||
mpf_pow_ui(@temp5, @temp2, 8UL) '' b^8 |
|||
mpf_mul(@temp4, @temp4, @temp5) '' 5.5 * b^8 |
|||
mpf_add(@temp3, @temp3, @temp4) '' 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 |
|||
mpf_mul_ui(@temp4, @temp2, 2UL) '' 2 * b |
|||
mpf_div(@temp5, @temp1, @temp4) '' a / (2 * b) |
|||
mpf_add(@temp3, @temp3, @temp5) '' 333.75 * b^6 + (11 * a^2 * b^2 - b^6 - 121 * b^4 - 2) * a^2 + 5.5 * b^8 + a / (2 * b) |
|||
mpf_clear(@temp1) : mpf_clear(@temp2) : mpf_clear(@temp4) : mpf_clear(@temp5) |
|||
mpf_clear(@temp6) : mpf_clear(@temp7) : mpf_clear(@temp8) |
|||
Return temp3 |
|||
End Function |
|||
Dim As Zstring * 60 z |
|||
Dim As __mpf_struct result, prev, prev2 |
|||
' We cache the two previous results to avoid recursive calls to v |
|||
For i As Integer = 1 To 100 |
|||
result = v(i, prev, prev2) |
|||
If (i >= 3 AndAlso i <= 8) OrElse i = 20 OrElse i = 30 OrElse i = 50 OrElse i = 100 Then |
|||
gmp_sprintf(@z,"%53.50Ff",@result) '' express result to 50 decimal places |
|||
Print "n ="; i , z |
|||
End If |
|||
prev2 = prev |
|||
prev = result |
|||
Next |
|||
mpf_clear(@prev) : mpf_clear(@prev2) '' note : prev = result |
|||
Dim As __mpf_struct e, balance, ii, temp |
|||
mpf_init(@e) : mpf_init(@balance) : mpf_init(@ii) : mpf_init(@temp) |
|||
mpf_set_str(@e, "2.71828182845904523536028747135266249775724709369995", 10) '' e to 50 decimal places |
|||
mpf_sub_ui(@balance, @e, 1UL) |
|||
For i As ULong = 1 To 25 |
|||
mpf_set_ui(@ii, i) |
|||
mpf_mul(@temp, @balance, @ii) |
|||
mpf_sub_ui(@balance, @temp, 1UL) |
|||
Next |
|||
Print |
|||
Print "Chaotic B/S balance after 25 years : "; |
|||
gmp_sprintf(@z,"%.16Ff",@balance) '' express balance to 16 decimal places |
|||
Print z |
|||
mpf_clear(@e) : mpf_clear(@balance) : mpf_clear(@ii) : mpf_clear(@temp) |
|||
Print |
|||
Dim rump As __mpf_struct |
|||
rump = f(77617.0, 33096.0) |
|||
gmp_sprintf(@z,"%.16Ff", @rump) '' express rump to 16 decimal places |
|||
Print "f(77617.0, 33096.0) = "; z |
|||
Print |
|||
Print "Press any key to quit" |
|||
Sleep</lang> |
|||
{{out}} |
|||
<pre> |
|||
n = 3 18.50000000000000000000000000000000000000000000000000 |
|||
n = 4 9.37837837837837837837837837837837837837837837837838 |
|||
n = 5 7.80115273775216138328530259365994236311239193083573 |
|||
n = 6 7.15441448097524935352789065386036202438123383819727 |
|||
n = 7 6.80678473692363298394175659627200908762327670780193 |
|||
n = 8 6.59263276870443839274200277636599482655298231773461 |
|||
n = 20 6.04355211018926886777747736409754013318771500000612 |
|||
n = 30 6.00678609303120575853055404795323970583307231443837 |
|||
n = 50 6.00017584662718718894561402074719546952373517709933 |
|||
n = 100 6.00000001931947792910408680340358571502435067543695 |
|||
Chaotic B/S balance after 25 years : 0.0399387296732302 |
|||
f(77617.0, 33096.0) = -0.8273960599468214 |
|||
</pre> |
|||
=={{header|Fōrmulæ}}== |
|||
In [http://wiki.formulae.org/Pathological_floating_point_problems this] page you can see the solution of this task. |
|||
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text ([http://wiki.formulae.org/Editing_F%C5%8Drmul%C3%A6_expressions more info]). Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for transportation effects more than visualization and edition. |
|||
The option to show Fōrmulæ programs and their results is showing images. Unfortunately images cannot be uploaded in Rosetta Code. |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
Line 2,600: | Line 2,600: | ||
{{out}} |
{{out}} |
||
<pre>-0.8273960599468214</pre> |
<pre>-0.8273960599468214</pre> |
||
=={{header|Perl 6}}== |
|||
{{works with|Rakudo|2016-01}} |
|||
The simple solution to doing calculations where floating point numbers exhibit pathological behavior is: don't do floating point calculations. :-) Perl 6 is just as susceptible to floating point error as any other C based language, however, it offers built-in rational Types; where numbers are represented as a ratio of two integers. For normal precision it uses Rats - accurate to 1/2^64, and for arbitrary precision, FatRats, whose denominators can grow as large as available memory. Rats don't require any special setup to use. Any decimal number within its limits of precision is automatically stored as a Rat. FatRats require explicit coercion and are "sticky". Any FatRat operand in a calculation will cause all further results to be stored as FatRats. |
|||
<lang perl6>say '1st: Convergent series'; |
|||
my @series = 2.FatRat, -4, { 111 - 1130 / $^v + 3000 / ( $^v * $^u ) } ... *; |
|||
for flat 3..8, 20, 30, 50, 100 -> $n {say "n = {$n.fmt("%3d")} @series[$n-1]"}; |
|||
say "\n2nd: Chaotic bank society"; |
|||
sub postfix:<!> (Int $n) { [*] 2..$n } # factorial operator |
|||
my $years = 25; |
|||
my $balance = sum map { 1 / FatRat.new($_!) }, 1 .. $years + 15; # Generate e-1 to sufficient precision with a Taylor series |
|||
put "Starting balance, \$(e-1): \$$balance"; |
|||
for 1..$years -> $i { $balance = $i * $balance - 1 } |
|||
printf("After year %d, you will have \$%1.16g in your account.\n", $years, $balance); |
|||
print "\n3rd: Rump's example: f(77617.0, 33096.0) = "; |
|||
sub f (\a, \b) { 333.75*b⁶ + a²*( 11*a²*b² - b⁶ - 121*b⁴ - 2 ) + 5.5*b⁸ + a/(2*b) } |
|||
say f(77617.0, 33096.0).fmt("%0.16g");</lang> |
|||
{{Out}} |
|||
<pre>1st: Convergent series |
|||
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 |
|||
2nd: Chaotic bank society |
|||
Starting balance, $(e-1): $1.7182818284590452353602874713526624977572470936999 |
|||
After year 25, you will have $0.0399387296732302 in your account. |
|||
3rd: Rump's example: f(77617.0, 33096.0) = -0.827396059946821 |
|||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Line 2,983: | Line 2,945: | ||
TASK 2: balance after 25 years = 0.039938729673230209 |
TASK 2: balance after 25 years = 0.039938729673230209 |
||
TASK 3: f(77617, 33096) = -54767/66192 = -0.827396059946821368</pre> |
TASK 3: f(77617, 33096) = -54767/66192 = -0.827396059946821368</pre> |
||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
{{works with|Rakudo|2016-01}} |
|||
The simple solution to doing calculations where floating point numbers exhibit pathological behavior is: don't do floating point calculations. :-) Perl 6 is just as susceptible to floating point error as any other C based language, however, it offers built-in rational Types; where numbers are represented as a ratio of two integers. For normal precision it uses Rats - accurate to 1/2^64, and for arbitrary precision, FatRats, whose denominators can grow as large as available memory. Rats don't require any special setup to use. Any decimal number within its limits of precision is automatically stored as a Rat. FatRats require explicit coercion and are "sticky". Any FatRat operand in a calculation will cause all further results to be stored as FatRats. |
|||
<lang perl6>say '1st: Convergent series'; |
|||
my @series = 2.FatRat, -4, { 111 - 1130 / $^v + 3000 / ( $^v * $^u ) } ... *; |
|||
for flat 3..8, 20, 30, 50, 100 -> $n {say "n = {$n.fmt("%3d")} @series[$n-1]"}; |
|||
say "\n2nd: Chaotic bank society"; |
|||
sub postfix:<!> (Int $n) { [*] 2..$n } # factorial operator |
|||
my $years = 25; |
|||
my $balance = sum map { 1 / FatRat.new($_!) }, 1 .. $years + 15; # Generate e-1 to sufficient precision with a Taylor series |
|||
put "Starting balance, \$(e-1): \$$balance"; |
|||
for 1..$years -> $i { $balance = $i * $balance - 1 } |
|||
printf("After year %d, you will have \$%1.16g in your account.\n", $years, $balance); |
|||
print "\n3rd: Rump's example: f(77617.0, 33096.0) = "; |
|||
sub f (\a, \b) { 333.75*b⁶ + a²*( 11*a²*b² - b⁶ - 121*b⁴ - 2 ) + 5.5*b⁸ + a/(2*b) } |
|||
say f(77617.0, 33096.0).fmt("%0.16g");</lang> |
|||
{{Out}} |
|||
<pre>1st: Convergent series |
|||
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 |
|||
2nd: Chaotic bank society |
|||
Starting balance, $(e-1): $1.7182818284590452353602874713526624977572470936999 |
|||
After year 25, you will have $0.0399387296732302 in your account. |
|||
3rd: Rump's example: f(77617.0, 33096.0) = -0.827396059946821 |
|||
</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
Line 3,143: | Line 3,144: | ||
{{out}}<pre>rump(77617, 33096) = -0.8273960599468214 |
{{out}}<pre>rump(77617, 33096) = -0.8273960599468214 |
||
</pre> |
</pre> |
||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
'''Muller's sequence''' |
'''Muller's sequence''' |