##-
# Author: Brian Tiffin
# Dedicated to the public domain
#
# Date started: December, 2016
# Modified: 2016-12-24/02:50-0500
##+
#
# 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, e - $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 annual 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