include xpllib;         \for Print

int Prime1;

func Semiprime(N);      \Returns 'true' if odd N >= 3 is semiprime
int  N, D, C;
[C:= 0;  D:= 3;
while D*D <= N do
    [while rem(N/D) = 0 do
        [if C = 2 then return false;
        C:= C+1;
        N:= N/D;
        ];
    D:= D+2;
    ];
Prime1:= N;
return C = 1;
];

int N, C, I, Goal, Prime2;
int FD, DC(10);         \final digit and digit counters
[Text(0, "First 50 Blum integers:^m^j");
N:= 3;  C:= 0;  Goal:= 100_000;
for I:= 0 to 9 do DC(I):= 0;
loop    [if Semiprime(N) then
            [if rem(Prime1/4) = 3 then
                [Prime2:= N/Prime1;
                if rem(Prime2/4) = 3 and Prime2 # Prime1 then
                    [C:= C+1;
                    if C <= 50 then
                        [Print("%5d", N);
                        if rem(C/10) = 0 then Print("\n");
                        ];
                    if C = 26_828 then
                        Print("\nThe  26,828th Blum integer is: %7,d\n", N);
                    if C = Goal then
                        [Print("The %6,dth Blum integer is: %7,d\n", Goal, N);
                        if Goal = 400_000 then quit;
                        Goal:= Goal + 100_000;
                        ];
                    FD:= rem(N/10);
                    DC(FD):= DC(FD)+1;
                    ];
                ];
        ];
    N:= N+2;
    ];
Print("\n% distribution of the first 400,000 Blum integers:\n");
for I:= 0 to 9 do
    if DC(I) > 0 then
        Print("%4.3f\% end in %d\n", float(DC(I)) / 4_000., I);
]