Pathological floating point problems: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎{{header|REXX}}: corrected typos, added a comment in the REXX section header.)
Line 305: Line 305:
f(a,b)= -0.82739605994682136814
f(a,b)= -0.82739605994682136814
</pre>
</pre>

=={{header|TI-83 BASIC}}==
Use the SEQ mode to enter the arithmetic progression. Note the way to set
u(1)=2
u(2)=-4
<lang ti83b>
nMin=1
u(n)=111-1130/u(n-1) + 3000/(u(n-1)*u(n-2))
u(nMin)={-4;2}
</lang>
The result is the wrong limit:
{{out}}
<pre>


u(20) : 100.055202
u(30) : 100
u(50) : 100
u(100) : 100

</pre>



=={{header|zkl}}==
=={{header|zkl}}==

Revision as of 11:43, 3 May 2016

Pathological floating point problems is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Most programmers are familiar with the inexactness of floating point calculations in a binary processor. The classic example being: 

0.1 + 0.2 =  0.30000000000000004

In many situations the amount of error in such calculations is very small and can be overlooked or eliminated with rounding.

There are pathological problems however, where seemingly simple, straight-forward calculations are extremely sensitive to even tiny amounts of imprecision.

Show how your language deals with such classes of problems.


A sequence that seems to converge to a wrong limit. Consider the sequence: 

   v1 = 2
   v2 = -4
   vn = 111 - 1130/vn-1 + 3000/(vn-1 * vn-2)

As n grows larger, the series should converge to 6 but small amounts of error will cause it to approach 100.

Display the values of the sequence where n = 3, 4, 5, 6, 7, 8, 20, 30, 50 & 100 to at least 16 decimal places. 

   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


The Chaotic Bank Society is offering a new investment account to their customers. You first deposit $e-1 where e is 2.7182818... the base of natural logarithms.

After each year, your account balance will be multiplied by the number of years that have passed, and $1 in service charges will be removed. So after 1 year, your balance will be multiplied by 1 and $1 will be removed for service charges, after 2 years your balance will be doubled and $1 removed, after 3 years your balance will be tripled and $1 removed, after 10 years, multiplied by 10 and $1 removed and so on. What will your balance be after 25 years? 

   Starting balance: $e-1
   Balance = (Balance * year) - 1 for 25 years
   Balance after 25 years: $0.0399387296732302


Siegfried Rump's example. Consider the following function, designed by Siegfried Rump in 1988. 

   f(a,b) = 333.75b⁶ + a²( 11a²b² - b⁶ - 121b⁴ - 2) + 5.5b⁸ + a/2b
   compute f(a,b) where a = 77617.0 and b = 33096.0.
   f(77617.0, 33096.0) = -0.827396059946821

Demonstrate how to solve at least one of the first two problems, or both, and the third if you're feeling particularly jaunty.


See also;



J

A sequence that seems to converge to a wrong limit.

Implementation of vn:

