#include "isprime.bas"
uses console

dim inc(7) as integer
inc(0) = 4: inc(1) = 2: inc(2) = 4
inc(3) = 2: inc(4) = 4: inc(5) = 6
inc(6) = 2: inc(7) = 6

function firstPrimeFactor(n as long) as long
    if n = 1 then return 1
    if n mod 3 = 0 then return 3
    if n mod 5 = 0 then return 5

    long k = 7
    int idx = 0

    while k * k <= n
        if mod(n, k) = 0 then return k
        k += inc(idx)
        idx = mod((idx + 1), 8)
    wend
    return n
end function

sub main()
    dim as long blum(49), counts(9)
    long bc = 0, i = 1, pct, p, q
    int j
    string s, t

    do
        p = firstPrimeFactor(i)
        if p mod 4 = 3 then
            q = i \ p
            if q <> p and q mod 4 = 3 and isPrime(q) then
                if bc < 50 then blum(bc) = i
                counts(i mod 10) = counts(i mod 10) + 1
                bc += 1

                if bc = 50 then
                    printl "First 50 Blum integers:"
                    for j = 0 to 49
                        s = str(blum(j))
                        while len(s) < 4: s = " " + s: wend
                        print s;
                        if mod((j + 1), 10) = 0 then printl
                    next
                    printl
                elseif bc = 26828 or bc mod 100000 = 0 then
                    s = str(bc)
                    while len(s) < 6: s = " " + s: wend
                    t = str(i)
                    while len(t) < 7: t = " " + t: wend
                    printl "The " + s + "th Blum integer is: " + t
                    if bc = 400000 then
                        printl cr "% distribution of the first 400,000 Blum integers:"

                        for j = 1 to 9 step 2
                            if j <> 5 then
                                pct = counts(j)/4000
                                s = str(pct)
                                while len(s) < 5: s = " " + s: wend
                                printl s + "% end in " + str(j)
                            end if
                        next
                        end
                    end if
                end if
            end if
        end if
        if mod(i, 5) = 3 then i += 4 else i += 2
    end do
end sub

main()

printl cr "Enter ..."
waitkey