41 lines
1.1 KiB
Text
41 lines
1.1 KiB
Text
GET "libhdr"
|
|
|
|
MANIFEST {
|
|
fib_sz = 1_000_000_000 // keep 9 significant digits for computing F(n)
|
|
}
|
|
|
|
LET msd(n) = VALOF {
|
|
UNTIL n < 10 DO
|
|
n /:= 10
|
|
RESULTIS n
|
|
}
|
|
|
|
LET fibonacci(n, tally) BE {
|
|
LET a, b, c, e = 0, 1, ?, 0
|
|
FOR i = 1 TO n {
|
|
TEST e = 0
|
|
THEN tally[msd(b)] +:= 1
|
|
ELSE tally[b / (fib_sz / 10)] +:= 1
|
|
c := a + b
|
|
IF c > fib_sz {
|
|
a := a / 10 - (a MOD 10 >= 5) // subtract, since condition evalutes to
|
|
b := b / 10 - (b MOD 10 >= 5) // eithr 0 or -1.
|
|
c := a + b
|
|
e +:= 1 // keep track of exponent, just 'cuz
|
|
}
|
|
a, b := b, c
|
|
}
|
|
}
|
|
|
|
LET start() = VALOF {
|
|
// expected value of benford: log10(d + 1/d)
|
|
LET expected = TABLE 0, 301, 176, 125, 97, 79, 67, 58, 51, 46
|
|
LET actual = VEC 9
|
|
FOR i = 0 TO 9 DO actual!i := 0
|
|
fibonacci(1000, actual)
|
|
writes("*nLeading digital distribution of the first 1,000 Fibonacci numbers*n")
|
|
writes("Digit*tActual*tExpected*n")
|
|
FOR i = 1 TO 9
|
|
writef("%i *t%0.3d *t%0.3d *n", i, actual!i, expected!i)
|
|
RESULTIS 0
|
|
}
|