53 lines
1.5 KiB
Plaintext
53 lines
1.5 KiB
Plaintext
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);
|
|
]
|