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.


#
# Pathological floating point problems
#
procedure main()
    sequence()
    chaotic()
end
 
#
# First task, sequence convergence
#
link printf
procedure sequence()
     local l := [2, -4]
     local iters := [3, 4, 5, 6, 7, 8, 20, 30, 50, 100, 200]
     local i, j, k  
     local n := 1
 
     write("Sequence convergence")
     # Demonstrate the convergence problem with various precision values
     every k := (100 | 300) do {
         n := 10^k
         write("\n", k, " digits of intermediate precision")
 
         # numbers are scaled up using large integer powers of 10
         every i := !iters do {
             l := [2 * n, -4 * n]
             printf("i: %3d", i)
 
             every j := 3 to i do {
                 # build out a list of intermediate passes
                 # order of scaling operations matters
                 put(l, 111 * n - (1130 * n * n / l[j - 1]) +
                        (3000 * n * n * n / (l[j - 1] * l[j - 2])))
             }
             # down scale the result to a real
             # some precision may be lost in the final display
             printf(" %20.16r\n", l[i] * 1.0 / n)
         }
     }
end
 
#
# Task 2, chaotic bank of Euler
#
procedure chaotic()
    local euler, e, scale, show, y, d
 
    write("\nChaotic Banking Society of Euler")
    # format the number for listing, string form, way overboard on digits
    euler :=
"2718281828459045235360287471352662497757247093699959574966967627724076630353_
  547594571382178525166427427466391932003059921817413596629043572900334295260_
  595630738132328627943490763233829880753195251019011573834187930702154089149_
  934884167509244761460668082264800168477411853742345442437107539077744992069_
  551702761838606261331384583000752044933826560297606737113200709328709127443_
  747047230696977209310141692836819025515108657463772111252389784425056953696_
  770785449969967946864454905987931636889230098793127736178215424999229576351_
  482208269895193668033182528869398496465105820939239829488793320362509443117_
  301238197068416140397019837679320683282376464804295311802328782509819455815_
  301756717361332069811250996181881593041690351598888519345807273866738589422_
  879228499892086805825749279610484198444363463244968487560233624827041978623_
  209002160990235304369941849146314093431738143640546253152096183690888707016_
  768396424378140592714563549061303107208510383750510115747704171898610687396_
  9655212671546889570350354"
 
    # precise math with long integers, string form just for pretty listing
    e := integer(euler)
 
    # 1000 digits after the decimal for scaling intermediates and service fee
    scale := 10^1000
 
    # initial deposit - $1
    d := e - scale
 
    # show balance with 16 digits
    show := 10^16
    write("Starting balance:       $", d * show / scale * 1.0 / show, "...")
 
    # wait 25 years, with only a trivial $1 service fee
    every y := 1 to 25 do {
        d := d * y - scale
    }
 
    # show final balance with 4 digits after the decimal (truncation)
    show := 10^4
    write("Balance after ", y, " years: $", d * show / scale * 1.0 / show)
end