' version 27-10-2016
' compile with: fbc -s console

#Define max 1000    ' total number of Fibonacci numbers
#Define max_sieve  15485863  ' should give 1,000,000

#Include Once "gmp.bi"   ' uses the GMP libary

Dim As ZString Ptr z_str
Dim As ULong n, d
ReDim As ULong digit(1 To 9)
Dim As Double expect, found

Dim As mpz_ptr fib1, fib2
fib1 = Allocate(Len(__mpz_struct)) : Mpz_init_set_ui(fib1, 0)
fib2 = Allocate(Len(__mpz_struct)) : Mpz_init_set_ui(fib2, 1)

digit(1) = 1  ' fib2
For n = 2 To max
    Swap fib1, fib2                   ' fib1 = 1, fib2 = 0
    mpz_add(fib2, fib1, fib2)         ' fib1 = 1, fib2 = 1 (fib1 + fib2)
    z_str = mpz_get_str(0, 10, fib2)
    d = Val(Left(*z_str, 1))          ' strip the 1 digit on the left off
    digit(d) = digit(d) +1
Next

mpz_clear(fib1) : DeAllocate(fib1)
mpz_clear(fib2) : DeAllocate(fib2)

Print
Print "First 1000 Fibonacci numbers"
Print "nr:  total     found   expected  difference"

For d = 1 To 9
    n = digit(d)
    found = n / 10
    expect = (Log(1 + 1 / d) / Log(10)) * 100
    Print Using " ##  #####  ###.## %   ###.## %    ##.### %"; _
                            d; n ; found; expect; expect - found
Next


ReDim digit(1 To 9)
ReDim As UByte sieve(max_sieve)

'For d = 4 To max_sieve Step 2
'    sieve(d) = 1
'Next
Print : Print "start sieve"
For d = 3 To sqr(max_sieve)
    If sieve(d) = 0 Then
        For n = d * d To max_sieve Step d * 2
            sieve(n) = 1
        Next
    End If
Next

digit(2) = 1 ' 2

Print "start collecting first digits"
For n = 3 To max_sieve Step 2
    If sieve(n) = 0 Then
        d = Val(Left(Trim(Str(n)), 1))
        digit(d) = digit(d) +1
    End If
Next

Dim As ulong total
For n = 1 To 9
    total = total + digit(n)
Next

Print
Print "First";total; " primes"
Print "nr:     total     found   expected   difference"

For d = 1 To 9
    n = digit(d)
    found = n / total * 100
    expect = (Log(1 + 1 / d) / Log(10)) * 100
    Print Using " ##  ########  ###.## %   ###.## %    ###.### %"; _
                                d; n ; found; expect; expect - found
Next


' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End