54 lines
1.6 KiB
Text
54 lines
1.6 KiB
Text
char Sieve;
|
|
proc MakeSieve(Size); \Make prime number sieve
|
|
int Size, Prime, I, K;
|
|
[Size:= Size/2; \ignore even numbers
|
|
Sieve:= MAlloc(Size+1); \(XPL0's heap only provides 64 MB)
|
|
for I:= 0 to Size do \set Sieve flags all true
|
|
Sieve(I):= true;
|
|
for I:= 0 to Size do
|
|
if Sieve(I) then \found a prime, which is equal to
|
|
[Prime:= I + I + 3; \ twice the index + 3
|
|
K:= I + Prime; \first multiple to strike off
|
|
while K <= Size do \strike off all multiples
|
|
[Sieve(K):= false;
|
|
K:= K + Prime;
|
|
];
|
|
];
|
|
];
|
|
|
|
func GetSig(N); \Return signature of N
|
|
\A "signature" is the count of each digit in N packed into a 32-bit word
|
|
int N, Sig;
|
|
[Sig:= 0;
|
|
repeat N:= N/10;
|
|
Sig:= Sig + 1<<(rem(0)*3);
|
|
until N = 0;
|
|
return Sig;
|
|
];
|
|
|
|
def Limit = 1_000_000_000;
|
|
int Cnt, N, N0, N1, Sig, Sig0, Sig1;
|
|
[MakeSieve(Limit);
|
|
Text(0, "Smallest members of first 25 Ormiston triples:^m^j");
|
|
Cnt:= 0; N0:= 0; N1:= 0; Sig0:= 0; Sig1:= 0; N:= 3;
|
|
Format(10, 0);
|
|
loop [if Sieve(N>>1-1) then \is prime
|
|
[Sig:= GetSig(N);
|
|
if Sig = Sig1 and Sig = Sig0 then
|
|
[Cnt:= Cnt+1;
|
|
if Cnt <= 25 then
|
|
[RlOut(0, float(N0));
|
|
if rem(Cnt/5) = 0 then CrLf(0);
|
|
];
|
|
];
|
|
Sig0:= Sig1; Sig1:= Sig;
|
|
N0:= N1; N1:= N;
|
|
];
|
|
if N >= Limit then
|
|
[IntOut(0, Cnt);
|
|
Text(0, " Ormiston triples before one billion.^m^j");
|
|
quit;
|
|
];
|
|
N:= N+2;
|
|
];
|
|
]
|