Pathological floating point problems: Difference between revisions

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

YREGS

YFPREGS

The divergence comes very soon.

{{out}}

===Task 1: Converging Sequence===

 

 

procedure Converging_Sequence is

Task_With_Short;

Task_With_Long;

 

{{out}}

===Task 2: Chaotic Bank Society===

 

 

procedure Chaotic_Bank is

Task_With_Short;

Task_With_Long;

 

{{out}}

===Task 3: Rump's Example===

 

procedure Rumps_example is

Put("Rump's Example, Long: ");

LIO.Put(C, Fore => 1, Aft => 16, Exp => 0); New_Line;

 

{{out}}

{{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

print( ( "25 year chaotic balance: ", fixed( chaotic balance, -22, 16 ), newline ) )

END

{{out}}

<pre>

 

awk code:

BEGIN {

do_task1()

}

 

 

This version doesn't include the arbitrary-precision libraries, so the program demonstrates the incorrect results:

===First two tasks===

{{libheader|GMP}}

#include<stdio.h>

#include<gmp.h>

return 0;

}

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>

The following sections source code must be located in a single file.

 

#define BANDED_ROWS

#define INCREASED_LIMITS

SetConsoleFormat(0);

}

 

===Task 1: Converging sequence===

'''See''' [[#VB.NET Task 1]]

 

{

public static IEnumerable<float> SequenceSingle()

}

}

 

{{out}}

'''See''' [[#VB.NET Task 2]]

 

{

public static IEnumerable<float> ChaoticBankSocietySingle()

}

}

 

{{out}}

'''See''' [[#VB.NET Task 3]]

 

{

public static float SiegfriedRumpSingle(float a, float b)

Console.WriteLine($" Exact: {br}");

}

 

{{out}}

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:

 

{{out}}

 

===Task 2: Chaotic Bank Society===

(defn bank [n m] (- (* n m) 1))

 

{{out}}

 

===Task 3: Rump's Example===

(+ (* (rationalize 333.75) (expt b 6))

(* (expt a 2)

(/ a (* 2 b))))

; Using BigInt numeric literal style to avoid integer overflow

 

{{out}}

-0.8273960599468214

</pre>

 

=={{header|Excel}}==

{{works with|Excel 2003|Excel 2015}}

'''A sequence that seems to converge to a wrong limit'''<br>

A2: -4

A3: =111-1130/A2+3000/(A2*A1)

A4: =111-1130/A3+3000/(A3*A2)

The result converges to the wrong limit!

{{out}}

=={{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

77617 33096 f "77617 33096 f = %.16f\n" printf ;

 

{{out}}

<pre>

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...

 

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.

WRITE (6,*) "Single precision:",SR4(77617.0,33096.0)

WRITE (6,*) "Double precision:",SR8(77617.0D0,33096.0D0) !Must match the types.

 

====Output====

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.

 

INTEGER YEAR !A stepper.

REAL*4 V !The balance.

CBSERIES = S !Convergence is ever-faster as N increases.

END FUNCTION CBSERIES !So this is easy.

 

And the output is (slightly decorated to show correct digits in bold):

 

=={{header|FreeBASIC}}==

 

' As FB's native types have only 64 bit precision at most we need to use the

Print

Print "Press any key to quit"

 

{{out}}

=={{header|Fōrmulæ}}==

 

 

 

 

=={{header|Go}}==

 

import (

f64, _ := s.Add(s.Add(s, t3), t4).Float64()

return f64

{{out}}

<pre>

Task 1 includes an extra step, 200 iterations, to demonstrate a closer convergence.

 

# Pathological floating point problems

#

show := 10^4

write("Balance after ", y, " years: $", d * show / scale * 1.0 / show)

 

{{out}}

Implementation of <code>vn</code>:

 

 

Example using IEEE-754 floating point:

 

3 9.3783783783783861

4 7.8011527377522611

30 100.0000000000000000

50 100.0000000000000000

 

Example using exact arithmetic:

 

3 9.3783783783783784

4 7.8011527377521614

30 6.0056486887714203

50 6.0001465345613879

 

'''The Chaotic Bank Society'''

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

 

(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.)

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.

 

 

That's clearly way too low, so let's try instead using <math>(1+6157974361033) \div 2265392166685</math> for ''e''

 

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

 

Now let's take a closer look using our rational approximation for ''e'':

 

0.0399387296732302

21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.30x

0.0000000000000000

21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.24x

 

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

 

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.

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:

 

term1 add term2 add term3 add term4

 

Example use:

 

 

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:

 

 

But given that b8 is

 

 

we're exceeding the limits of our representation here, if we're using 64 bit IEEE-754 floating point arithmetic.

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;

 

siegfriedRump();

}

{{out}}

<pre>Wrong Convergence Sequence

===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

end

end;

 

{{out}}

<pre>18.5

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):

def e: 1|exp;

reduce range(0;n) as $i ({}; new_balance($i) )

 

{{out}}

0

The task could be completed even using ```Rational{BigFloat}``` type.

 

# arbitrary precision

setprecision(2000)

333.75b ^ 6 + a ^ 2 * ( 11a ^ 2 * b ^ 2 - b ^ 6 - 121b ^ 4 - 2 ) + 5.5b ^ 8 + a / 2b

 

 

{{out}}

 

=={{header|Kotlin}}==

 

import java.math.*

f += c8 * a.pow(2, con480)

println("\nf(77617.0, 33096.0) is ${"%18.16f".format(f)}")

 

{{out}}

</pre>

 

 

Task 1:

v[2] = -4;

v[n_] := Once[111 - 1130/v[n - 1] + 3000/(v[n - 1]*v[n - 2])]

{{out}}

<pre>{18.500000000000000000000000000000000000000000000000000000000000000000000000000000,

 

Task 2:

{{out}}<pre>0.0399387296732302089036714552104</pre>

 

Task 3:

{{out}}<pre>-0.827396059946821368141165095479816291999033115784384819917814842</pre>

 

=={{header|PARI/GP}}==

 

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)

----

Task 2: Define function balance(deposit,years):

 

Output ''balance(exp(1), 25)'':

----

Task 3: Define function f(a,b):

 

Output:<pre>f(77617.0,33096.0): -0.827396059946821368141165...</pre>

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);

($nu,$de) = $s[$n-1] =~ m#^(\d+)/(\d+)#;;

printf "n = %3d %18.15f\n", $n, $nu/$de;

{{out}}

<pre>n = 3 18.500000000000000

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;

print "Starting balance, \$(e-1): \$$balance\n";

for $i (1..$years) { $balance = $i * $balance - 1 }

{{out}}

<pre>Starting balance, $(e-1): $1.718281828459045235360287471352662497757247

 

==== Rump's example ====

 

$a = 77617;

$b = 33096;

{{out}}

<pre>-0.8273960599468214</pre>

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.

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}}

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).

Identical output

 

=={{header|PicoLisp}}==

(de task1 (N)

(cache '(NIL) N

(*/ 5.5 B8 1.0)

(*/ 1.0 A (*/ 2.0 B 1.0)) ) ) )

