with javascript_semantics
constant LIMIT = 1e7, N = floor((floor(LIMIT/3)-1)/4)+1

function Sieve4n_3_Primes()
    sequence sieve = repeat(0,N), p4n3 = {}
    for idx=1 to N do
        if sieve[idx]=0 then
            integer n = idx*4-1
            p4n3 &= n
            if idx+n>N then
                // collect the rest
                for j=idx+1 to N do
                    if sieve[j]=0 then
                        p4n3 &= 4*j-1
                    end if
                end for
                exit
            end if
            for j=idx+n to N by n do
                sieve[j] = 1
            end for
        end if
    end for
    return p4n3
end function

sequence p4n3 = Sieve4n_3_Primes(),
         BlumField = repeat(false,LIMIT),
         Blum50 = {}, counts = repeat(0,10)

for idx,n in p4n3 do
    for bj in p4n3 from idx+1 do
        atom k = n*bj
        if k>LIMIT then exit end if
        BlumField[k] = true
    end for
end for
integer count = 0
for n,k in BlumField do
    if k then
        if count<50 then Blum50 &= n end if
        counts[remainder(n,10)] += 1
        count += 1
        if count=50 then
            printf(1,"First 50 Blum integers:\n%s\n",{join_by(Blum50,1,10," ",fmt:="%3d")})
        elsif count=26828 or remainder(count,1e5)=0 then
            printf(1,"The %,7d%s Blum integer is: %,9d\n", {count,ord(count),n})
            if count=4e5 then exit end if
        end if
    end if
end for
printf(1,"\nPercentage distribution of the first 400,000 Blum integers:\n")
for i,n in counts do
    if n then
        printf(1,"  %6.3f%% end in %d\n", {n/4000, i})
    end if
end for