Rosetta Code

Pathological floating point problems: Difference between revisions

Page Discussion

Pathological floating point problems (view source)

Revision as of 17:09, 14 March 2020

17 bytes added ,  4 years ago
Rename Perl 6 -> Raku, alphabetize, minor clean-up
(Rename Perl 6 -> Raku, alphabetize, minor clean-up)
Line 1,303:
77617 33096 f = -0.8273960599468214
</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 &mdash;i.e. XML, JSON&mdash; 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}}==
Line 1,699 ⟶ 1,573:
24**************** 48072289.9364179447293282 0.0415975491869292 0.16679020D-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 &mdash;i.e. XML, JSON&mdash; 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}}==
Line 2,600:
{{out}}
<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}}==
Line 2,983 ⟶ 2,945:
TASK 2: balance after 25 years = 0.039938729673230209
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}}==
Line 3,143 ⟶ 3,144:
{{out}}<pre>rump(77617, 33096) = -0.8273960599468214
</pre>
 
=={{header|Sidef}}==
'''Muller's sequence'''
Thundergnat
10,327

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/264655"