Ormiston triples
An Ormiston triple is three consecutive prime numbers which are anagrams, i.e. contain the same decimal digits but in a different order.
Ormiston triples is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- Definition
- Example
The three consecutive primes (11117123, 11117213, 11117321) are an Ormiston triple.
- Task
- Find and show the smallest member of the first 25 Ormiston triples.
- Find and show the count of Ormiston triples up to one billion.
- Stretch
- Find and show the count of Ormiston triples up to ten billion.
- Reference
- Related task
C
Run time about 83 seconds.
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stdint.h>
#include <locale.h>
bool *sieve(uint64_t limit) {
uint64_t i, p;
limit++;
// True denotes composite, false denotes prime.
bool *c = calloc(limit, sizeof(bool)); // all false by default
c[0] = true;
c[1] = true;
for (i = 4; i < limit; i += 2) c[i] = true;
p = 3; // Start from 3.
while (true) {
uint64_t p2 = p * p;
if (p2 >= limit) break;
for (i = p2; i < limit; i += 2 * p) c[i] = true;
while (true) {
p += 2;
if (!c[p]) break;
}
}
return c;
}
typedef struct {
char digs[20];
int count;
} digits;
digits getDigits(uint64_t n) {
if (n == 0) return (digits){ {0}, 1 };
digits d;
d.count = 0;
while (n > 0) {
d.digs[d.count++] = n % 10;
n = n / 10;
}
return d; // note digits are in reverse order
}
int main() {
const uint64_t limit = 10000000000;
uint64_t i, j, pc = 0, p1, p2, p3, key1, key2, key3;
int k, count, count2;
digits d;
bool *c = sieve(limit);
for (i = 0; i < limit; ++i) {
if (!c[i]) ++pc;
}
uint64_t *primes = (uint64_t *)malloc(pc * sizeof(uint64_t));
for (i = 0, j = 0; i < limit; ++i) {
if (!c[i]) primes[j++] = i;
}
free(c);
uint64_t orm25[25];
int counts[2];
j = limit/10;
for (i = 0; i < pc-2; ++i) {
p1 = primes[i];
p2 = primes[i+1];
p3 = primes[i+2];
if ((p2 - p1) % 18 || (p3 - p2) % 18) continue;
key1 = 1;
d = getDigits(p1);
for (k = 0; k < d.count; ++k) key1 *= primes[d.digs[k]];
key2 = 1;
d = getDigits(p2);
for (k = 0; k < d.count; ++k) key2 *= primes[d.digs[k]];
if (key1 != key2) continue;
key3 = 1;
d = getDigits(p3);
for (k = 0; k < d.count; ++k) key3 *= primes[d.digs[k]];
if (key2 == key3) {
if (count < 25) orm25[count] = p1;
if (p1 >= j) {
counts[count2++] = count;
j *= 10;
}
++count;
}
}
counts[count2] = count;
printf("Smallest members of first 25 Ormiston triples:\n");
setlocale(LC_NUMERIC, "");
for (i = 0; i < 25; ++i) {
printf("%'10ld ", orm25[i]);
if (!((i+1) % 5)) printf("\n");
}
printf("\n");
j = limit/10;
for (i = 0; i < 2; ++i) {
printf("%'d Ormiston triples before %'ld\n", counts[i], j);
j *= 10;
printf("\n");
}
free(primes);
return 0;
}
- Output:
Smallest members of first 25 Ormiston triples: 11,117,123 12,980,783 14,964,017 32,638,213 32,964,341 33,539,783 35,868,013 44,058,013 46,103,237 48,015,013 50,324,237 52,402,783 58,005,239 60,601,237 61,395,239 74,699,789 76,012,879 78,163,123 80,905,879 81,966,341 82,324,237 82,523,017 83,279,783 86,050,781 92,514,341 368 Ormiston triples before 1,000,000,000 4,925 Ormiston triples before 10,000,000,000
C++
#include <array>
#include <iostream>
#include <primesieve.hpp>
class ormiston_triple_generator {
public:
ormiston_triple_generator() {
for (int i = 0; i < 2; ++i) {
primes_[i] = pi_.next_prime();
digits_[i] = get_digits(primes_[i]);
}
}
std::array<uint64_t, 3> next_triple() {
for (;;) {
uint64_t prime = pi_.next_prime();
auto digits = get_digits(prime);
bool is_triple = digits == digits_[0] && digits == digits_[1];
uint64_t prime0 = primes_[0];
primes_[0] = primes_[1];
primes_[1] = prime;
digits_[0] = digits_[1];
digits_[1] = digits;
if (is_triple)
return {prime0, primes_[0], primes_[1]};
}
}
private:
static std::array<int, 10> get_digits(uint64_t n) {
std::array<int, 10> result = {};
for (; n > 0; n /= 10)
++result[n % 10];
return result;
}
primesieve::iterator pi_;
std::array<uint64_t, 2> primes_;
std::array<std::array<int, 10>, 2> digits_;
};
int main() {
ormiston_triple_generator generator;
int count = 0;
std::cout << "Smallest members of first 25 Ormiston triples:\n";
for (; count < 25; ++count) {
auto primes = generator.next_triple();
std::cout << primes[0] << ((count + 1) % 5 == 0 ? '\n' : ' ');
}
std::cout << '\n';
for (uint64_t limit = 1000000000; limit <= 10000000000; ++count) {
auto primes = generator.next_triple();
if (primes[2] > limit) {
std::cout << "Number of Ormiston triples < " << limit << ": "
<< count << '\n';
limit *= 10;
}
}
}
- Output:
Smallest members of first 25 Ormiston triples: 11117123 12980783 14964017 32638213 32964341 33539783 35868013 44058013 46103237 48015013 50324237 52402783 58005239 60601237 61395239 74699789 76012879 78163123 80905879 81966341 82324237 82523017 83279783 86050781 92514341 Number of Ormiston triples < 1000000000: 368 Number of Ormiston triples < 10000000000: 4925
F#
This task uses Ormiston_pairs#F#
// Ormiston triples. Nigel Galloway: February 3rd., 2023
let oTriples n=n|>Seq.pairwise|>Seq.filter(fun((n,i),(g,l))->i=g)|>Seq.map(fun((n,i),(g,l))->(n,g,l))
primes32()|>oPairs|>oTriples|>Seq.take 25|>Seq.iter(fun(n,_,_)->printf "%d " n); printfn ""
printfn $"<100 million: %d{primes32()|>Seq.takeWhile((>)100000000)|>oPairs|>oTriples|>Seq.length}"
printfn $"<1 billion: %d{primes32()|>Seq.takeWhile((>)1000000000)|>oPairs|>oTriples|>Seq.length}"
- Output:
11117123 12980783 14964017 32638213 32964341 33539783 35868013 44058013 46103237 48015013 50324237 52402783 58005239 60601237 61395239 74699789 76012879 78163123 80905879 81966341 82324237 82523017 83279783 86050781 92514341 <100 million: 25 <1 billion: 368
Go
This runs in about 54 seconds on my Core i7 machine.
package main
import (
"fmt"
"rcu"
)
func main() {
const limit = 1e10
primes := rcu.Primes(limit)
var orm25 []int
j := int(1e9)
count := 0
var counts []int
for i := 0; i < len(primes)-2; i++ {
p1 := primes[i]
p2 := primes[i+1]
p3 := primes[i+2]
if (p2-p1)%18 != 0 || (p3-p2)%18 != 0 {
continue
}
key1 := 1
for _, dig := range rcu.Digits(p1, 10) {
key1 *= primes[dig]
}
key2 := 1
for _, dig := range rcu.Digits(p2, 10) {
key2 *= primes[dig]
}
if key1 != key2 {
continue
}
key3 := 1
for _, dig := range rcu.Digits(p3, 10) {
key3 *= primes[dig]
}
if key2 == key3 {
if count < 25 {
orm25 = append(orm25, p1)
}
if p1 >= j {
counts = append(counts, count)
j *= 10
}
count++
}
}
counts = append(counts, count)
fmt.Println("Smallest members of first 25 Ormiston triples:")
for i := 0; i < 25; i++ {
fmt.Printf("%8v ", orm25[i])
if (i+1)%5 == 0 {
fmt.Println()
}
}
fmt.Println()
j = int(1e9)
for i := 0; i < len(counts); i++ {
fmt.Printf("%s Ormiston triples before %s\n", rcu.Commatize(counts[i]), rcu.Commatize(j))
j *= 10
fmt.Println()
}
}
- Output:
Smallest members of first 25 Ormiston triples: 11117123 12980783 14964017 32638213 32964341 33539783 35868013 44058013 46103237 48015013 50324237 52402783 58005239 60601237 61395239 74699789 76012879 78163123 80905879 81966341 82324237 82523017 83279783 86050781 92514341 368 Ormiston triples before 1,000,000,000 4,925 Ormiston triples before 10,000,000,000
Haskell
import Data.Numbers.Primes (primes)
import Data.List (sort)
--------------------- ORMISTON TRIPLES -------------------
ormistons :: [(Integer, Integer, Integer)]
ormistons =
concat $ zipWith3
(\(dx, x) (dy, y) (dz, z)
-> [(x, y, z) | dx == dy && dx == dz])
primeDigits
(tail primeDigits)
(drop 2 primeDigits)
primeDigits :: [(Integer, Integer)]
primeDigits = ((,) =<< read . sort . show) <$> primes
--------------------------- TEST -------------------------
main :: IO ()
main = do
putStrLn "First 25 Ormistons:"
mapM_ print $ take 25 ormistons
putStrLn "\nCount of Ormistons up to 10^8:"
let limit = 10^8
print $ length $
takeWhile (\(_, _, c) -> c <= limit) ormistons
- Output:
First 25 Ormistons: (11117123,11117213,11117321) (12980783,12980837,12980873) (14964017,14964071,14964107) (32638213,32638231,32638321) (32964341,32964413,32964431) (33539783,33539837,33539873) (35868013,35868031,35868103) (44058013,44058031,44058103) (46103237,46103273,46103327) (48015013,48015031,48015103) (50324237,50324273,50324327) (52402783,52402837,52402873) (58005239,58005293,58005329) (60601237,60601273,60601327) (61395239,61395293,61395329) (74699789,74699879,74699897) (76012879,76012897,76012987) (78163123,78163213,78163231) (80905879,80905897,80905987) (81966341,81966413,81966431) (82324237,82324273,82324327) (82523017,82523071,82523107) (83279783,83279837,83279873) (86050781,86050817,86050871) (92514341,92514413,92514431) Count of Ormistons up to 10^8: 25
J
Taking the laziest approach here, we'll use the definition of isorm from the Ormiston pairs task:
omt=: (#~ isorm) _4 p: (#~ isorm) i.&.(p:inv) 1e9
#omt NB. number of ormiston triples less than a billion
368
5 5$omt NB. first prime of the first 25 triples.
11117123 12980783 14964017 32638213 32964341
33539783 35868013 44058013 46103237 48015013
50324237 52402783 58005239 60601237 61395239
74699789 76012879 78163123 80905879 81966341
82324237 82523017 83279783 86050781 92514341
Julia
Takes about 2 minutes, without the benefit of the mod 18 optimization.
using Primes
const dig = zeros(Int, 10)
const primes10B = primes(10_000_000_000)
const sort10B = map(n -> evalpoly(10, sort!(digits!(dig, n))), primes10B)
const ormiston_indices = filter(i -> sort10B[i] == sort10B[i + 1] == sort10B[i + 2],
firstindex(sort10B):lastindex(sort10B) - 2)
println("First 25 Ormiston triples:")
for i in 1:25
println(primes10B[ormiston_indices[i]:ormiston_indices[i]+2])
end
println("\nOrmiston triples before 1 billion: ",
count(t -> primes10B[t] < 1_000_000_000, ormiston_indices))
println("Ormiston triples before 10 billion: ", length(ormiston_indices))
- Output:
First 25 Ormiston triples: [11117123, 11117213, 11117321] [12980783, 12980837, 12980873] [14964017, 14964071, 14964107] [32638213, 32638231, 32638321] [32964341, 32964413, 32964431] [33539783, 33539837, 33539873] [35868013, 35868031, 35868103] [44058013, 44058031, 44058103] [46103237, 46103273, 46103327] [48015013, 48015031, 48015103] [50324237, 50324273, 50324327] [52402783, 52402837, 52402873] [58005239, 58005293, 58005329] [60601237, 60601273, 60601327] [61395239, 61395293, 61395329] [74699789, 74699879, 74699897] [76012879, 76012897, 76012987] [78163123, 78163213, 78163231] [80905879, 80905897, 80905987] [81966341, 81966413, 81966431] [82324237, 82324273, 82324327] [82523017, 82523071, 82523107] [83279783, 83279837, 83279873] [86050781, 86050817, 86050871] [92514341, 92514413, 92514431] Ormiston triples before 1 billion: 368 Ormiston triples before 10 billion: 4925
Pascal
Free Pascal
checking for used digits like Ormiston_triples#Go and removed mod 18 calculation by using an array of boolean [0..990] with every 18 marked as true.Not as effective as I expected .
program Ormiston;
{$IFDEF FPC}{$MODE DELPHI} {$OPTIMIZATION ON,ALL}{$ENDIF}
{$IFDEF WINDOWS}{$APPLICATION CONSOLE}{$ENDIF}
uses
sysutils,strUtils;
//********* segmented sieve of erathostenes *********
{segmented sieve of Erathostenes using only odd numbers}
{using presieved sieve of small primes, to reduce the most time consuming}
const
smlPrimes :array [0..10] of Byte = (2,3,5,7,11,13,17,19,23,29,31);
maxPreSievePrimeNum = 7;
maxPreSievePrime = 17;//smlPrimes[maxPreSievePrimeNum];
cSieveSize = 4*16384;//<= High(Word)+1 // Level I Data Cache
type
tSievePrim = record
svdeltaPrime:word;//diff between actual and new prime
svSivOfs:word; //Offset in sieve
svSivNum:LongWord;//1 shl (1+16+32) = 5.6e14
end;
tpSievePrim = ^tSievePrim;
var
//sieved with primes 3..maxPreSievePrime.here about 255255 Byte
{$ALIGN 32}
preSieve :array[0..3*5*7*11*13*17-1] of Byte;//must be > cSieveSize
{$ALIGN 32}
Sieve :array[0..cSieveSize-1] of Byte;
{$ALIGN 32}
//prime = FoundPrimesOffset + 2*FoundPrimes[0..FoundPrimesCnt]
FoundPrimes : array[0..12252] of Word;
{$ALIGN 32}
sievePrimes : array[0..1077863] of tSievePrim;
FoundPrimesOffset : Uint64;
FoundPrimesCnt,
FoundPrimesIdx,
FoundPrimesTotal,
SieveNum,
SieveMaxIdx,
preSieveOffset,
LastInsertedSievePrime :NativeUInt;
procedure CopyPreSieveInSieve; forward;
procedure CollectPrimes; forward;
procedure sieveOneSieve; forward;
procedure Init0Sieve; forward;
procedure SieveOneBlock; forward;
procedure preSieveInit;
var
i,pr,j,umf : NativeInt;
Begin
fillchar(preSieve[0],SizeOf(preSieve),#1);
i := 1;
pr := 3;// starts with pr = 3
umf := 1;
repeat
IF preSieve[i] =1 then
Begin
pr := 2*i+1;
j := i;
repeat
preSieve[j] := 0;
inc(j,pr);
until j> High(preSieve);
umf := umf*pr;
end;
inc(i);
until (pr = maxPreSievePrime)OR(umf>High(preSieve)) ;
preSieveOffset := 0;
end;
function InsertSievePrimes(PrimPos:NativeInt):NativeInt;
var
j :NativeUINt;
i,pr : NativeUInt;
begin
i := 0;
//ignore first primes already sieved with
if SieveNum = 0 then
i := maxPreSievePrimeNum;
pr :=0;
j := Uint64(SieveNum)*cSieveSize*2-LastInsertedSievePrime;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i]*2+1;
svdeltaPrime := pr+j;
j := pr;
end;
inc(PrimPos);
for i := i+1 to FoundPrimesCnt-1 do
Begin
IF PrimPos > High(sievePrimes) then
BREAK;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i]*2+1;
svdeltaPrime := (pr-j);
j := pr;
end;
inc(PrimPos);
end;
LastInsertedSievePrime :=Uint64(SieveNum)*cSieveSize*2+pr;
result := PrimPos;
end;
procedure CalcSievePrimOfs(lmt:NativeUint);
//lmt High(sievePrimes)
var
i,pr : NativeUInt;
sq : Uint64;
begin
pr := 0;
i := 0;
repeat
with sievePrimes[i] do
Begin
pr := pr+svdeltaPrime;
IF sqr(pr) < (cSieveSize*2) then
Begin
svSivNum := 0;
svSivOfs := (pr*pr-1) DIV 2;
end
else
Begin
SieveMaxIdx := i;
pr := pr-svdeltaPrime;
BREAK;
end;
end;
inc(i);
until i > lmt;
for i := i to lmt do
begin
with sievePrimes[i] do
Begin
pr := pr+svdeltaPrime;
sq := sqr(pr);
svSivNum := sq DIV (2*cSieveSize);
svSivOfs := ( (sq - Uint64(svSivNum)*(2*cSieveSize))-1)DIV 2;
end;
end;
end;
procedure sievePrimesInit;
var
i,j,pr,PrimPos:NativeInt;
Begin
LastInsertedSievePrime := 0;
preSieveOffset := 0;
SieveNum :=0;
CopyPreSieveInSieve;
//normal sieving of first sieve
i := 1; // start with 3
repeat
while Sieve[i] = 0 do
inc(i);
pr := 2*i+1;
inc(i);
j := ((pr*pr)-1) DIV 2;
if j > High(Sieve) then
BREAK;
repeat
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
until false;
CollectPrimes;
PrimPos := InsertSievePrimes(0);
LastInsertedSievePrime := FoundPrimes[PrimPos]*2+1;
IF PrimPos < High(sievePrimes) then
Begin
Init0Sieve;
sieveOneBlock;
repeat
sieveOneBlock;
dec(SieveNum);
PrimPos := InsertSievePrimes(PrimPos);
inc(SieveNum);
until PrimPos > High(sievePrimes);
end;
Init0Sieve;
end;
procedure Init0Sieve;
begin
FoundPrimesTotal :=0;
preSieveOffset := 0;
SieveNum :=0;
CalcSievePrimOfs(High(sievePrimes));
end;
procedure CopyPreSieveInSieve;
var
lmt : NativeInt;
Begin
lmt := preSieveOffset+cSieveSize;
lmt := lmt-(High(preSieve)+1);
IF lmt<= 0 then
begin
Move(preSieve[preSieveOffset],Sieve[0],cSieveSize);
if lmt <> 0 then
inc(preSieveOffset,cSieveSize)
else
preSieveOffset := 0;
end
else
begin
Move(preSieve[preSieveOffset],Sieve[0],cSieveSize-lmt);
Move(preSieve[0],Sieve[cSieveSize-lmt],lmt);
preSieveOffset := lmt
end;
end;
procedure sieveOneSieve;
var
sp:tpSievePrim;
pSieve :pByte;
i,j,pr,sn,dSievNum :NativeUint;
Begin
pr := 0;
sn := sieveNum;
sp := @sievePrimes[0];
pSieve := @Sieve[0];
For i := SieveMaxIdx downto 0 do
with sp^ do
begin
pr := pr+svdeltaPrime;
IF svSivNum = sn then
Begin
j := svSivOfs;
repeat
pSieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
dSievNum := j DIV cSieveSize;
svSivOfs := j-dSievNum*cSieveSize;
svSivNum := sn+dSievNum;
end;
inc(sp);
end;
i := SieveMaxIdx+1;
repeat
if i > High(SievePrimes) then
BREAK;
with sp^ do
begin
if svSivNum > sn then
Begin
SieveMaxIdx := I-1;
Break;
end;
pr := pr+svdeltaPrime;
j := svSivOfs;
repeat
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
dSievNum := j DIV cSieveSize;
svSivOfs := j-dSievNum*cSieveSize;
svSivNum := sn+dSievNum;
end;
inc(i);
inc(sp);
until false;
end;
procedure CollectPrimes;
//extract primes to FoundPrimes
var
pSieve : pbyte;
pFound : pWord;
i,idx : NativeUint;
Begin
FoundPrimesOffset := SieveNum*2*cSieveSize;
FoundPrimesIdx := 0;
pFound :=@FoundPrimes[0];
i := 0;
idx := 0;
IF SieveNum = 0 then
//include small primes used to pre-sieve
Begin
repeat
pFound[idx]:= (smlPrimes[idx]-1) DIV 2;
inc(idx);
until smlPrimes[idx]>maxPreSievePrime;
i := (smlPrimes[idx] -1) DIV 2;
end;
//grabbing the primes without if then -> reduces time extremly
//primes are born to let branch-prediction fail.
pSieve:= @Sieve[Low(Sieve)];
repeat
//store every value until a prime aka 1 is found
pFound[idx]:= i;
inc(idx,pSieve[i]);
inc(i);
until i>High(Sieve);
FoundPrimesCnt:= idx;
inc(FoundPrimesTotal,Idx);
end;
procedure SieveOneBlock;
begin
CopyPreSieveInSieve;
sieveOneSieve;
CollectPrimes;
inc(SieveNum);
end;
function Nextprime:Uint64;
Begin
result := FoundPrimes[FoundPrimesIdx]*2+1+FoundPrimesOffset;
if (FoundPrimesIdx=0) AND (sievenum = 1) then
inc(result);
inc(FoundPrimesIdx);
If FoundPrimesIdx>= FoundPrimesCnt then
SieveOneBlock;
end;
function PosOfPrime: Uint64;inline;
Begin
result := FoundPrimesTotal-FoundPrimesCnt+FoundPrimesIdx;
end;
function TotalCount :Uint64;inline;
begin
result := FoundPrimesTotal;
end;
function SieveSize :LongInt;inline;
Begin
result := 2*cSieveSize;
end;
function SieveStart:Uint64;inline;
Begin
result := (SieveNum-1)*2*cSieveSize;
end;
procedure InitPrime;
Begin
Init0Sieve;
SieveOneBlock;
end;
//********* segmented sieve of erathostenes *********
const
Limit= 10*1000*1000*1000;
type
tDigits10 = array[0..15] of byte;
td10_UsedDgts2 = array[0..3] of Uint32;
td10_UsedDgts3 = array[0..1] of Uint64;
tpd10_UsedDgts3 = ^td10_UsedDgts3;
procedure OutIn(cnt,p1:NativeInt);
Begin
write(Numb2USA(IntToStr(p1)):13,' ');
if cnt MOD 5 = 0 then
writeln;
end;
function OutByPot10(cnt,prLimit:NativeInt):NativeInt;
Begin
writeln(Numb2USA(IntToStr(cnt)):12,' Ormiston triples before ',Numb2USA(IntToStr(prLimit)):14);
result := 10*prLimit;
end;
function Convert2Digits10(p:NativeUint):Uint64;inline;
const
smlPrimes : array[0..9] of integer = (2,3,5,7,11,13,17,19,23,29);
var
r : NativeUint;
begin
result := 1;
repeat
r := p DIV 10;
result := result*smlPrimes[p-10*r];
p := r;
until r = 0;
end;
{$align 32}
var
IsMod18 : array[0..(1000 DIV 18)*18] of boolean;
procedure OneRun;
var
{$align 32}
p1 : Uint64;
pr1,pr2,pr3,prLimit :nativeUInt;
cnt,prCnt : NativeUint;
begin
pr1 := 0;
pr2 := 0;
prCnt := 0;
prLimit := 10*1000*1000;
repeat
inc(prCnt);
pr3 := pr2;
pr2 := pr1;
pr1 := nextprime;
until pr1 > prLimit;
prLimit *= 10;
cnt := 0;
repeat
inc(prCnt);
if IsMod18[pr2-pr3] then
begin
p1 := Convert2Digits10(pr3);
if Convert2Digits10(pr2) = p1 then
if IsMod18[pr1-pr2] then
begin
if Convert2Digits10(pr1) =p1 then
begin
inc(cnt);
IF cnt <= 25 then
OutIn(cnt,pr1);
end
end
end;
if pr1 >=prLimit then
prlimit:= OutByPot10(cnt,prlimit);
pr3 := pr2;
pr2 := pr1;
pr1 := nextprime;
until pr2 > limit;
writeln(' prime count ',prCnt);
writeln(' last primes ',pr3:12,pr2:12,pr1:12);
end;
var
j :Int32;
Begin
preSieveInit;
sievePrimesInit;
InitPrime;
j := 0;
while j < High(IsMod18)do
Begin
IsMod18[j] := true;
j += 18;
end;
OneRun;
end.
- @TIO.RUN:
11,117,321 12,980,873 14,964,107 32,638,321 32,964,431
33,539,873 35,868,103 44,058,103 46,103,327 48,015,103
50,324,327 52,402,873 58,005,329 60,601,327 61,395,329
74,699,897 76,012,987 78,163,231 80,905,987 81,966,431
82,324,327 82,523,107 83,279,873 86,050,871 92,514,431
25 Ormiston triples before 100,000,000
368 Ormiston triples before 1,000,000,000
4,925 Ormiston triples before 10,000,000,000
prime count 455052513
last primes 9999999967 10000000019 10000000033
Real time: 15.780 s CPU share: 99.20 %
Phix
-- -- demo\rosetta\Ormiston_triplets.exw -- ================================== -- -- Uses a segmented sieve, which is about half the speed of get_primes_le(), but uses far less memory. -- If permited, get_primes_le(1e10) would generate a result of 455,052,511 primes, more than 32 bit -- can cope with, and use over 6GB of ram, and take about 11mins 44s, that is on this box at least, -- whereas this processes them on-the-fly, and only uses about 6MB of memory (ie 0.1% of 6GB). -- with javascript_semantics atom t0 = time() procedure ormiston_triplets(atom limit) // Generate primes using the segmented sieve of Eratosthenes. // credit: https://gist.github.com/kimwalisch/3dc39786fab8d5b34fee integer segment_size = floor(sqrt(limit)), count = 0, i = 3, s = 3, triplen = 1 atom p1 = 2, p2, n = 3, nc = min(1e9,limit), low = 0, t1 = time()+1 sequence isprime = repeat(true,segment_size+1), primes = {}, multiples = {}, orm25 = repeat(0,25) while low<=limit do sequence sieve = repeat(true,segment_size+1) if time()>t1 then progress("Processing %,d/%,d (%3.2f%%)\r",{low,limit,(low/limit)*100}) t1 = time()+1 end if // current segment = [low, high] atom high = min(low+segment_size,limit) // generate sieving primes using simple sieve of Eratosthenes while i*i<=min(high,segment_size) do if isprime[i+1] then for j=i*i to segment_size by i do isprime[j+1] = false end for end if i += 2 end while // initialize sieving primes for segmented sieve while s*s<=high do if isprime[s+1] then primes &= s multiples &= s*s-low end if s += 2 end while // sieve the current segment for mi,j in multiples do integer k = primes[mi]*2 while j<segment_size do sieve[j+1] = false j += k end while multiples[mi] = j - segment_size end for while n<=high do if sieve[n-low+1] then // n is a prime if triplen=1 then if remainder(n-p1,18)=0 and sort(sprint(p1))=sort(sprint(n)) then p2 = n triplen = 2 else p1 = n end if elsif triplen=2 and remainder(n-p2,18)=0 and sort(sprint(p2))=sort(sprint(n)) then -- triplet found! if p1>=nc then string e = elapsed_short(time()-t0) progress("%,d Ormiston triplets before %,d (%s)\n", {count, nc, e}) nc *= 10 end if count += 1 if count<=25 then orm25[count] = sprintf("%d",{p1}) if count=25 then printf(1,"Smallest members of first 25 Ormiston triplets:\n%s\n",join_by(orm25,1,5)) end if end if -- overlapping (and leave triplen set to 2): p1 = p2 p2 = n -- (for disjoint-only just set triplen to 0) else p1 = n triplen = 1 end if end if n += 2 end while low += segment_size end while string e = elapsed_short(time()-t0) progress("%,d Ormiston triplets before %,d (%s)\n", {count, nc, e}) end procedure ormiston_triplets(iff(platform()=JS?1e8:1e9))
- Output:
With limit upped to 1e10
Smallest members of first 25 Ormiston triplets: 11117123 12980783 14964017 32638213 32964341 33539783 35868013 44058013 46103237 48015013 50324237 52402783 58005239 60601237 61395239 74699789 76012879 78163123 80905879 81966341 82324237 82523017 83279783 86050781 92514341 368 Ormiston triplets before 1,000,000,000 (56s) 4,925 Ormiston triplets before 10,000,000,000 (21:32)
Note that running this under pwa/p2js shows 25<1e8 in 8s, a limit of 1e9 would get you a blank screen for 1min 25s
I have not the patience to see whether JavaScript would actually cope with 1e10, but it should (with a blank screen for at least half an hour).
Python
Using Python bindings for the primesieve C++ library, on my machine, this takes about 58 seconds to find Ormiston triples up to one billion, and just over 9 minutes up to 10 billion.
import textwrap
from itertools import pairwise
from typing import Iterator
from typing import List
import primesieve
def primes() -> Iterator[int]:
it = primesieve.Iterator()
while True:
yield it.next_prime()
def triplewise(iterable):
for (a, _), (b, c) in pairwise(pairwise(iterable)):
yield a, b, c
def is_anagram(a: int, b: int, c: int) -> bool:
return sorted(str(a)) == sorted(str(b)) == sorted(str(c))
def up_to_one_billion() -> int:
count = 0
for triple in triplewise(primes()):
if is_anagram(*triple):
count += 1
if triple[2] >= 1_000_000_000:
break
return count
def up_to_ten_billion() -> int:
count = 0
for triple in triplewise(primes()):
if is_anagram(*triple):
count += 1
if triple[2] >= 10_000_000_000:
break
return count
def first_25() -> List[int]:
rv: List[int] = []
for triple in triplewise(primes()):
if is_anagram(*triple):
rv.append(triple[0])
if len(rv) >= 25:
break
return rv
if __name__ == "__main__":
print("Smallest members of first 25 Ormiston triples:")
print(textwrap.fill(" ".join(str(i) for i in first_25())), "\n")
print(up_to_one_billion(), "Ormiston triples before 1,000,000,000")
print(up_to_ten_billion(), "Ormiston triples before 10,000,000,000")
- Output:
Smallest members of first 25 Ormiston triples: 11117123 12980783 14964017 32638213 32964341 33539783 35868013 44058013 46103237 48015013 50324237 52402783 58005239 60601237 61395239 74699789 76012879 78163123 80905879 81966341 82324237 82523017 83279783 86050781 92514341 368 Ormiston triples before 1,000,000,000 4925 Ormiston triples before 10,000,000,000
Raku
use Lingua::EN::Numbers;
use List::Divvy;
my @primes = lazy (^∞).hyper.grep( &is-prime ).map: { $_ => .comb.sort.join };
my @Ormistons = @primes.kv.map: { $^value.key if $^value.value eq @primes[$^key+1].value eq @primes[$^key+2].value};
say "First twenty-five Ormiston triples:";
say @Ormistons[^25].batch(5)».join(', ').join: "\n";
say '';
say +@Ormistons.&before( *[0] > $_ ) ~ " Ormiston triples before " ~ .Int.&cardinal for 1e8, 1e9;
- Output:
First twenty-five Ormiston triples: 11117123, 12980783, 14964017, 32638213, 32964341 33539783, 35868013, 44058013, 46103237, 48015013 50324237, 52402783, 58005239, 60601237, 61395239 74699789, 76012879, 78163123, 80905879, 81966341 82324237, 82523017, 83279783, 86050781, 92514341 25 Ormiston triples before one hundred million 368 Ormiston triples before one billion
Wren
This uses a segmented sieve to get up to 10 billion without running out of memory (bools in Wren require 8 bytes rather than one) and takes just over 13 minutes to achieve the stretch goal.
Limiting the search to a billion, takes about 73 seconds (68 seconds using our 'standard' sieve).
import "./math" for Int
import "./fmt" for Fmt
var limit = 1e10
var primes = Int.segmentedSieve(limit, 8)
var orm25 = []
var j = 1e9
var count = 0
var counts = []
for (i in 0...primes.count-2) {
var p1 = primes[i]
var p2 = primes[i+1]
var p3 = primes[i+2]
if ((p2 - p1) % 18 != 0 || (p3 - p2) % 18 != 0) continue
var key1 = 1
for (dig in Int.digits(p1)) key1 = key1 * primes[dig]
var key2 = 1
for (dig in Int.digits(p2)) key2 = key2 * primes[dig]
if (key1 != key2) continue
var key3 = 1
for (dig in Int.digits(p3)) key3 = key3 * primes[dig]
if (key2 == key3) {
if (count < 25) orm25.add(p1)
if (p1 >= j) {
counts.add(count)
j = j * 10
}
count = count + 1
}
}
counts.add(count)
System.print("Smallest members of first 25 Ormiston triples:")
Fmt.tprint("$,10d ", orm25, 5)
System.print()
j = 1e9
for (i in 0...counts.count) {
Fmt.print("$,d Ormiston triples before $,d", counts[i], j)
j = j * 10
System.print()
}- Output:
Smallest members of first 25 Ormiston triples: 11,117,123 12,980,783 14,964,017 32,638,213 32,964,341 33,539,783 35,868,013 44,058,013 46,103,237 48,015,013 50,324,237 52,402,783 58,005,239 60,601,237 61,395,239 74,699,789 76,012,879 78,163,123 80,905,879 81,966,341 82,324,237 82,523,017 83,279,783 86,050,781 92,514,341 368 Ormiston triples before 1,000,000,000 4,925 Ormiston triples before 10,000,000,000
XPL0
Runs in 87 seconds on Pi4. Beware that MAlloc works differently on other computers.
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;
];
]- Output:
Smallest members of first 25 Ormiston triples: 11117123 12980783 14964017 32638213 32964341 33539783 35868013 44058013 46103237 48015013 50324237 52402783 58005239 60601237 61395239 74699789 76012879 78163123 80905879 81966341 82324237 82523017 83279783 86050781 92514341 368 Ormiston triples before one billion.