{{out}}

<pre>

Using rational numbers via standard library <code>fractions</code>

 

 

def muller_seq(n:int) -> float:

 

for n in [3, 4, 5, 6, 7, 8, 20, 30, 50, 100]:

 

{{Out}}

Using <code>decimal</code> numbers with a high precision

 

 

def bank(years:int) -> float:

return(float(decimal_balance))

 

 

{{Out}}

 

'''but, still incorrectly diverging after some time, aprox. 250 years'''

print(year, '->', bank(year))

 

{{Out}}

Using rational numbers via standard library <code>fractions</code>

 

 

def rump(generic_a, generic_b) -> float:

 

print("rump(77617, 33096) = ", rump(77617.0, 33096.0))

 

{{Out}}

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))

(displayln "EXACT (Rational) NUMBERS")

(parameterize ([current-do-exact-calculations? #t])

 

{{out}}

{{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.

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]"};

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) }

{{Out}}

<pre>1st: Convergent series

<br>calculations, &nbsp; as well how many fractional decimal digits &nbsp; (past the decimal point) &nbsp; 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.*/

/*display digs past the dec. point───┐ */

if wordpos(n, #)\==0 then say 'v.'left(n, w) "=" format(v.n, , show)

{{out|output|text=&nbsp; when using the default inputs:}}

<pre>

With 150 decimal digits, the balance becomes negative after &nbsp; '''96''' &nbsp; years.

e=2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713,

||8217852516642742746639193200305992181741359662904357290033429526059563073813232862794,

||3490763233829880753195251019011573834187930702154089149934884167509244761460668082264,

||8001684774118537423454424371075390777449920695517027618386062613313845830007520449338

parse arg digs show y . /*obtain optional arguments from the CL*/

if digs=='' | digs=="," then digs= d /*Not specified? Then use the default.*/

end /*n*/

@@@= 'With ' d " decimal digits, the balance after " y ' years is: '

{{out|output|text=&nbsp; when using the default inputs:}}

<pre>

 

===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.*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

f: procedure; parse arg a,b; return 333.75* b**6 + a**2 * (11* a**2* b**2 - b**6,

{{out|output|text=&nbsp; when using the default inputs:}}

<pre>

 

=={{header|Ring}}==

# Project : Pathological floating point problems

 

ok

next

Output:

<pre>

===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])) }

