XIncludeFile "isprime.pb"

Structure PrimeHelper
  inc.i[8]
  index.i
EndStructure

Procedure.i firstPrimeFactor(n.q)
  If n = 1 : ProcedureReturn 1 : EndIf
  If n % 3 = 0 : ProcedureReturn 3 : EndIf
  If n % 5 = 0 : ProcedureReturn 5 : EndIf

  Define helper.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

  Define k.q = 7
  While k * k <= n
    If n % k = 0 : ProcedureReturn k : EndIf
    k + helper\inc[helper\index]
    helper\index = (helper\index + 1) % 8
  Wend
  ProcedureReturn n
EndProcedure

OpenConsole()
Define Dim blum.q(49)
Define.q bc = 0, i = 1
Define Dim counts.q(9)

Repeat
  Define p.q = firstPrimeFactor(i)
  If p % 4 = 3
    Define q.q = i / p
    If q <> p And q % 4 = 3 And isPrime(q)
      If bc < 50 : blum(bc) = i : EndIf
      counts(i % 10) + 1
      bc + 1

      If bc = 50
        PrintN("First 50 Blum integers:")
        For j = 0 To 49
          Print(" " + RSet(Str(blum(j)), 3))
          If (j + 1) % 10 = 0 : PrintN("") : EndIf
        Next
        PrintN("")
      ElseIf bc = 26828 Or bc % 100000 = 0
        PrintN("The " + RSet(Str(bc), 6) + "th Blum integer is: " + RSet(Str(i), 7))
        If bc = 400000
          PrintN(#CRLF$ + "% distribution of the first 400,000 Blum integers:")
          For j = 1 To 9 Step 2
            If j <> 5
              PrintN(RSet(StrF(counts(j)/4000, 3), 5) + "% end in " + Str(j))
            EndIf
          Next
          Break
        EndIf
      EndIf
    EndIf
  EndIf
  If i % 5 = 3
    i + 4
  Else
    i + 2
  EndIf
ForEver

PrintN(#CRLF$ + "Press ENTER to exit"): Input()