89 lines
3.0 KiB
Plaintext
89 lines
3.0 KiB
Plaintext
#
|
|
# 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
|