puts "%3d -> %0.16f" % [n, ar[n]]

end

{{Out}}

<pre> 3 -> 18.5000000000000000

===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}

{{Out}}

<pre>Bank balance after 25 years = 0.03993872967323021

===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

 

{{out}}<pre>rump(77617, 33096) = -0.8273960599468214

</pre>

=={{header|Sidef}}==

'''Muller's sequence'''

var (u, v) = (2, -4)

(n-2).times { (u, v) = (v, 111 - 1130/v + 3000/(v * u)) }

[(3..8)..., 20, 30, 50, 100].each {|n|

printf("n = %3d -> %s\n", n, series(n))

{{out}}

<pre>

 

'''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) }

{{out}}

<pre>

 

'''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)

}

 

{{out}}

<pre>

=={{header|Stata}}==

'''Task 1'''

set obs 100

gen n=_n

replace v=111-1130/v[_n-1]+3000/(v[_n-1]*v[_n-2]) in 3/l

format %20.16f v

 

'''Output'''

'''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

 

'''Output'''

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

 

'''Output'''

We can check the hexadecimal representation of both numbers. Note they differ only in the last bit:

 

 

'''Output'''

+1.b7e151628aed3X+000

</pre>

 

=={{header|TI-83 BASIC}}==

u(1)=2

u(2)=-4

u(n)=111-1130/u(n-1) + 3000/(u(n-1)*u(n-2))

The result converges to the wrong limit!

{{out}}

Main() procedure and output formatting:

 

Imports System.Numerics

Imports Numerics

SetConsoleFormat(0)

End Sub

 

<!--Anchors for the C# section to link to-->

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

Next

End Sub

 

{{out}}

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

Next

End Sub

 

{{out}}

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,

Console.WriteLine($" Exact: {br}")

End Sub

 

{{out|note=debug build}}

{{libheader|Wren-big}}

{{libheader|Wren-fmt}}

 

var LIMIT = 100

var c8 = c5 * a.pow(2) * b.pow(2) - b.pow(6) - c6 * b.pow(4) - 2

f = f + c8 * a.pow(2)

 

{{out}}

 

f(77617.0, 33096.0) is -0.8273960599468214

<pre>

 

=={{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)));

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

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.

</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);

balance=[1..years].reduce(fcn(balance,i){ balance*i - tenEN },balance);

balance=tostr(balance,en,2);

{{out}}

<pre>

</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);

tostr(r,66,32)

}

{{out}}

<pre>