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.

This task's purpose is to 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.



Task 1

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



Task 2

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



Task 3, extra credit

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

  f(a,b) = 333.75b6 + a2( 11a2b2 - b6 - 121b4 - 2 ) + 5.5b8 + 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;
   Floating-Point Arithmetic   Section 1.3.2 Difficult problems.


 
program Pathological_floating_point_problems;
 
{$APPTYPE CONSOLE}
 
uses
  System.SysUtils,
  Velthuis.BigIntegers,
  Velthuis.BigRationals,
  Velthuis.BigDecimals;
 
type
  TBalance = record
    e, d: BigInteger;
  end;
 
procedure Sequence(Precision: Integer);
var
  v, v1, c111, c1130, c3000, t2, t3: BigRational;
  n: integer;
begin
  BigDecimal.DefaultPrecision := Precision;
  v1 := 2;
  v := -4;
  n := 2;
  c111 := BigRational(111);
  c1130 := BigRational(1130);
  c3000 := BigRational(3000);
 
  var r :=
    function(): BigRational
    begin
      t3 := v * v1;
      t3 := BigRational.Create(BigDecimal.Divide(c3000, t3));
      t2 := BigRational.Create(BigDecimal.Divide(c1130, v));
      result := c111 - t2;
      result := result + t3;
    end;
 
  writeln('  n  sequence value');
 
  for var x in [3, 4, 5, 6, 7, 8, 20, 30, 50, 100] do
  begin
    while n < x do
    begin
      var tmp := BigRational.Create(v);
      v := BigRational.Create(r());
      v1 := BigRational.Create(tmp);
      inc(n);
    end;
 
    var f := double(v);
    writeln(format('%3d %19.16f', [n, f]));
  end;
end;
 
procedure Bank();
var
  balance: TBalance;
  m, one, ef, df: BigInteger;
  e, b: BigDecimal;
begin
  balance.e := 1;
  balance.d := -1;
  one := BigInteger.One;
 
  for var y := 1 to 25 do
  begin
    m := y;
    balance.e := m * balance.e;
    balance.d := m * balance.d;
    balance.d := balance.d - one;
  end;
 
  e := BigDecimal.Create('2.71828182845904523536028747135');
 
  ef := balance.e;
  df := balance.d;
 
  b := e * ef;
  b := b + df;
 
  writeln(format('Bank balance after 25 years: $%.2f', [b.AsDouble]));
end;
 
function f(a, b: Double): Double;
var
  a1, a2, b1, b2, b4, b6, b8, two, t1, t21, t2, t3, t4, t23: BigDecimal;
begin
  var fp :=
    function(x: double): BigDecimal
    begin
      Result := BigDecimal.Create(x);
    end;
 
  a1 := fp(a);
  b1 := fp(b);
  a2 := a1 * a1;
  b2 := b1 * b1;
  b4 := b2 * b2;
  b6 := b2 * b4;
  b8 := b4 * b4;
 
  two := fp(2);
  t1 := fp(333.75);
  t1 := t1 * b6;
  t21 := fp(11);
  t21 := t21 * a2 * b2;
  t23 := fp(121);
  t23 := t23 * b4;
 
  t2 := t21 - b6;
  t2 := (t2 - t23) - two;
  t2 := a2 * t2;
 
  t3 := fp(5.5);
  t3 := t3 * b8;
 
  t4 := two * b1;
  t4 := a1 / t4;
 
  var s := t1 + t2;
  s := s + t3;
  s := s + t4;
  result := s.AsDouble;
end;
 
procedure Rump();
var
  a, b: Double;
begin
  a := 77617;
  b := 33096;
  writeln(format('Rump f(%g, %g): %g', [a, b, f(a, b)]));
end;
 
{ TBigRationalHelper }
 
 
begin
  for var Precision in [100, 200] do
  begin
    writeln(#10'Precision: ', Precision);
    sequence(Precision);
    bank();
    rump();
  end;
 
  readln;
end.