<lang J> vn=: 111 +(_1130 % _1&{) + (3000 % _1&{ * _2&{)</lang>

Example using IEEE-754 floating point:

<lang J> 3 21j16 ":"1] 3 4 5 6 7 8 20 30 50 100 ([,.{) (,vn)^:100(2 _4) 

 3   9.3783783783783861
 4   7.8011527377522611
 5   7.1544144809765555
 6   6.8067847369419638
 7   6.5926327689743687
 8   6.4494659378910058
20  99.9934721906444960
30 100.0000000000000000
50 100.0000000000000000

100 100.0000000000000000</lang>

Example using exact arithmetic:

<lang J> 3 21j16 ":"1] 3 4 5 6 7 8 20 30 50 100 ([,.{) (,vn)^:100(2 _4x) 

 3   9.3783783783783784
 4   7.8011527377521614
 5   7.1544144809752494
 6   6.8067847369236330
 7   6.5926327687044384
 8   6.4494659337902880
20   6.0360318810818568
30   6.0056486887714203
50   6.0001465345613879

100 6.0000000160995649</lang>

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.

<lang J> e=: +/%1,*/\1+i.100x 

  81j76":e
  2.7182818284590452353602874713526624977572470936999595749669676277240766303535
  21j16":+`*/,_1,.(1+i.-25),e
  0.0399387296732302</lang>

(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 for e, to represent the limit of what can be expressed using 64 bit IEEE 754 floating point.

<lang J> 31j16":+`*/,_1,.(1+i.-25),6157974361033%2265392166685x _2053975868590.1852178761057505</lang>

That's clearly way too low, so let's try instead using for e <lang J> 31j16":+`*/,_1,.(1+i.-25),6157974361034%2265392166685x 

4793054977300.3491517765983265</lang>

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:

<lang J> 1x1 2.71828 

  +`*/,_1,.(1+i.-25),1x1

_2.24237e9</lang>

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

<lang J> 21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.40x 

  0.0399387296732302
  21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.30x
  0.0399387277260840
  21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.26x
  0.0384615384615385
  21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.25x
  0.0000000000000000
  21j16":+`*/,_1,.(1+i.-25),+/%1,*/\1+i.24x
 _1.0000000000000000</lang>

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?

<lang J> 41j36":+/%1,*/\1+i.26x 

  2.718281828459045235360287471257428715
  41j36":+/%1,*/\1+i.25x
  2.718281828459045235360287468777832452
  41j36":+/%1,*/\1+i.24x
  2.718281828459045235360287404308329608</lang>

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.

Siegfried Rump's example.

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:

<lang J>rump=:4 :0 

 NB. enforce exact arithmetic
 add=. +&x:
 sub=. -&x:
 mul=. *&x:
 div=. %&x:

 a=. x
 a2=. a mul a

 b=. y
 b2=. b mul b
 b4=. b2 mul b2
 b6=. b2 mul b4
 b8=. b4 mul b4

 c333_75=. 1335 div 4 NB. 333.75
 term1=. c333_75 mul b6

 t11a2b2=. 11 mul a2 mul b2
 tnb6=. 0 sub b6
 tn121b4=. 0 sub 121 mul b4
 term2=. a2*(t11a2b2 + tnb6 + tn121b4 sub 2)

 c5_5=. 11 div 2 NB. 5.5
 term3=. c5_5 mul b8

 term4=. a div 2 mul b

 term1 add term2 add term3 add term4

)</lang>

Example use:

<lang J> 21j16": 77617 rump 33096 

 _0.8273960599468214</lang>

Note that replacing the definitions of add, sub, div, mul with implementations which promote to floating point gives a very different result:

<lang J> 77617 rump 33096 _1.39061e21</lang>

But given that b8 is

<lang J> 33096^8 1.43947e36</lang>

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

Perl 6

Works with: Rakudo version 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, which can grow as large as available memory. Rats don't require any special 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 banking 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>

Output:
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 banking 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

REXX

The REXX language uses character-based arithmetic.   So effectively, it looks, feels, and tastes like decimal floating point
(implemented in software).

So, the only (minor) problem is how many decimal digits should be used to solve these pathological floating point problems.

A little extra boilerplate code was added to support the specification of how many decimal digits that should be used for the
calculations,   as well how many decimal digits   (past the decimal point)   should be displayed.

A sequence that seems to converge to a wrong limit

<lang rexx>/*REXX pgm (pathological FP problem): 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.*/ if show== | show=="," then show= 20 /* " " " " " " */ numeric digits digs /*have REXX use "digs" decimal digits. */

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

fin=word(#, words(#) ) /*determine the last (largest) idx num.*/ w=length(fin) /*determine the length of largest index*/ v.1= 2 /*the value of the first V element. */ v.2=-4 /* " " " " second " " */ 

      do n=3  to fin;  nm1=n-1;      nm2=n-2    /*compute some values of the V elements*/
      v.n=111 - 1130/v.nm1 + 3000/(v.nm1*v.nm2) /*   "      a  value  of  a  " element.*/
      if wordpos(n, #)\==0  then say  'v.'left(n, w)      "="       format(v.n, , show)
      end   /*n*/                               /*display SHOW digits past the dec. pt.*/
                                                /*stick a fork in it,  we're all done. */</lang>

output &when using the default inputs: 

v.4   = 9.37837837837837837838
v.5   = 7.80115273775216138329
v.6   = 7.15441448097524935353
v.7   = 6.80678473692363298394
v.8   = 6.59263276870443839274
v.9   = 6.44946593379028796887
v.20  = 6.04355211018926886778
v.30  = 6.00678609303120575853
v.50  = 6.00017584662718718895
v.100 = 6.00000001931947792910

The Chaotic Bank Society

<lang rexx>/*REXX pgm (pathological FP problem): the chaotic bank society offering a new investment*/ parse arg digs show . /*obtain optional arguments from the CL*/ if digs== | digs=="," then digs=150 /*Not specified? Then use the default.*/ if show== | show=="," then show= 20 /* " " " " " " */ numeric digits digs /*have REXX use "digs" decimal digits. */ $=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138 $=$ - 1 /*and subtract one 'cause that's that. */ 

                                                /* [↑]  value of newly opened account. */
      do n=1 for 25                             /*compute the value of the account/year*/
      $=$*n - 1                                 /*   "     "    "    "  "  account now.*/
      end   /*n*/
                                                /*display SHOW digits past the dec. pt.*/

say 'Balance after 25 years:' format($, , show) /*display balance after 1/4 century. */ 

                                                /*stick a fork in it,  we're all done. */</lang>

output &when using the default inputs: 

Balance after 25 years: 0.03993872967323020890

Siegfried Rump's example

<lang rexx>/*REXX pgm (pathological FP problem): the Siegfried Rump's example (problem dated 1988).*/ parse arg digs show . /*obtain optional arguments from the CL*/ if digs== | digs=="," then digs=150 /*Not specified? Then use the default.*/ if show== | show=="," then show= 20 /* " " " " " " */ numeric digits digs /*have REXX use "digs" decimal digits. */ a= 77617.0 /*initialize A to it's defined value.*/ b= 33096.0 /* " B " " " " */ 

                                                /*display SHOW digits past the dec. pt.*/

say 'f(a,b)=' format( f(a,b), , show) /*display result from the F function.*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ f: procedure; parse arg a,b; a2=a**2; b2=b**2; b4=b2**2; b6=b4*b2; b8=b4**2 

    return  333.75*b6   +   a2*(11*a2*b2  - b6  - 121*b4  - 2)   +   5.5*b8   +   a/(2*b)</lang>

output &when using the default inputs: 

f(a,b)= -0.82739605994682136814

TI-83 BASIC

Use the SEQ mode to enter the arithmetic progression. Note the way to set 

 u(1)=2
 u(2)=-4

<lang ti83b> 

  nMin=1
  u(n)=111-1130/u(n-1) + 3000/(u(n-1)*u(n-2))
  u(nMin)={-4;2}

</lang> The result is the wrong limit:

Output:


u(20)  : 100.055202
u(30)  : 100
u(50)  : 100
u(100) : 100


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. <lang zkl>Series:=Walker(fcn(vs){ // just keep appending new values to a list 

  vs.append(111.0 - 1130.0/vs[-1] + 3000.0/(vs[-1]*vs[-2])) }.fp(List(2,-4)));

series:=Series.drop(100).value;</lang> We'll use the convenient formula given in the referenced paper to create a fraction with big ints <lang zkl>var BN=Import("zklBigNum"), ten2n=BN(10).pow(64);

fcn u(n){ // use formula to create a fraction of big ints 

  const B=-3, Y=4;
  N:=BN(6).pow(n+1)*B + BN(5).pow(n+1)*Y;
  D:=BN(6).pow(n)*B   + BN(5).pow(n)*Y;
  tostr(N*ten2n/D,64,32)

}

fcn tostr(bn,m,r){ // convert big int (*10^m) to float string with len r remainder, flakey 

  str,d:=bn.toString(), str.len()-m;
  if(d<0) String(".","0"*-d,str[0,r]);
  else    String(str[0,d],".",str[d,r]);

}

println("1st: Convergent series"); 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();

}</lang>

Output:

Note that, at n=100, we still have diverged (at the 15th place) from the Perl6 solution and 12th place from the J solution. 

1st: Convergent series
n =  3;  18.50000000000000000000  18.50000000000000000000000000000000
n =  4;   9.37837837837837895449  9.37837837837837837837837837837837
n =  5;   7.80115273775216877539  7.80115273775216138328530259365994
n =  6;   7.15441448097533339023  7.15441448097524935352789065386036
n =  7;   6.80678473692481134094  6.80678473692363298394175659627200
n =  8;   6.59263276872179204702  6.59263276870443839274200277636599
n = 20;  98.34950312216535905918  6.04355211018926886777747736409754
n = 30;  99.99999999999893418590  6.00678609303120575853055404795323
n = 50; 100.00000000000000000000  6.00017584662718718894561402074719
n =100; 100.00000000000000000000  6.00000001931947792910408680340358

Chaotic banking society is just nasty so we use a five hundred digit e (the e:= text is one long line). <lang zkl>println("\n2nd: Chaotic banking society"); e:="271828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723069697720931014169283681902551510865746377211125238978442505695369677078544996996794686445490598793163688923009879312"; var en=(e.len()-1), tenEN=BN(10).pow(en); years,balance:=25, BN(e).sub(tenEN); // in place math 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));</lang>

Output:
2nd: Chaotic banking society
After year 25, you will have $.039 in your account.

For Rump's example, multiple the formula by 10ⁿ so we can use integer math. <lang zkl>fcn rump(a,b){ b=BN(b); 

  b2,b4,b6,b8:=b.pow(2),b.pow(4),b.pow(6),b.pow(8);
  a2:=BN(a).pow(2);
  r:=( b6*33375 + a2*(a2*b2*11 - b6 - b4*121 - 2)*100 + b8*550 )*ten2n;
  r+=BN(a)*ten2n*100/(2*b);
  tostr(r,66,32)

} println("\n3rd: Rump's example: f(77617.0, 33096.0) = ",rump(77617,33096));</lang>

Output:
3rd: Rump's example: f(77617.0, 33096.0) = -.82739605994682136814116509547981
Retrieved from "https://rosettacode.org/wiki/Pathological_floating_point_problems?oldid=264542"