#include "isprime.bas"

Type PrimeHelper
    inc(7) As Integer
    idx As Integer
End Type

Function initPrimeHelper() As PrimeHelper
    Dim helper As PrimeHelper
    helper.inc(0) = 4 : helper.inc(1) = 2 : helper.inc(2) = 4
    helper.inc(3) = 2 : helper.inc(4) = 4 : helper.inc(5) = 6
    helper.inc(6) = 2 : helper.inc(7) = 6
    helper.idx = 0
    Return helper
End Function

Function firstPrimeFactor(n As Integer) As Integer
    If n = 1 Then Return 1
    If n Mod 3 = 0 Then Return 3
    If n Mod 5 = 0 Then Return 5

    Dim helper As PrimeHelper = initPrimeHelper()
    Dim k As Integer = 7

    While k * k <= n
        If n Mod k = 0 Then Return k
        k += helper.inc(helper.idx)
        helper.idx = (helper.idx + 1) Mod 8
    Wend

    Return n
End Function

Sub main()
    Dim As Integer blum(49), counts(9)
    Dim As Integer bc = 0, i = 1, p, q
    Dim As Integer j

    Dim As Double t0 = Timer
    Do
        p = firstPrimeFactor(i)
        If p Mod 4 = 3 Then
            q = i \ p
            If q <> p Andalso q Mod 4 = 3 Andalso isPrime(q) Then
                If bc < 50 Then blum(bc) = i
                counts(i Mod 10) += 1
                bc += 1

                If bc = 50 Then
                    Print "First 50 Blum integers:"
                    For j = 0 To 49
                        Print Using "####"; blum(j);
                        If (j + 1) Mod 10 = 0 Then Print
                    Next
                    Print
                Elseif bc = 26828 Orelse bc Mod 100000 = 0 Then
                    Print Using "The ###,###th Blum integer is: #,###,###"; bc; i

                    If bc = 400000 Then
                        Print !"\n% distribution of the first 400,000 Blum integers:"
                        For j = 1 To 9 Step 2
                            If j <> 5 Then Print Using "  ##.###% end in #"; (counts(j)/4000); j
                        Next
                        Exit Do
                    End If
                End If
            End If
        End If
        i += Iif(i Mod 5 = 3, 4, 2)
    Loop
    Print Chr(10); Timer - t0; " sec."
End Sub

main()

Sleep