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

Balance: Num := Ada.Numerics.E - 1.0;

begin

Integer'Image(After) & " digits");

Ada.Text_IO.Put_Line("year balance");

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>

v_100 : 6.000000019319477929060000000000000000000000000000000000000000000000000000000000

</pre>

 

=={{header|Clojure}}==

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

</pre>

 

 

 

 

 

 

 

 

 

 

{{out}}

<pre>

 

 

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

24**************** 48072289.9364179447293282 0.0415975491869292 0.16679020D-17 12

25**************** 1201807247.4104485511779785 0.0399387296732302 0.11269608D-17 12

 

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

 

=={{header|jq}}==

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

 

=={{header|Julia}}==

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

 

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>

 

 

 

n = 100 6.000000019319477929104086803403585715024350675436952458072592750856521767230266

 

 

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

TASK 2: balance after 25 years = 0.039938729673230209

TASK 3: f(77617, 33096) = -54767/66192 = -0.827396059946821368</pre>

 

=={{header|REXX}}==

 

A little extra boilerplate code was added to support the specification of how many decimal digits that should be used for the

===A sequence that seems to converge to a wrong limit===

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

numeric digits digs /*have REXX use "digs" decimal digits. */

#= 2 4 5 6 7 8 9 20 30 50 100 /*the indices to display value of V.n */

v.1= 2 /*the value of the first V element. */

v.2=-4 /* " " " " second " " */

<pre>

v.4 = 9.37837837837837837838

 

===The Chaotic Bank Society===

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

if show=='' | show=="," then show= 20 /* " " " " " " */

if y=='' | y=="," then y= 25 /* " " " " " " */

numeric digits digs /*have REXX use "digs" decimal digits. */

/* [↑] value of newly opened account. */

do n=1 for y /*compute the value of the account/year*/

end /*n*/

<pre>

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

exit /*stick a fork in it, we're all done. */

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

<pre>

f(a,b)= -0.82739605994682136814

 

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

u(50) : 100

u(100) : 100

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

<pre>

1st: Convergent series

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