Pathological floating point problems: Difference between revisions
Pathological floating point problems (view source)
Revision as of 01:07, 28 August 2022
, 1 year agosyntax highlighting fixup automation
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
|||
Line 81:
in precision. A lot of more precise formats have been added since.<br>
'''A sequence that seems to converge to a wrong limit'''<br>
<
PATHOFP CSECT
USING PATHOFP,R13
Line 133:
YREGS
YFPREGS
END PATHOFP</
The divergence comes very soon.
{{out}}
Line 159:
===Task 1: Converging Sequence===
<
procedure Converging_Sequence is
Line 204:
Task_With_Short;
Task_With_Long;
end Converging_Sequence;</
{{out}}
Line 233:
===Task 2: Chaotic Bank Society===
<
procedure Chaotic_Bank is
Line 268:
Task_With_Short;
Task_With_Long;
end Chaotic_Bank;</
{{out}}
Line 309:
===Task 3: Rump's Example===
<
procedure Rumps_example is
Line 327:
Put("Rump's Example, Long: ");
LIO.Put(C, Fore => 1, Aft => 16, Exp => 0); New_Line;
end Rumps_example; </
{{out}}
Line 337:
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
In Algol 68G, we can specify the precision of LONG LONG REAL values
<
# task 1 #
BEGIN
Line 383:
print( ( "25 year chaotic balance: ", fixed( chaotic balance, -22, 16 ), newline ) )
END
END</
{{out}}
<pre>
Line 419:
awk code:
<
BEGIN {
do_task1()
Line 461:
}
</syntaxhighlight>
This version doesn't include the arbitrary-precision libraries, so the program demonstrates the incorrect results:
Line 508:
===First two tasks===
{{libheader|GMP}}
<syntaxhighlight lang="c">
#include<stdio.h>
#include<gmp.h>
Line 603:
return 0;
}
</syntaxhighlight>
The reason I included the trivial case was the discovery that the value of 0.3 is stored inexactly even by GMP if 0.1 and 0.2 are set via the mpf_set_d function, a point observed also during the solution of the [[Currency]] task. Thus mpf_set_str has been used to set the values, for the 2nd task, a great learning was not to convert the values into floating points unless they have to be printed out. It's a polynomial with out a single floating point coefficient or exponent, a clear sign that the values must be treated as rationals for the highest accuracy. Thus there are two implementations for this task in the above code, one uses floating points, the other rationals. There is still a loss of accuracy even when floating points are used, probably during the conversion of a rational to a float.
<pre>
Line 643:
The following sections source code must be located in a single file.
<
#define BANDED_ROWS
#define INCREASED_LIMITS
Line 711:
SetConsoleFormat(0);
}
}</
===Task 1: Converging sequence===
'''See''' [[#VB.NET Task 1]]
<
{
public static IEnumerable<float> SequenceSingle()
Line 833:
}
}
}</
{{out}}
Line 853:
'''See''' [[#VB.NET Task 2]]
<
{
public static IEnumerable<float> ChaoticBankSocietySingle()
Line 946:
}
}
}</
{{out}}
Line 995:
'''See''' [[#VB.NET Task 3]]
<
{
public static float SiegfriedRumpSingle(float a, float b)
Line 1,116:
Console.WriteLine($" Exact: {br}");
}
}</
{{out}}
Line 1,136:
In Clojure, rational numbers are a first class data type! This allow us to avoid these floating point calculation problems. As long as the operations all involve integers, rational numbers will automatically be used behind the scenes. If for some twisted reason you really wanted to introduce the error, you could change a number to a floating point in the equation to force floating point calculations instead.
===Task 1: Converging Sequence===
<
(def values [3 4 5 6 7 8 20 30 50 100])
; print decimal values:
(pprint (sort (zipmap values (map double (map #(nth converge-to-six (dec %)) values)))))
; print rational values:
(pprint (sort (zipmap values (map #(nth converge-to-six (dec %)) values))))</
{{out}}
Line 1,173:
===Task 2: Chaotic Bank Society===
<
(defn bank [n m] (- (* n m) 1))
(double (reduce bank (- e-ratio 1) (range 1 26)))</
{{out}}
Line 1,187:
===Task 3: Rump's Example===
<
(+ (* (rationalize 333.75) (expt b 6))
(* (expt a 2)
Line 1,194:
(/ a (* 2 b))))
; Using BigInt numeric literal style to avoid integer overflow
(double (rump 77617 33096N))</
{{out}}
Line 1,205:
===Task 1: Muller's sequence===
A) Using BigDecimal "div" (similar to Ruby's "quo" method), default iterations 100, here 132 min for last stable output.
<
ar = [0.to_big_d, 2.to_big_d, -4.to_big_d]
Line 1,214:
puts "%3d -> %0.16f" % [n, ar[n]]
end
</syntaxhighlight>
{{Out}}
<pre> 3 -> 18.5000000000000000
Line 1,229:
B) Using BigRationals.
<
ar = [0, 2, -4].map(&.to_big_r)
Line 1,238:
puts "%3d -> %0.16f" % [n, ar[n]]
end
</syntaxhighlight>
{{Out}}
<pre> 3 -> 18.5000000000000000
Line 1,254:
===Task 2: Chaotic Bank Society===
C) Unlike Ruby, had to manually create "E" to sufficient precision.
<
def e(precision)
Line 1,272:
balance = e(50) - 1
1.upto(25) { |y| balance = (balance * y) - 1 }
puts "Bank balance after 25 years = #{balance.to_f}"</
{{Out}}
<pre>Bank balance after 25 years = 0.03993872967323021</pre>
D) Or from Factor, use large rational approximation for "E = 106246577894593683 / 39085931702241241".
<
e = 106246577894593683.to_big_d.div(39085931702241241.to_big_d)
Line 1,283:
balance = e - 1
1.upto(25) { |y| balance = (balance * y) - 1 }
puts "Bank balance after 25 years = #{balance.to_f}"</
{{Out}}
<pre>Bank balance after 25 years = 0.03993873004901714</pre>
Line 1,289:
===Task 3: Rump's example===
Using BigRationals.
<
def rump(a, b)
Line 1,296:
end
puts "rump(77617, 33096) = #{rump(77617, 33096).to_f}"</
{{out}}<pre>rump(77617, 33096) = -0.8273960599468213
</pre>
Line 1,306:
{{libheader| Velthuis.BigDecimals}}
{{Trans|Go}}
<syntaxhighlight lang="delphi">
program Pathological_floating_point_problems;
Line 1,455:
readln;
end.</
{{out}}
<pre>Precision: 100
Line 1,490:
{{works with|Excel 2003|Excel 2015}}
'''A sequence that seems to converge to a wrong limit'''<br>
<syntaxhighlight lang="text"> A1: 2
A2: -4
A3: =111-1130/A2+3000/(A2*A1)
A4: =111-1130/A3+3000/(A3*A2)
...</
The result converges to the wrong limit!
{{out}}
Line 1,538:
=={{header|Factor}}==
These problems are straightforward due to Factor's rational numbers. One needs only take care not to introduce floating point values to the calculations.
<
math.ranges sequences ;
IN: rosetta-code.pathological
Line 1,568:
77617 33096 f "77617 33096 f = %.16f\n" printf ;
MAIN: pathological-demo</
{{out}}
<pre>
Line 1,600:
Here, no such attempt is made. In the spirit of Formula Translation, this is a direct translation of the specified formulae into Fortran, with single and double precision results on display. There is no REAL*16 option, nor the REAL*10 that some systems allow to correspond to the eighty-bit floating-point format supported by the floating-point processor. The various integer constants cause no difficulty and I'm not bothering with writing them as <integer>.0 - the compiler can deal with this. The constants with fractional parts happen to be exactly represented in binary so there is no fuss over 333.75 and 333.75D0 whereas by contrast 0.1 and 0.1D0 are not equal. Similarly, there is no attempt to rearrange the formulae, for instance to have <code>A**2 * B**2</code> replaced by <code>(A*B)**2</code>, nor worry over <code>B**8</code> where 33096**8 = 1.439E36 and the largest possible single-precision number is 3.4028235E+38, in part because arithmetic within an expression can be conducted with a greater dynamic range. Most of all, no attention has been given to the subtractions...
This would be F77 style Fortran, except for certain conveniences offered by F90, especially the availability of generic functions such as EXP whose type is determined by the type of its parameter, rather than having to use EXP and DEXP for single and double precision respectively, or else... The END statement for subroutines and functions names the routine being ended, a useful matter to have checked. <
REAL*4 VN,VNL1,VNL2 !The exact precision and dynamic range
REAL*8 WN,WNL1,WNL2 !Depends on the format's precise usage of bits.
Line 1,669:
WRITE (6,*) "Single precision:",SR4(77617.0,33096.0)
WRITE (6,*) "Double precision:",SR8(77617.0D0,33096.0D0) !Must match the types.
END</
====Output====
Line 1,792:
A simple function CBSERIES handles the special case deposit. The only question is how many terms of the series are required to produce a value accurate to the full precision in use. Thanks to the enquiry function EPSILON(x) offered by F90, the smallest number such that ''1 + eps'' differs from ''1'' for the precision of ''x'' is available without the need for cunning programming; this is a constant. An alternative form might be that EPSILON(X) returned the smallest number that, added to X, produced a result different from X in floating-point arithmetic of the precision of X - but this would not be a constant. Since the terms of the series are rapidly diminishing (and all are positive) a new term may be too small to affect the sum; this happens when S + T = S, or 1 + T/S = 1 + eps, thus the test in CBSERIES of T/S >= EPSILON(S) checks that the term affected the sum so that the loop stops for the first term that does not.
A misthimk had TINY(S) instead of EPSILON(S), and this demonstrates again the importance of providing output that shows the actual behaviour of a scheme and comparing it to expectations, since it showed that over a hundred terms were being calculated and the last term was tiny. Routine TINY(S) reports the smallest possible floating-point number in the precision of its parameter, which is not what is wanted! EPSILON(S) is tiny, but not so tiny as TINY(S). 2·220446049250313E-016 instead of 2·225073858507201E-308. <
INTEGER YEAR !A stepper.
REAL*4 V !The balance.
Line 1,827:
CBSERIES = S !Convergence is ever-faster as N increases.
END FUNCTION CBSERIES !So this is easy.
END SUBROUTINE CBS !Madness! </
And the output is (slightly decorated to show correct digits in bold):
Line 1,861:
=={{header|FreeBASIC}}==
<
' As FB's native types have only 64 bit precision at most we need to use the
Line 1,958:
Print
Print "Press any key to quit"
Sleep</
{{out}}
Line 1,987:
=={{header|Go}}==
<
import (
Line 2,080:
f64, _ := s.Add(s.Add(s, t3), t4).Float64()
return f64
}</
{{out}}
<pre>
Line 2,102:
Task 1 includes an extra step, 200 iterations, to demonstrate a closer convergence.
<
# Pathological floating point problems
#
Line 2,189:
show := 10^4
write("Balance after ", y, " years: $", d * show / scale * 1.0 / show)
end</
{{out}}
Line 2,235:
Implementation of <code>vn</code>:
<
Example using IEEE-754 floating point:
<
3 9.3783783783783861
4 7.8011527377522611
Line 2,249:
30 100.0000000000000000
50 100.0000000000000000
100 100.0000000000000000</
Example using exact arithmetic:
<
3 9.3783783783783784
4 7.8011527377521614
Line 2,263:
30 6.0056486887714203
50 6.0001465345613879
100 6.0000000160995649</
'''The Chaotic Bank Society'''
Line 2,269:
Let's start this example by using exact arithmetic, to make sure we have the right algorithm. We'll go a bit overboard, in representing ''e'', so we don't have to worry too much about that.
<
81j76":e
2.7182818284590452353602874713526624977572470936999595749669676277240766303535
21j16":+`*/,_1,.(1+i.-25),e
0.0399387296732302</
(Aside: here, we are used the same mechanism for adding -1 to ''e'' that we are using to add -1 to the product of the year number and the running balance.)
Line 2,279:
Next, we will use <math>6157974361033 \div 2265392166685</math> for ''e'', to represent the limit of what can be expressed using 64 bit IEEE 754 floating point.
<
_2053975868590.1852178761057505</
That's clearly way too low, so let's try instead using <math>(1+6157974361033) \div 2265392166685</math> for ''e''
<
4793054977300.3491517765983265</
So, our problem seems to be that there's no way we can express enough bits of ''e'', using 64 bit IEEE-754 floating point arithmetic. Just to confirm:
<
2.71828
+`*/,_1,.(1+i.-25),1x1
_2.24237e9</
Now let's take a closer look using our rational approximation for ''e'':
<
0.0399387296732302
21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.30x
Line 2,304:
0.0000000000000000
21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.24x
_1.0000000000000000</
Things go haywire when our approximation for ''e'' uses the same number of terms as our bank's term. So, what does that look like, in terms of precision?
<
2.718281828459045235360287471257428715
41j36":+/%1,*/\1+i.25x
2.718281828459045235360287468777832452
41j36":+/%1,*/\1+i.24x
2.718281828459045235360287404308329608</
In other words, we go astray when our approximation for ''e'' is inaccurate in the 26th position after the decimal point. But IEEE-754 floating point arithmetic can only represent approximately 16 decimal digits of precision.
Line 2,321:
Again, we use exact arithmetic to see if we have the algorithm right. That said, we'll also do this in small steps, to make sure we're being exact every step of the way, and to keep from building overly long lines:
<
NB. enforce exact arithmetic
add=. +&x:
Line 2,351:
term1 add term2 add term3 add term4
)</
Example use:
<
_0.8273960599468214</
Note that replacing the definitions of <code>add</code>, <code>sub</code>, <code>div</code>, <code>mul</code> with implementations which promote to floating point gives a very different result:
<
_1.39061e21</
But given that b8 is
<
1.43947e36</
we're exceeding the limits of our representation here, if we're using 64 bit IEEE-754 floating point arithmetic.
Line 2,373:
Uses BigRational class: [[Arithmetic/Rational/Java]]. For comparison purposes, each task is also implemented with standard 64-bit floating point numbers.
<
import java.math.RoundingMode;
Line 2,469:
siegfriedRump();
}
}</
{{out}}
<pre>Wrong Convergence Sequence
Line 2,522:
===v series===
The following implementation illustrates how a cache can be used in jq to avoid redundant computations. A JSON object is used as the cache.
<
def v:
# Input: cache
Line 2,538:
end
end;
. as $m | {} | v_($m) | .[($m|tostring)] ; </
Example:<
{{out}}
<pre>18.5
Line 2,557:
To avoid the pathological issues, the following uses symbolic arithmetic, with {"e": m, "c": n} representing (e*m + n).
<
# Input: { e: m, c: n } representing m*e + n
def new_balance(n):
Line 2,567:
def e: 1|exp;
reduce range(0;n) as $i ({}; new_balance($i) )
| (.e * e) + .c;</
Example:<syntaxhighlight lang
{{out}}
0
Line 2,578:
The task could be completed even using ```Rational{BigFloat}``` type.
<
# arbitrary precision
setprecision(2000)
Line 2,616:
333.75b ^ 6 + a ^ 2 * ( 11a ^ 2 * b ^ 2 - b ^ 6 - 121b ^ 4 - 2 ) + 5.5b ^ 8 + a / 2b
println("\nTask 3 - Siegfried Rump's example:\nf(77617.0, 33096.0) = ", @sprintf "%.20f" f(big(77617.0), big(33096.0)))</
{{out}}
Line 2,634:
=={{header|Kotlin}}==
<
import java.math.*
Line 2,675:
f += c8 * a.pow(2, con480)
println("\nf(77617.0, 33096.0) is ${"%18.16f".format(f)}")
}</
{{out}}
Line 2,785:
=={{header|Mathematica}}/{{header|Wolfram Language}}==
Task 1:
<
v[2] = -4;
v[n_] := Once[111 - 1130/v[n - 1] + 3000/(v[n - 1]*v[n - 2])]
N[Map[v, {3, 4, 5, 6, 7, 8, 20, 30, 50, 100}], 80]</
{{out}}
<pre>{18.500000000000000000000000000000000000000000000000000000000000000000000000000000,
Line 2,802:
Task 2:
<
{{out}}<pre>0.0399387296732302089036714552104</pre>
Task 3:
<
f[77617, 33096]</
{{out}}<pre>-0.827396059946821368141165095479816291999033115784384819917814842</pre>
Line 2,824:
To avoid duplication of code, we used generics.
<
import decimal
import bignum
Line 2,938:
echo &"f({A}, {B}) float = ", rumpFloat
echo &"f({A}, {B}) decimal = ", rumpDecimal.format(1, 16)
echo &"f({A}, {B}) rational = ", rumpRational.format(1, 16)</
{{out}}
Line 2,967:
Task 1: Define recursive function V(n):
<
In order to work set precision to at least 200 digits:
<pre>\p 200: realprecision = 211 significant digits (200 digits displayed)
Line 2,976:
----
Task 2: Define function balance(deposit,years):
<
Output ''balance(exp(1), 25)'':
Line 2,982:
----
Task 3: Define function f(a,b):
<
Output:<pre>f(77617.0,33096.0): -0.827396059946821368141165...</pre>
Line 2,992:
The constants in the equation must be represented with a decimal point (even just <code>111.</code>) so that
they are treated as rationals, not integers.
<
@s = qw(2, -4);
Line 3,002:
($nu,$de) = $s[$n-1] =~ m#^(\d+)/(\d+)#;;
printf "n = %3d %18.15f\n", $n, $nu/$de;
}</
{{out}}
<pre>n = 3 18.500000000000000
Line 3,019:
The value of 𝑒 is imported from a module, but could also be calculated, as is done in
[[Calculating_the_value_of_e#Perl|Calculating the value of e]]
<
$e = bexp(1,43);
$years = 25;
Line 3,026:
print "Starting balance, \$(e-1): \$$balance\n";
for $i (1..$years) { $balance = $i * $balance - 1 }
printf "After year %d, you will have \$%1.15g in your account.\n", $years, $balance;</
{{out}}
<pre>Starting balance, $(e-1): $1.718281828459045235360287471352662497757247
Line 3,032:
==== Rump's example ====
<
$a = 77617;
$b = 33096;
printf "%0.16f\n", 333.75*$b**6 + $a**2 * ( 11*$a**2 * $b**2 - $b**6 - 121*$b**4 - 2) + 5.5*$b**8 + $a/(2*$b);</
{{out}}
<pre>-0.8273960599468214</pre>
Line 3,045:
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.
<!--<
<span style="color: #008080;">include</span> <span style="color: #000000;">builtins</span><span style="color: #0000FF;">\</span><span style="color: #000000;">bigatom</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Task 1\n"</span><span style="color: #0000FF;">)</span>
Line 3,090:
<span style="color: #000000;">f_ab</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ba_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ba_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ba_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pb633375</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pa2mid</span><span style="color: #0000FF;">),</span><span style="color: #000000;">pb855</span><span style="color: #0000FF;">),</span><span style="color: #000000;">ba_divide</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ba_multiply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)))</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"f(%d, %d) = %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ba_sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%.15B"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">f_ab</span><span style="color: #0000FF;">)})</span>
<!--</
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}}
Line 3,121:
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).
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.0"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (mpfr_set_default_prec[ision] has been renamed)</span>
Line 3,199:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"f(%d, %d) = %s\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">15</span><span style="color: #0000FF;">)})</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">t1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t7</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_free</span><span style="color: #0000FF;">({</span><span style="color: #000000;">t1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t7</span><span style="color: #0000FF;">})</span>
<!--</
Identical output
=={{header|PicoLisp}}==
<
(de task1 (N)
(cache '(NIL) N
Line 3,253:
(*/ 5.5 B8 1.0)
(*/ 1.0 A (*/ 2.0 B 1.0)) ) ) )
(prinl "Rump's example: " (round (task3 77617.0 33096.0) 20))</
{{out}}
<pre>
Line 3,276:
Using rational numbers via standard library <code>fractions</code>
<
def muller_seq(n:int) -> float:
Line 3,287:
for n in [3, 4, 5, 6, 7, 8, 20, 30, 50, 100]:
print("{:4d} -> {}".format(n, muller_seq(n)))</
{{Out}}
Line 3,306:
Using <code>decimal</code> numbers with a high precision
<
def bank(years:int) -> float:
Line 3,321:
return(float(decimal_balance))
print("Bank balance after 25 years = ", bank(25))</
{{Out}}
Line 3,328:
'''but, still incorrectly diverging after some time, aprox. 250 years'''
<
print(year, '->', bank(year))
</syntaxhighlight>
{{Out}}
Line 3,351:
Using rational numbers via standard library <code>fractions</code>
<
def rump(generic_a, generic_b) -> float:
Line 3,362:
print("rump(77617, 33096) = ", rump(77617.0, 33096.0))
</syntaxhighlight>
{{Out}}
Line 3,375:
The examples below use real numbers, and the <code>x</code> function is used to transform them to floats, if desired, with the function <code>exact->inexact</code>.
<
(define current-do-exact-calculations? (make-parameter exact->inexact))
Line 3,422:
(displayln "EXACT (Rational) NUMBERS")
(parameterize ([current-do-exact-calculations? #t])
(all-tests)))</
{{out}}
Line 3,462:
{{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. :-) Raku 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.
<syntaxhighlight lang="raku"
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]"};
Line 3,476:
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");</
{{Out}}
<pre>1st: Convergent series
Line 3,505:
<br>calculations, as well how many fractional decimal digits (past the decimal point) should be displayed.
===A sequence that seems to converge to a wrong limit===
<
parse arg digs show . /*obtain optional arguments from the CL*/
if digs=='' | digs=="," then digs= 150 /*Not specified? Then use the default.*/
Line 3,519:
/*display digs past the dec. point───┐ */
if wordpos(n, #)\==0 then say 'v.'left(n, w) "=" format(v.n, , show)
end /*n*/ /*stick a fork in it, we're all done. */</
{{out|output|text= when using the default inputs:}}
<pre>
Line 3,543:
With 150 decimal digits, the balance becomes negative after '''96''' years.
<
e=2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713,
||8217852516642742746639193200305992181741359662904357290033429526059563073813232862794,
Line 3,560:
end /*n*/
@@@= 'With ' d " decimal digits, the balance after " y ' years is: '
say @@@ '$'format($, , show) / 1 /*stick a fork in it, we're all done. */</
{{out|output|text= when using the default inputs:}}
<pre>
Line 3,567:
===Siegfried Rump's example===
<
parse arg digs show . /*obtain optional arguments from the CL*/
if digs=='' | digs=="," then digs=150 /*Not specified? Then use the default.*/
Line 3,579:
/*──────────────────────────────────────────────────────────────────────────────────────*/
f: procedure; parse arg a,b; return 333.75* b**6 + a**2 * (11* a**2* b**2 - b**6,
- 121*b**4 - 2) + 5.5*b**8 + a / (2*b)</
{{out|output|text= when using the default inputs:}}
<pre>
Line 3,586:
=={{header|Ring}}==
<
# Project : Pathological floating point problems
Line 3,599:
ok
next
</syntaxhighlight>
Output:
<pre>
Line 3,617:
===Task 1: Muller's sequence===
Ruby numbers have a "quo" division method, which returns a rational (a fraction) when possible, avoiding Float inaccuracy.
<
100.times{ar << (111 - 1130.quo(ar[-1])+ 3000.quo(ar[-1]*ar[-2])) }
Line 3,623:
puts "%3d -> %0.16f" % [n, ar[n]]
end
</syntaxhighlight>
{{Out}}
<pre> 3 -> 18.5000000000000000
Line 3,638:
===Task 2: The Chaotic Bank Society===
Using BigDecimal provides a way to specify the number of digits for E. 50 seems to be sufficient.
<
balance = BigMath.E(50) - 1
1.upto(25){|y| balance = balance * y - 1}
puts "Bank balance after 25 years = #{balance.to_f}"</
{{Out}}
<pre>Bank balance after 25 years = 0.03993872967323021
Line 3,648:
===Task 3: Rump's example===
Rationals again.
<
a, b = a.to_r, b.to_r
333.75r * b**6 + a**2 * ( 11 * a**2 * b**2 - b**6 - 121 * b**4 - 2 ) + 5.5r * b**8 + a / (2 * b)
end
puts "rump(77617, 33096) = #{rump(77617, 33096).to_f}"</
{{out}}<pre>rump(77617, 33096) = -0.8273960599468214
</pre>
Line 3,659:
=={{header|Sidef}}==
'''Muller's sequence'''
<
var (u, v) = (2, -4)
(n-2).times { (u, v) = (v, 111 - 1130/v + 3000/(v * u)) }
Line 3,667:
[(3..8)..., 20, 30, 50, 100].each {|n|
printf("n = %3d -> %s\n", n, series(n))
}</
{{out}}
<pre>
Line 3,683:
'''The Chaotic Bank Society'''
<
var balance = (1 .. years+15 -> sum_by {|n| 1 / n! })
say "Starting balance, $(e-1): $#{balance}"
for i in (1..years) { balance = (i*balance - 1) }
printf("After year %d, you will have $%1.16g in your account.\n", years, balance)</
{{out}}
<pre>
Line 3,695:
'''Siegfried Rump's example'''
<
(333.75 * b**6) + (a**2 * ((11 * a**2 * b**2) -
b**6 - (121 * b**4) - 2)) + (5.5 * b**8) + a/(2*b)
}
say f(77617.0, 33096.0)</
{{out}}
<pre>
Line 3,708:
=={{header|Stata}}==
'''Task 1'''
<
set obs 100
gen n=_n
Line 3,716:
replace v=111-1130/v[_n-1]+3000/(v[_n-1]*v[_n-2]) in 3/l
format %20.16f v
list if inlist(n,3,4,5,6,7,8,20,30,50,100), noobs</
'''Output'''
Line 3,739:
'''Task 2'''
<
set obs 26
gen year=_n-1
gen balance=exp(1)-1 in 1
replace balance=year*balance[_n-1]-1 in 2/l
list balance if year==25, noobs</
'''Output'''
Line 3,758:
If we replace exp(1)-1 by its value computed in higher precision by other means, the result is different, owing to the extreme sensibility of the computation.
<
set obs 26
gen year=_n-1
gen balance=1.71828182845904523 in 1
replace balance=year*balance[_n-1]-1 in 2/l
list balance if year==25, noobs</
'''Output'''
Line 3,777:
We can check the hexadecimal representation of both numbers. Note they differ only in the last bit:
<
di %21x 1.71828182845904523</
'''Output'''
Line 3,792:
Using mkrd's Swift-BigInt library.
<
@inlinable
public func power(_ n: Self) -> Self {
Line 3,888:
print("\nRump's function")
print("rump(77617.0, 33096.0) as Float \(rumpFloat); as Double = \(rumpDouble); as Decimal = \(rumpDecimal); as BDouble = \(rumpBigDouble.decimalExpansion(precisionAfterDecimalPoint: 16, rounded: false))")</
{{out}}
Line 3,918:
u(1)=2
u(2)=-4
<
u(n)=111-1130/u(n-1) + 3000/(u(n-1)*u(n-2))
u(nMin)={-4;2}</
The result converges to the wrong limit!
{{out}}
Line 3,959:
Main() procedure and output formatting:
<
Imports System.Numerics
Imports Numerics
Line 4,012:
SetConsoleFormat(0)
End Sub
End Module</
<!--Anchors for the C# section to link to-->
Line 4,018:
Somewhat predictably, Single fairs the worst and Double slightly better. Decimal lasts past iteration 20, and BigRational remains exactly precise.
<
Iterator Function SequenceSingle() As IEnumerable(Of Single)
' n, n-1, and n-2
Line 4,124:
Next
End Sub
End Module</
{{out}}
Line 4,145:
Decimal appears to be doing well by year 25, but begins to degenerate ''the very next year'' and soon overflows.
<
Iterator Function ChaoticBankSocietySingle() As IEnumerable(Of Single)
Dim balance As Single = Math.E - 1
Line 4,247:
Next
End Sub
End Module</
{{out}}
Line 4,300:
Each of these three tasks are susceptible to this phenomenon; the outputs of this program for floating-point arithmetic should not be expected to be the same on different systems.
<
Function SiegfriedRumpSingle(a As Single, b As Single) As Single
Dim a2 = a * a,
Line 4,408:
Console.WriteLine($" Exact: {br}")
End Sub
End Module</
{{out|note=debug build}}
Line 4,442:
{{libheader|Wren-big}}
{{libheader|Wren-fmt}}
<
import "/fmt" for Fmt
Line 4,476:
var c8 = c5 * a.pow(2) * b.pow(2) - b.pow(6) - c6 * b.pow(4) - 2
f = f + c8 * a.pow(2)
System.print("\nf(77617.0, 33096.0) is %(f.toDecimal(16))")</
{{out}}
Line 4,586:
=={{header|zkl}}==
zkl doesn't have a big rational or big float library (as of this writing) but does have big ints (via GNU GMP). It does have 64 bit doubles.
<
vs.append(111.0 - 1130.0/vs[-1] + 3000.0/(vs[-1]*vs[-2])) }.fp(List(2,-4)));
series:=Series.drop(100).value;</
We'll use the convenient formula given in the referenced paper to create a fraction with big ints
<
fcn u(n){ // use formula to create a fraction of big ints
Line 4,608:
foreach n in (T(3,4,5,6,7,8,20,30,50,100)){
"n =%3d; %3.20F %s".fmt(n,series[n-1],u(n-1)).println();
}</
{{out}}
Note that, at n=100, we still have diverged (at the 15th place) from the Raku solution and 12th place from the J solution.
Line 4,625:
</pre>
Chaotic banking society is just nasty so we use a five hundred digit e (the e:= text is one long line).
<
e:="271828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723069697720931014169283681902551510865746377211125238978442505695369677078544996996794686445490598793163688923009879312";
var en=(e.len()-1), tenEN=BN(10).pow(en);
Line 4,631:
balance=[1..years].reduce(fcn(balance,i){ balance*i - tenEN },balance);
balance=tostr(balance,en,2);
println("After year %d, you will have $%s in your account.".fmt(years,balance));</
{{out}}
<pre>
Line 4,638:
</pre>
For Rump's example, multiple the formula by 10ⁿ so we can use integer math.
<
b2,b4,b6,b8:=b.pow(2),b.pow(4),b.pow(6),b.pow(8);
a2:=BN(a).pow(2);
Line 4,645:
tostr(r,66,32)
}
println("\n3rd: Rump's example: f(77617.0, 33096.0) = ",rump(77617,33096));</
{{out}}
<pre>
|