Ormiston triples
You are encouraged to solve this task according to the task description, using any language you may know.
- Definition
An Ormiston triple is three consecutive prime numbers which are anagrams, i.e. contain the same decimal digits but in a different order.
- 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
Without external libraries
This example uses a segmented prime iterator, and therefore does not require any external libraries.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
class SegmentedPrimeIterator {
public:
SegmentedPrimeIterator(const uint64_t& aLimit) {
square_root = std::sqrt(aLimit);
index = 0;
low = 0;
high = square_root;
small_sieve(square_root);
}
uint64_t next() {
if ( index == primes.size() ) {
index = 0;
segmented_sieve();
}
return primes[index++];
}
private:
void segmented_sieve() {
low += square_root;
high += square_root;
std::vector<bool> marked_prime(square_root, true);
for ( uint64_t i = 0; i < small_primes.size(); ++i ) {
uint64_t low_limit = ( low / small_primes[i] ) * small_primes[i];
if ( low_limit < low ) {
low_limit += small_primes[i];
}
for ( uint64_t j = low_limit; j < high; j += small_primes[i] ) {
marked_prime[j - low] = false;
}
}
primes.clear();
for ( uint64_t i = low; i < high; ++i ) {
if ( marked_prime[i - low] ) {
primes.emplace_back(i);
}
}
}
void small_sieve(const uint32_t& square_root) {
std::vector<bool> marked_prime(square_root + 1, true);
for ( uint32_t p = 2; p * p <= square_root; ++p ) {
if ( marked_prime[p] ) {
for ( uint32_t i = p * p; i <= square_root; i += p ) {
marked_prime[i] = false;
}
}
}
for ( uint32_t p = 2; p <= square_root; ++p ) {
if ( marked_prime[p] ) {
primes.emplace_back(p);
}
}
small_primes.insert(small_primes.end(), primes.begin(), primes.end());
}
uint32_t index, square_root;
uint64_t low, high;
std::vector<uint64_t> primes;
std::vector<uint64_t> small_primes;
};
std::vector<uint32_t> digits(uint64_t number) {
std::vector<uint32_t> result;
while ( number > 0 ) {
result.emplace_back(number % 10);
number /= 10;
}
std::sort(result.begin(), result.end());
return result;
}
bool haveSameDigits(uint64_t first, uint64_t second) {
return digits(first) == digits(second);
}
int main() {
const uint64_t limit = 10'000'000'000;
SegmentedPrimeIterator primes(limit);
std::vector<uint64_t> ormistons;
uint64_t lower_limit = limit / 10;
uint32_t count = 0;
std::vector<uint32_t> counts;
uint64_t p1 = 0, p2 = 0, p3 = 0;
while ( p3 < limit ) {
p1 = p2; p2 = p3; p3 = primes.next();
if ( ( p2 - p1 ) % 18 != 0 || ( p3 - p2 ) % 18 != 0 ) {
continue;
}
if ( ! haveSameDigits(p1, p2) ) {
continue;
}
if ( haveSameDigits(p2, p3) ) {
if ( count < 25 ) {
ormistons.emplace_back(p1);
}
if ( p1 >= lower_limit ) {
counts.emplace_back(count);
lower_limit *= 10;
}
count += 1;
}
}
counts.emplace_back(count);
std::cout << "The smallest member of the first 25 Ormiston triples:" << std::endl;
for ( uint64_t i = 0; i < ormistons.size(); ++i ) {
std::cout << std::setw(10) << ormistons[i] << ( i % 5 == 4 ? "\n" : "" );
}
std::cout << std::endl;
lower_limit = limit / 10;
for ( uint64_t i = 0; i < counts.size(); ++i ) {
std::cout << counts[i] << " Ormiston triples less than " << lower_limit << std::endl;
lower_limit *= 10;
}
}
- Output:
The smallest member of the 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 less than 1000000000 4925 Ormiston triples less than 10000000000
Delphi
I went a little bit overboard on this problem. Delphi-6 only has a 32-bit memory space, so it would be hard to find triples to 10-billion when the memory only goes up to 2 billion. To solve the problem, I created a bit-boolean array object that stores 8 booleans per byte. I made this into a an object that simulates the functioning of an array, so access to array items is identical accessing a regular boolean array:
PrimeArray[I]:=True;
The whole problem is built as a series of hierchical objects:
1. Bit-Boolean Array Object. This object is responsilbe for accessing bits in each byte. It uses Delphi's array property features, which allows objects to be treated as arrays accessed with the same syntax as normal arrays:
PrimeArray[I]:=PrimeArray[I+1];
2. Sieve Prime Object. This object builds and maintains a list of booliean flags indicating whether a number is prime or not. Since all even numbers are guaranteed not prime, the array only stores odd values. This cuts the space requirements by 50%. However, to the outside world the, array of primes appears to have all the even and odd values. The object uses Delphi's array property feature that allows array syntax to be used, but calls a subroutine to actually actuall read or write the data. The subroutine returns false for every even number and handles special cases such 0, 1 and 2.
3. Ormiston Triple Iterator Object. This object creats a prime table and then allows the user to iterate through all the Ormiston Triples. This makes it easy to perform different task such reading the first 25 triples or rapidly scanning through 1-billion or 10-billion triples to count the one that satisfy the Ormiston criteria.
The program takes 17 seconds to scan 1-billion primes and 10 minutes to test 10 billion.
{Bit boolean - because it stores 8 bools per bytye, it will}
{handle up to 16 gigabyte in a 32 bit programming environment}
type TBitBoolArray = class(TObject)
private
FSize: int64;
ByteArray: array of Byte;
function GetValue(Index: int64): boolean;
procedure WriteValue(Index: int64; const Value: boolean);
function GetSize: int64;
procedure SetSize(const Value: int64);
protected
public
property Value[Index: int64]: boolean read GetValue write WriteValue; default;
constructor Create;
property Count: int64 read GetSize write SetSize;
procedure Clear(Value: boolean);
end;
{ TBitBoolArray }
const BitArray: array [0..7] of byte = ($01, $02, $04, $08, $10, $20, $40, $80);
function TBitBoolArray.GetValue(Index: int64): boolean;
begin
{Note: (Index and 7) is faster than (Index mod 8)}
Result:=(ByteArray[Index shr 3] and BitArray[Index and 7])<>0;
end;
procedure TBitBoolArray.WriteValue(Index: int64; const Value: boolean);
var Inx: int64;
begin
Inx:=Index shr 3;
{Note: (Index and 7) is faster than (Index mod 8)}
if Value then ByteArray[Inx]:=ByteArray[Inx] or BitArray[Index and 7]
else ByteArray[Inx]:=ByteArray[Inx] and not BitArray[Index and 7]
end;
constructor TBitBoolArray.Create;
begin
SetLength(ByteArray,0);
end;
function TBitBoolArray.GetSize: int64;
begin
Result:=FSize;
end;
procedure TBitBoolArray.SetSize(const Value: int64);
var Len: int64;
begin
FSize:=Value;
{Storing 8 items per byte}
Len:=Value div 8;
{We need one more to fill partial bits}
if (Value mod 8)<>0 then Inc(Len);
SetLength(ByteArray,Len);
end;
procedure TBitBoolArray.Clear(Value: boolean);
var Fill: byte;
begin
if Value then Fill:=$FF else Fill:=0;
FillChar(ByteArray[0],Length(ByteArray),Fill);
end;
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
{Sieve object the generates and holds prime values}
{Enable this flag if you need primes past 2 billion.
The flag signals the code to use bit-booleans arrays
which can contain up to 8 x 4 gigabytes = 32 gig booleans.}
{$define BITBOOL}
type TPrimeSieve = class(TObject)
private
{$ifdef BITBOOL}
PrimeArray: TBitBoolArray;
{$else}
PrimeArray: array of boolean;
{$endif}
FArraySize: int64;
function GetPrime(Index: int64): boolean;
procedure Clear;
protected
procedure DoSieve;
property ArraySize: int64 read FArraySize;
public
BitBoolean: boolean;
constructor Create;
destructor Destroy; override;
procedure Intialize(Size: int64);
property Prime[Index: int64]: boolean read GetPrime; default;
function NextPrime(Start: int64): int64;
end;
procedure TPrimeSieve.Clear;
begin
{$ifdef BITBOOL}
PrimeArray.Clear(True);
{$else}
FillChar(PrimeArray[0],Length(PrimeArray),True);
{$endif}
end;
constructor TPrimeSieve.Create;
begin
{$ifdef BITBOOL}
PrimeArray:=TBitBoolArray.Create;
BitBoolean:=True;
{$else}
BitBoolean:=False;
{$endif}
end;
destructor TPrimeSieve.Destroy;
begin
{$ifdef BITBOOL}
PrimeArray.Free;
{$endif}
inherited;
end;
procedure TPrimeSieve.DoSieve;
{Load flags with true/false to flag that number is prime}
{Note: does not store even values, because except for 2, all primes are even}
{Starts storing flags at Index=3, so reading/writing routines compensate}
{Uses for-loops for boolean arrays and while-loops for bitbooleans arrays}
{$ifdef BITBOOL}
var Offset, I, K: int64;
{$else}
var Offset, I, K: cardinal;
{$endif}
begin
Clear;
{$ifdef BITBOOL}
I:=0;
while I<ArraySize do
{$else}
for I:=0 to ArraySize-1 do
{$endif}
begin
if PrimeArray[I] then
begin
Offset:= I + I + 3;
K:= I + Offset;
while K <=(ArraySize-1) do
begin
PrimeArray[K]:= False;
K:= K + Offset;
end;
end;
{$ifdef BITBOOL} Inc(I); {$endif}
end;
end;
function TPrimeSieve.GetPrime(Index: int64): boolean;
{Get a prime flag from array - compensates}
{ for 0,1,2 and even numbers not being stored}
begin
if Index in [0,1,2] then Result:=True
else if (Index and 1)=0 then Result:=false
else Result:=PrimeArray[(Index div 2)-1];
end;
function TPrimeSieve.NextPrime(Start: int64): int64;
{Get next prime after Start}
begin
Result:=Start+1;
while Result<=((ArraySize-1) * 2) do
begin
if Self.Prime[Result] then break;
Inc(Result);
end;
end;
procedure TPrimeSieve.Intialize(Size: int64);
{Set array size and do Sieve to load flag array with}
begin
FArraySize:=Size div 2;
{$ifdef BITBOOL}
PrimeArray.Count:=FArraySize;
{$else}
SetLength(PrimeArray,FArraySize);
{$endif}
DoSieve;
end;
{-------------------------------------------------------------------------------}
type TTripleInfo = record
Prime1,Prime2,Prime3: int64;
Count: int64;
end;
{Iterator for Ormiston Triple}
type TOrm3Iterator = class(TObject)
FInfo: TTripleInfo;
PS: TPrimeSieve;
private
function IsOrmistonTriple(P1, P2, P3: int64): boolean;
function EncodeNumber(N: int64): int64;
protected
public
procedure Reset;
procedure SetSize(Size: int64);
function GetNext(Limit: int64; var Info: TTripleInfo): boolean;
constructor Create;
destructor Destroy; override;
end;
procedure TOrm3Iterator.Reset;
{Restart iterator}
begin
FInfo.Count:=0;
FInfo.Prime1:=1; FInfo.Prime2:=3; FInfo.Prime3:=5;
end;
procedure TOrm3Iterator.SetSize(Size: int64);
begin
PS.Intialize(Size);
Reset;
end;
constructor TOrm3Iterator.Create;
begin
PS:=TPrimeSieve.Create;
{Start with trivial prime set}
SetSize(100);
end;
destructor TOrm3Iterator.Destroy;
begin
PS.Free;
inherited;
end;
function TOrm3Iterator.EncodeNumber(N: int64): int64;
{Encode N by counting digits 0..9 into nibbles}
{Get the product of the integers in a number}
var T: integer;
const NibMap: array [0..9] of int64 = (
{0} $1, {1} $10, {2} $100, {3} $1000, {4} $10000, {5} $100000,
{6} $1000000, {7} $10000000, {8} $100000000, {9} $1000000000);
begin
Result:=0;
repeat
begin
T:=N mod 10;
N:=N div 10;
Result:=Result + NibMap[T];
end
until N<1;
end;
function TOrm3Iterator.IsOrmistonTriple(P1,P2,P3: int64): boolean;
var Pd1,Pd2,Pd3: int64;
begin
Result:=False;
{Optimization - difference in primes should be multiple of 18}
if (((P2 - P1) mod 18)<>0) or (((p3 - p2) mod 18)<>0) then exit;
Pd1:=EncodeNumber(P1);
Pd2:=EncodeNumber(P2);
if Pd1<>Pd2 then exit;
Pd3:=EncodeNumber(P3);
Result:=Pd2=Pd3;
end;
function TOrm3Iterator.GetNext(Limit: int64; var Info: TTripleInfo): boolean;
{Iterate to next Ormiston Pair - automatically stop at prime>Limit}
{Returns false if it hits limit - true if it found Next Ormiston Pair}
begin
Result:=False;
while true do
begin
{Get next set of primes}
FInfo.Prime1:=FInfo.Prime2;
FInfo.Prime2:=FInfo.Prime3;
FInfo.Prime3:=PS.NextPrime(FInfo.Prime3);
{Abort if 3rd prime is ove limit}
if FInfo.Prime3>=Limit then break;
{Test if it is an Ormiston triple}
if IsOrmistonTriple(FInfo.Prime1,FInfo.Prime2,FInfo.Prime3) then
begin
{Return info on triple}
Inc(FInfo.Count);
Info:=FInfo;
Result:=True;
break;
end;
end;
end;
procedure ShowOrmistonTriple(Memo: TMemo);
var I: integer;
var S: string;
var OI: TOrm3Iterator;
var Info: TTripleInfo;
const Limit = 10000000000;
const DisMod = Limit div 10000000;
const Bill = Limit div 1000000000;
var NS: string;
procedure DisplayTitle;
begin
NS:=IntToStr(Bill)+'-Billion';
Memo.Lines.Add('====== Find Ormiston Triples ======');
Memo.Lines.Add('First 25 plus number to '+NS);
S:='Bit-Boolean Array: ';
if OI.PS.BitBoolean then S:=S+'Yes' else S:=S+'No';
Memo.Lines.Add(S);
end;
begin
{Create iterator}
OI:=TOrm3Iterator.Create;
try
DisplayTitle;
Memo.Lines.Add('Sieving Primes');
OI.SetSize(Limit);
Memo.Lines.Add('Finding first 25 Triples');
{Iterate throug 1st 25 tuples}
for I:=1 to 25 do
begin
OI.GetNext(High(Int64),Info);
Memo.Lines.Add(Format('%3d - (%6D %6D %6D) ',[I,Info.Prime1,Info.Prime2,Info.Prime3]));
end;
Memo.Lines.Add('Count='+IntToStr(Info.Count));
Memo.Lines.Add('Counting triples to '+NS);
{Iterate to limit number of triples}
while OI.GetNext(Limit,Info) do
if (Info.Count mod DisMod)=0 then Memo.Lines.Add('Count='+IntToStr(Info.Count));
Memo.Lines.Add(NS+'='+IntToStr(Info.Count));
finally OI.Free; end;
end;
- Output:
====== Find Ormiston Triples ====== First 25 plus number to 1-Billion Bit-Boolean Array: No Sieving Primes Finding first 25 Triples 1 - (11117123 11117213 11117321) 2 - (12980783 12980837 12980873) 3 - (14964017 14964071 14964107) 4 - (32638213 32638231 32638321) 5 - (32964341 32964413 32964431) 6 - (33539783 33539837 33539873) 7 - (35868013 35868031 35868103) 8 - (44058013 44058031 44058103) 9 - (46103237 46103273 46103327) 10 - (48015013 48015031 48015103) 11 - (50324237 50324273 50324327) 12 - (52402783 52402837 52402873) 13 - (58005239 58005293 58005329) 14 - (60601237 60601273 60601327) 15 - (61395239 61395293 61395329) 16 - (74699789 74699879 74699897) 17 - (76012879 76012897 76012987) 18 - (78163123 78163213 78163231) 19 - (80905879 80905897 80905987) 20 - (81966341 81966413 81966431) 21 - (82324237 82324273 82324327) 22 - (82523017 82523071 82523107) 23 - (83279783 83279837 83279873) 24 - (86050781 86050817 86050871) 25 - (92514341 92514413 92514431) Count=25 Counting triples to 1-Billion Count=100 Count=200 Count=300 1-Billion=368 Elapsed Time: 17.303 Sec. ====== Find Ormiston Triples ====== First 25 plus number to 10-Billion Bit-Boolean Array: Yes Sieving Primes Finding first 25 Triples 1 - (11117123 11117213 11117321) 2 - (12980783 12980837 12980873) 3 - (14964017 14964071 14964107) 4 - (32638213 32638231 32638321) 5 - (32964341 32964413 32964431) 6 - (33539783 33539837 33539873) 7 - (35868013 35868031 35868103) 8 - (44058013 44058031 44058103) 9 - (46103237 46103273 46103327) 10 - (48015013 48015031 48015103) 11 - (50324237 50324273 50324327) 12 - (52402783 52402837 52402873) 13 - (58005239 58005293 58005329) 14 - (60601237 60601273 60601327) 15 - (61395239 61395293 61395329) 16 - (74699789 74699879 74699897) 17 - (76012879 76012897 76012987) 18 - (78163123 78163213 78163231) 19 - (80905879 80905897 80905987) 20 - (81966341 81966413 81966431) 21 - (82324237 82324273 82324327) 22 - (82523017 82523071 82523107) 23 - (83279783 83279837 83279873) 24 - (86050781 86050817 86050871) 25 - (92514341 92514413 92514431) Count=25 Counting triples to 10-Billion Count=1000 Count=2000 Count=3000 Count=4000 10-Billion=4925 Elapsed Time: 09:52.882 min
EasyLang
fastfunc isprim num .
if num mod 3 = 0
return 0
.
i = 5
while i <= sqrt num
if num mod i = 0
return 0
.
i += 2
if num mod i = 0
return 0
.
i += 4
.
return 1
.
func nextprim n .
repeat
n += 2
until isprim n = 1
.
return n
.
func digs n .
while n > 0
r += pow 10 (n mod 10)
n = n div 10
.
return r
.
print "Smallest member of the first 25 Ormiston triples:"
b = 2
a = 3
repeat
c = b
dc = db
b = a
db = da
a = nextprim a
until a > 1000000000
da = digs a
if da = db and da = dc
cnt += 1
if cnt <= 25
write c & " "
.
.
.
print "Ormiston triples before 1 billion: " & cntF#
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
FreeBASIC
Function GetSig(Byval N As Uinteger) As Uinteger 'Return signature of N
'A "signature" is the count of each digit in N packed into a 32-bit word
Dim As Uinteger Sig = 0
Do While N > 0
Sig += 1 Shl (N Mod 10)
N \= 10
Loop
Return Sig
End Function
Dim As Uinteger limite = 1e9
Dim Shared As Boolean Sieve(limite)
Sub MakeSieve(Size As Uinteger) 'Make prime number sieve
Dim As Uinteger Prime, i, K
Size \= 2 'ignore even numbers
For i = 0 To Size 'set Sieve flags all true
Sieve(i) = True
Next
For i = 0 To Size
If Sieve(i) Then 'found a prime, which is equal to
Prime = i * 2 + 3 ' twice the index + 3
K = i + Prime 'first multiple to strike off
While K <= Size 'strike off all multiples
Sieve(K) = False
K += Prime
Wend
End If
Next
End Sub
MakeSieve(limite)
Print "Smallest members of first 25 Ormiston triples:"
Dim As Uinteger Cnt = 0, N0 = 0, N1 = 0, Sig0 = 0, Sig1 = 0, N = 3, Sig
Dim As Double t0 = Timer
Do
If Sieve(N \ 2 - 1) Then ' is prime
Sig = GetSig(N)
If Sig = Sig0 Andalso Sig = Sig1 Then
Cnt += 1
If Cnt <= 25 Then
Print Using "#,###,###,###"; N0;
If Cnt Mod 5 = 0 Then Print
End If
End If
Sig0 = Sig1: Sig1 = Sig
N0 = N1: N1 = N
End If
If N >= limite Then Print Cnt; " Ormiston triples before one billion." : Exit Do
N += 2
Loop
Sleep
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
Java
This example uses a segmented prime iterator to complete the stretch task in 60 seconds, without external libraries.
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
public final class OrmistonTriples {
public static void main(String[] aArgs) {
final long limit = 10_000_000_000L;
SegmentedPrimeIterator primes = new SegmentedPrimeIterator(limit);
List<Long> ormistons = new ArrayList<Long>();
long lowerLimit = limit / 10;
int count = 0;
List<Integer> counts = new ArrayList<Integer>();
long p1 = 0, p2 = 0, p3 = 0;
while ( p3 < limit ) {
p1 = p2; p2 = p3; p3 = primes.next();
if ( ( p2 - p1 ) % 18 != 0 || ( p3 - p2 ) % 18 != 0 ) {
continue;
}
if ( ! haveSameDigits(p1, p2) ) {
continue;
}
if ( haveSameDigits(p2, p3) ) {
if ( count < 25 ) {
ormistons.add(p1);
}
if ( p1 >= lowerLimit ) {
counts.add(count);
lowerLimit *= 10;
}
count += 1;
}
}
counts.add(count);
System.out.println("The smallest member of the first 25 Ormiston triples:");
for ( int i = 0; i < ormistons.size(); i++ ) {
System.out.print(String.format("%10d%s", ormistons.get(i), ( i % 5 == 4 ? "\n" : "" )));
}
System.out.println();
lowerLimit = limit / 10;
for ( int i = 0; i < counts.size(); i++ ) {
System.out.println(counts.get(i) + " Ormiston triples less than " + lowerLimit);
lowerLimit *= 10;
}
}
private static boolean haveSameDigits(long aOne, long aTwo) {
return digits(aOne).equals(digits(aTwo));
}
private static List<Integer> digits(long aNumber) {
List<Integer> result = new ArrayList<Integer>();
while ( aNumber > 0 ) {
result.add((int) ( aNumber % 10 ));
aNumber /= 10;
}
Collections.sort(result);
return result;
}
}
final class SegmentedPrimeIterator {
public SegmentedPrimeIterator(long aLimit) {
squareRoot = (int) Math.sqrt(aLimit);
high = squareRoot;
smallSieve(squareRoot);
}
public long next() {
if ( index == primes.size() ) {
index = 0;
segmentedSieve();
}
return primes.get(index++);
}
private void segmentedSieve() {
low += squareRoot;
high += squareRoot;
boolean[] markedPrime = new boolean[squareRoot];
Arrays.fill(markedPrime, true);
for ( int i = 0; i < smallPrimes.size(); i++ ) {
long lowLimit = ( low / smallPrimes.get(i) ) * smallPrimes.get(i);
if ( lowLimit < low ) {
lowLimit += smallPrimes.get(i);
}
for ( long j = lowLimit; j < high; j += smallPrimes.get(i) ) {
markedPrime[(int) ( j - low )] = false;
}
}
primes.clear();
for ( long i = low; i < high; i++ ) {
if ( markedPrime[(int) ( i - low )] ) {
primes.add(i);
}
}
}
private void smallSieve(int aSquareRoot) {
boolean[] markedPrime = new boolean[aSquareRoot + 1];
Arrays.fill(markedPrime, true);
for ( int p = 2; p * p <= aSquareRoot; p++ ) {
if ( markedPrime[p] ) {
for ( int i = p * p; i <= aSquareRoot; i += p ) {
markedPrime[i] = false;
}
}
}
for ( int p = 2; p <= aSquareRoot; p++ ) {
if ( markedPrime[p] ) {
primes.add((long) p);
}
}
smallPrimes.addAll(primes);
}
private static int index, squareRoot;
private static long low, high;
private static List<Long> primes = new ArrayList<Long>();
private static List<Long> smallPrimes = new ArrayList<Long>();
}
- Output:
The smallest member of the 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 less than 1000000000 4925 Ormiston triples less than 10000000000
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
Nim
import std/[algorithm, bitops, math, strformat, strutils]
type Sieve = object
data: seq[byte]
func `[]`(sieve: Sieve; idx: Positive): bool =
## Return value of element at index "idx".
let idx = idx shr 1
let iByte = idx shr 3
let iBit = idx and 7
result = sieve.data[iByte].testBit(iBit)
func `[]=`(sieve: var Sieve; idx: Positive; val: bool) =
## Set value of element at index "idx".
let idx = idx shr 1
let iByte = idx shr 3
let iBit = idx and 7
if val: sieve.data[iByte].setBit(iBit)
else: sieve.data[iByte].clearBit(iBit)
func initSieve(lim: Positive): Sieve =
## Initialize a sieve from 2 to "lim".
result.data = newSeq[byte]((lim + 16) shr 4)
result[1] = true
for n in countup(3, sqrt(lim.toFloat).int, 2):
if not result[n]:
for k in countup(n * n, lim, 2 * n):
result[k] = true
func isPrime(sieve: Sieve; n: int): bool =
## Return true if "n" is prime.
result = if (n and 1) == 0: n == 2 else: not sieve[n]
func nextPrime(sieve: Sieve; n: int): int =
## Return next prime greater than "n".
result = n
while true:
inc result
if sieve.isPrime(result):
return
func digits(n: Positive): seq[byte] =
## Return the sorted list of digits of "n".
var n = n.Natural
while n != 0:
result.add byte(n mod 10)
n = n div 10
result.sort()
proc main() =
const N = 10_000_000_000
let sieve = initSieve(N)
echo "Smallest member of the first 25 Ormiston triples:"
var count = 0
var limit = N div 10
var p1 = 2
while true:
if p1 >= limit:
echo &"Number of Ormiston pairs below {insertSep($limit)}: {count}"
limit *= 10
if limit > N: break
# Check p1 and p2.
let p2 = sieve.nextPrime(p1)
if (p2 - p1) mod 18 != 0:
p1 = p2
continue
# Check p2 and p3.
let p3 = sieve.nextPrime(p2)
if (p3 - p2) mod 18 != 0:
p1 = p3 # Skip p2.
continue
# Check p1.digits and p2.digits.
let d1 = p1.digits
if p2.digits != d1:
p1 = p2
continue
# Check p1.digits and p3.digits.
if p3.digits != d1:
p1 = p3 # Skip p2.
continue
# Ormiston triple found.
inc count
if count <= 25:
stdout.write &"{p1:8}"
stdout.write if count mod 5 == 0: '\n' else: ' '
if count == 25: echo()
# Try next.
p1 = p2
main()
- Output:
Smallest member of the 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 pairs below 1_000_000_000: 368 Number of Ormiston pairs below 10_000_000_000: 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 %
Perl
use strict;
use warnings;
use feature 'say';
use ntheory <primes vecfirstidx>;
my(@O,$pairs);
my @primes = @{ primes(1,1e8) };
my @A = map { join '', sort split '', $_ } @primes;
for (1..$#primes-2) { push @O, $_ if $A[$_] eq $A[$_+1] and $A[$_] eq $A[$_+2] }
say "First 25 Ormiston triples:";
$pairs .= sprintf "%8d, ", $primes[$_] for @O[0..24];
$pairs =~ s/, $//;
say $pairs =~ s/.{50}\K/\n/gr;
for (
[1e8, 'one hundred million'],
[1e9, 'one billion'],
) {
my($limit,$text) = @$_;
my $i = vecfirstidx { $primes[$_] >= $limit } @O;
printf "%3d Ormiston triples before %s\n", $i, $text;
}
- Output:
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 25 Ormiston triples before one hundred million 368 Ormiston triples before one billion
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).
faster version (win64 only)
Uses that primesieve dll I found, almost no memory at all, like less than 10MB, and is over twelve times faster than the above.
requires("1.0.3") requires(WINDOWS) requires(64,true) include builtins/primesieve.e atom t0 = time(), t1 = time()+1, p = primesieve_next_prime(), p1 = p, p2, count = 0, nc = 1e9 sequence orm25 = repeat(0,25) procedure showt() if count=25 then progress("Smallest members of first 25 Ormiston triplets:\n") progress("%s\n",join_by(orm25,1,5)) else string e = elapsed_short(time()-t0) progress("%,6d Ormiston triplets before %,d (%s)\n",{count,nc,e}) nc *= 10 end if end procedure while p<1e10 do p2 = p1 p1 = p p = primesieve_next_prime() if remainder(p-p1,18)=0 and remainder(p1-p2,18)=0 then string s = sort(sprint(p)) if sort(sprint(p1))=s and sort(sprint(p2))=s then if p>=nc then showt() end if count += 1 if count<=25 then orm25[count] = sprintf("%d",{p2}) if count=25 then showt() end if end if end if elsif time()>t1 then progress("%,d\r",{p}) t1 = time()+1 end if end while showt()
- Output:
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 (11s) 4,925 Ormiston triplets before 10,000,000,000 (1:45)
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
Ruby
require 'prime'
def all_anagram?(arr)
sorted = arr.first.digits.sort
arr[1..].all?{|a| a.digits.sort == sorted}
end
res = Prime.each.each_cons(3).lazy.filter_map{|ar| ar.first if all_anagram?(ar)}.first(25)
res.each_slice(5){|slice| puts slice.join(", ")}
n = 1E9
res = Prime.each(n).each_cons(3).count {|ar| all_anagram?(ar)}
puts "\nThere are #{res} Ormiston triples below #{n}"
- 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 There are 368 Ormiston triples below 1000000000.0
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, 128 * 1024)
var orm25 = []
var j = limit/10
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 = limit/10
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
Faster version
This uses our bindings to the C++ library primesieve and, as we're now iterating through the primes rather than sieving for them, uses very little memory compared to the first version.
It's also far quicker - 4.4 seconds to search up to 1 billion and 43.4 seconds to search up to 10 billion.
import "./psieve" for Primes
import "./math" for Int
import "./fmt" for Fmt
var limit = 1e10
var digitPrimes = Int.primeSieve(30)
var it = Primes.iter()
var orm25 = []
var j = limit/10
var count = 0
var counts = []
var p1 = it.next
var p2 = it.next
var p3 = it.next
while (true) {
p1 = p2
p2 = p3
p3 = it.next
if ((p2 - p1) % 18 != 0 || (p3 - p2) % 18 != 0) continue
var key1 = 1
for (dig in Int.digits(p1)) key1 = key1 * digitPrimes[dig]
var key2 = 1
for (dig in Int.digits(p2)) key2 = key2 * digitPrimes[dig]
if (key1 != key2) continue
var key3 = 1
for (dig in Int.digits(p3)) key3 = key3 * digitPrimes[dig]
if (key2 == key3) {
if (count < 25) orm25.add(p1)
if (p1 >= j) {
counts.add(count)
if (j == limit) break
j = j * 10
}
count = count + 1
}
}
System.print("Smallest members of first 25 Ormiston triples:")
Fmt.tprint("$,10d ", orm25, 5)
System.print()
j = limit/10
for (i in 0...counts.count) {
Fmt.print("$,d Ormiston triples before $,d", counts[i], j)
j = j * 10
System.print()
}
- Output:
Same as first version.
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.