Ormiston pairs
Ormiston pairs 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.
An Ormiston pair is two consecutive prime numbers which are anagrams, i.e. contain the same decimal digits but in a different order.
(1913, 1931) is the first such pair.
- Task
- Find and show the first 30 Ormiston pairs.
- Find and show the count of Ormiston pairs up to one million.
- Stretch
- Find and show the count of Ormiston pairs up to ten million.
- See also
ALGOL 68
When running this with ALGOL 68G, you will need to specify a large heap size with e.g., -heap 256M on the ALGOL 68G command.
Uses the "signature" idea from the XPL0 sample and the "difference is 0 MOD 18" filter from the Wren sample.
Also shows the number of prime pairs whose difference is 0 MOD 18 but are not Ormiston pairs.
BEGIN # find some Orimiston pairs - pairs of primes where the first and next #
# prime are anagrams #
PR read "primes.incl.A68" PR # include prime utilities #
INT max prime = 10 000 000; # maximum number we will consider #
INT max digits = BEGIN # count the digits of max prime #
INT v := 1;
INT d := 1;
WHILE ( v *:= 10 ) < max prime DO d +:= 1 OD;
d
END;
[ 0 : 9 ]LONG INT dp; # table of max digit powers for signature creation #
dp[ 0 ] := 1; FOR i TO UPB dp DO dp[ i ] := max digits * dp[ i - 1 ] OD;
[]BOOL prime = PRIMESIEVE max prime;
# construct a list of the primes up to the maximum prime to consider #
[]INT prime list = EXTRACTPRIMESUPTO max prime FROMPRIMESIEVE prime;
# splits n into its digits, returning the sum of their counts, each #
# scaled by the digit power of max digits #
PROC get digits = ( INT n )LONG INT:
BEGIN
INT v := n;
LONG INT result := dp[ v MOD 10 ];
WHILE ( v OVERAB 10 ) > 0 DO
result +:= dp[ v MOD 10 ]
OD;
result
END # get digits # ;
# count the Ormiston pairs #
INT o count := 0;
INT n count := 0;
INT p10 := 100 000;
FOR i TO UPB prime list - 1 DO
INT p1 = prime list[ i ];
INT p2 = prime list[ i + 1 ];
IF ( p2 - p1 ) MOD 18 = 0 THEN
# p2 and p1 might be anagrams #
IF get digits( p1 ) /= get digits( p2 ) THEN
# not an Ormiston pair afterall #
n count +:= 1
ELSE
# p1 and p2 are an Ormiston pair #
o count +:= 1;
IF o count <= 30 THEN
print( ( " (", whole( p1, -5 ), ", ", whole( p2, -5 ), ")"
, IF o count MOD 3 = 0 THEN newline ELSE " " FI
)
)
ELIF p1 >= p10 THEN
print( ( whole( o count - 1, -9 )
, " Ormiston pairs below "
, whole( p10, 0 )
, newline
)
);
p10 *:= 10
FI
FI
FI
OD;
print( ( whole( o count, -9 ), " Ormiston pairs below ", whole( max prime, 0 ), newline ) );
print( ( whole( n count, -9 ), " non-Ormiston ""0 MOD 18"" pairs bwlow ", whole( max prime, 0 ) ) )
END- Output:
( 1913, 1931) (18379, 18397) (19013, 19031)
(25013, 25031) (34613, 34631) (35617, 35671)
(35879, 35897) (36979, 36997) (37379, 37397)
(37813, 37831) (40013, 40031) (40213, 40231)
(40639, 40693) (45613, 45631) (48091, 48109)
(49279, 49297) (51613, 51631) (55313, 55331)
(56179, 56197) (56713, 56731) (58613, 58631)
(63079, 63097) (63179, 63197) (64091, 64109)
(65479, 65497) (66413, 66431) (74779, 74797)
(75913, 75931) (76213, 76231) (76579, 76597)
40 Ormiston pairs below 100000
382 Ormiston pairs below 1000000
3722 Ormiston pairs below 10000000
53369 non-Ormiston "0 MOD 18" pairs bwlow 10000000
C++
#include <array>
#include <iomanip>
#include <iostream>
#include <utility>
#include <primesieve.hpp>
class ormiston_pair_generator {
public:
ormiston_pair_generator() { prime_ = pi_.next_prime(); }
std::pair<uint64_t, uint64_t> next_pair() {
for (;;) {
uint64_t prime = prime_;
auto digits = digits_;
prime_ = pi_.next_prime();
digits_ = get_digits(prime_);
if (digits_ == digits)
return std::make_pair(prime, prime_);
}
}
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_;
uint64_t prime_;
std::array<int, 10> digits_;
};
int main() {
ormiston_pair_generator generator;
int count = 0;
std::cout << "First 30 Ormiston pairs:\n";
for (; count < 30; ++count) {
auto [p1, p2] = generator.next_pair();
std::cout << '(' << std::setw(5) << p1 << ", " << std::setw(5) << p2
<< ')' << ((count + 1) % 3 == 0 ? '\n' : ' ');
}
std::cout << '\n';
for (uint64_t limit = 1000000; limit <= 1000000000; ++count) {
auto [p1, p2] = generator.next_pair();
if (p1 > limit) {
std::cout << "Number of Ormiston pairs < " << limit << ": " << count
<< '\n';
limit *= 10;
}
}
}
- Output:
First 30 Ormiston pairs: ( 1913, 1931) (18379, 18397) (19013, 19031) (25013, 25031) (34613, 34631) (35617, 35671) (35879, 35897) (36979, 36997) (37379, 37397) (37813, 37831) (40013, 40031) (40213, 40231) (40639, 40693) (45613, 45631) (48091, 48109) (49279, 49297) (51613, 51631) (55313, 55331) (56179, 56197) (56713, 56731) (58613, 58631) (63079, 63097) (63179, 63197) (64091, 64109) (65479, 65497) (66413, 66431) (74779, 74797) (75913, 75931) (76213, 76231) (76579, 76597) Number of Ormiston pairs < 1000000: 382 Number of Ormiston pairs < 10000000: 3722 Number of Ormiston pairs < 100000000: 34901 Number of Ormiston pairs < 1000000000: 326926
F#
This task uses Extensible Prime Generator (F#)
// Ormiston pairs. Nigel Galloway: January 31st., 2023
let fG(n,g)=let i=Array.zeroCreate<int>10
let rec fG n g=if g<10 then i[g]<-n i[g] 1 else i[g%10]<-n i[g%10] 1; fG n (g/10)
fG (+) n; fG (-) g; Array.forall ((=)0) i
let oPairs n=n|>Seq.pairwise|>Seq.filter fG
primes32()|>oPairs|>Seq.take 30|>Seq.iter(printf "%A "); printfn ""
printfn $"<1 million: %d{primes32()|>Seq.takeWhile((>)1000000)|>oPairs|>Seq.length}"
printfn $"<10 million: %d{primes32()|>Seq.takeWhile((>)10000000)|>oPairs|>Seq.length}"
printfn $"<100 million: %d{primes32()|>Seq.takeWhile((>)100000000)|>oPairs|>Seq.length}"
printfn $"<1 billion: %d{primes32()|>Seq.takeWhile((>)1000000000)|>oPairs|>Seq.length}"
- Output:
(1913, 1931) (18379, 18397) (19013, 19031) (25013, 25031) (34613, 34631) (35617, 35671) (35879, 35897) (36979, 36997) (37379, 37397) (37813, 37831) (40013, 40031) (40213, 40231) (40639, 40693) (45613, 45631) (48091, 48109) (49279, 49297) (51613, 51631) (55313, 55331) (56179, 56197) (56713, 56731) (58613, 58631) (63079, 63097) (63179, 63197) (64091, 64109) (65479, 65497) (66413, 66431) (74779, 74797) (75913, 75931) (76213, 76231) (76579, 76597) <1 million: 382 <10 million: 3722 <100 million: 34901 <1 billion: 326926
Factor
USING: grouping io kernel lists lists.lazy math math.parser
math.primes.lists math.statistics prettyprint sequences ;
: ormistons ( -- list )
lprimes dup cdr lzip
[ first2 [ >dec histogram ] same? ] lfilter ;
"First 30 Ormiston pairs:" print
30 ormistons ltake list>array 5 group simple-table. nl
ormistons [ first 1e6 < ] lwhile llength pprint bl
"Ormiston pairs less than a million." print
- Output:
First 30 Ormiston pairs:
{ 1913 1931 } { 18379 18397 } { 19013 19031 } { 25013 25031 } { 34613 34631 }
{ 35617 35671 } { 35879 35897 } { 36979 36997 } { 37379 37397 } { 37813 37831 }
{ 40013 40031 } { 40213 40231 } { 40639 40693 } { 45613 45631 } { 48091 48109 }
{ 49279 49297 } { 51613 51631 } { 55313 55331 } { 56179 56197 } { 56713 56731 }
{ 58613 58631 } { 63079 63097 } { 63179 63197 } { 64091 64109 } { 65479 65497 }
{ 66413 66431 } { 74779 74797 } { 75913 75931 } { 76213 76231 } { 76579 76597 }
382 Ormiston pairs less than a million.
Go
package main
import (
"fmt"
"rcu"
)
func main() {
const limit = 1e9
primes := rcu.Primes(limit)
var orm30 [][2]int
j := int(1e5)
count := 0
var counts []int
for i := 0; i < len(primes)-1; i++ {
p1 := primes[i]
p2 := primes[i+1]
if (p2-p1)%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 {
if count < 30 {
orm30 = append(orm30, [2]int{p1, p2})
}
if p1 >= j {
counts = append(counts, count)
j *= 10
}
count++
}
}
counts = append(counts, count)
fmt.Println("First 30 Ormiston pairs:")
for i := 0; i < 30; i++ {
fmt.Printf("%5v ", orm30[i])
if (i+1)%3 == 0 {
fmt.Println()
}
}
fmt.Println()
j = int(1e5)
for i := 0; i < len(counts); i++ {
fmt.Printf("%s Ormiston pairs before %s\n", rcu.Commatize(counts[i]), rcu.Commatize(j))
j *= 10
}
}
- Output:
First 30 Ormiston pairs: [ 1913 1931] [18379 18397] [19013 19031] [25013 25031] [34613 34631] [35617 35671] [35879 35897] [36979 36997] [37379 37397] [37813 37831] [40013 40031] [40213 40231] [40639 40693] [45613 45631] [48091 48109] [49279 49297] [51613 51631] [55313 55331] [56179 56197] [56713 56731] [58613 58631] [63079 63097] [63179 63197] [64091 64109] [65479 65497] [66413 66431] [74779 74797] [75913 75931] [76213 76231] [76579 76597] 40 Ormiston pairs before 100,000 382 Ormiston pairs before 1,000,000 3,722 Ormiston pairs before 10,000,000 34,901 Ormiston pairs before 100,000,000 326,926 Ormiston pairs before 1,000,000,000
J
For this, we would like to be able to test if a prime number is the first value in an Ormiston pair:
isorm=: -:&(/:~)&":&> 4 p: ]
We could also use a routine to organize pairs of numbers as moderate width lines of text:
fmtpairs=: {{ names <@([,',',])&":/"1 y}}
Then the task becomes:
fmtpairs (,. 4 p:]) p:30{.I. isorm i.&.(p:inv) 1e6
1913,1931 18379,18397 19013,19031 25013,25031
34613,34631 35617,35671 35879,35897 36979,36997
37379,37397 37813,37831 40013,40031 40213,40231
40639,40693 45613,45631 48091,48109 49279,49297
51613,51631 55313,55331 56179,56197 56713,56731
58613,58631 63079,63097 63179,63197 64091,64109
65479,65497 66413,66431 74779,74797 75913,75931
76213,76231 76579,76597
+/isorm i.&.(p:inv) 1e6 NB. number of Ormiston pairs less than 1e6
382
+/isorm i.&.(p:inv) 1e7 NB. number of Ormiston pairs less than 1e7
3722
Pascal
Free Pascal
update the digits by adding difference.
program Ormiston;
{$IFDEF FPC}{$MODE DELPHI} {$OPTIMIZATION ON,ALL}{$ENDIF}
{$IFDEF WINDOWS}{$APPLICATION CONSOLE}{$ENDIF}
uses
sysutils,strUtils;
const
Limit= 100*1000*1000;
type
tPrimesSieve = array of boolean;
tpPrimes = pBoolean;
var
PrimeSieve : tPrimesSieve;
//*********** Sieve of erathostenes
procedure ClearAll;
begin
setlength(PrimeSieve,0);
end;
function BuildWheel(pPrimes:tpPrimes;lmt:Uint64): Uint64;
var
wheelSize, wpno, pr, pw, i, k: NativeUint;
wheelprimes: array[0..15] of byte;
begin
pr := 1;//the mother of all numbers 1 ;-)
pPrimes[1] := True;
WheelSize := 1;
wpno := 0;
repeat
Inc(pr);
//pw = pr projected in wheel of wheelsize
pw := pr;
if pw > wheelsize then
Dec(pw, wheelsize);
if pPrimes[pw] then
begin
k := WheelSize + 1;
//turn the wheel (pr-1)-times
for i := 1 to pr - 1 do
begin
Inc(k, WheelSize);
if k < lmt then
move(pPrimes[1], pPrimes[k - WheelSize], WheelSize)
else
begin
move(pPrimes[1], pPrimes[k - WheelSize], Lmt - WheelSize * i);
break;
end;
end;
Dec(k);
if k > lmt then
k := lmt;
wheelPrimes[wpno] := pr;
pPrimes[pr] := False;
Inc(wpno);
WheelSize := k;//the new wheelsize
//sieve multiples of the new found prime
i := pr;
i := i * i;
while i <= k do
begin
pPrimes[i] := False;
Inc(i, pr);
end;
end;
until WheelSize >= lmt;
//re-insert wheel-primes 1 still stays prime
while wpno > 0 do
begin
Dec(wpno);
pPrimes[wheelPrimes[wpno]] := True;
end;
result := pr;
end;
procedure Sieve(pPrimes:tpPrimes;lmt:Uint64);
var
sieveprime, fakt, i: UInt64;
begin
sieveprime := BuildWheel(pPrimes,lmt);
repeat
repeat
Inc(sieveprime);
until pPrimes[sieveprime];
fakt := Lmt div sieveprime;
while Not(pPrimes[fakt]) do
dec(fakt);
if fakt < sieveprime then
BREAK;
i := (fakt + 1) mod 6;
if i = 0 then
i := 4;
repeat
pPrimes[sieveprime * fakt] := False;
repeat
Dec(fakt, i);
i := 6 - i;
until (fakt < sieveprime) OR pPrimes[fakt];
if fakt < sieveprime then
BREAK;
until False;
until False;
pPrimes[1] := False;//remove 1
end;
procedure InitAndGetPrimes;
begin
setlength(PrimeSieve,Limit+2);
Sieve(@PrimeSieve[0],Limit);
PrimeSieve[Limit+1] := true;
end;
//*********** End Sieve of erathostenes
type
tDigits10 = record
d10_UsedDgts : array[0..11] of byte;
d10_dgts :array[0..15] of byte;
d10_maxIdx : Uint32;
end;
procedure Out_Digits10(const p0:tDigits10);
var
i : Int32;
begin
with p0 do
begin
i := d10_maxIdx;
while i >= 0 do
begin
write(d10_dgts[i]);
dec(i);
end;
write(' ');
For i := 0 to 9 do
write(d10_UsedDgts[i]);
writeln;
end;
end;
procedure OutIn(cnt,p1,p2:NativeInt);
Begin
write('[',Numb2USA(IntToStr(p1)):6,'|',Numb2USA(IntToStr(p2)):6,']');
if cnt MOD 5 = 0 then
writeln;
end;
function OutByPot10(cnt,prLimit:NativeInt):NativeInt;
Begin
writeln(Numb2USA(IntToStr(cnt)):12,' Ormiston pairs before ',Numb2USA(IntToStr(prLimit)):14);
result := 10*prLimit;
end;
function Convert2Digits10(p:Uint32):tDigits10;
var
i,r,m10: Uint32;
begin
fillchar(result,SizeOf(result),#0);
with result do
begin
i := 0;
repeat
r := p DIV 10;
m10 := p-10*r;
d10_dgts[i] := m10;
inc(d10_UsedDgts[m10]);
inc(i);
p := r;
until r = 0;
end;
end;
function Add2Digits10(delta:Uint32;const p0:tDigits10):tDigits10;
var
i,r,m10,dgt: Uint32;
begin
result := p0;
with result do
begin
i := 0;
repeat
m10 := d10_dgts[i];
//remove old digit
dec(d10_UsedDgts[m10]);
r := delta DIV 10;
dgt := m10 + delta -10*r;
if dgt >= 10 then
begin
dgt -= 10;
r += 1;
end;
d10_dgts[i] := dgt;
//insert new digit
inc(d10_UsedDgts[dgt]);
delta := r;
inc(i);
until (i> d10_maxIdx) OR (delta = 0);
if (i> d10_maxIdx) AND (delta>0) then
begin
d10_maxIdx += 1;
r := delta DIV 10;
dgt := delta -10*r;
d10_dgts[i] := dgt;
inc(d10_UsedDgts[dgt]);
end;
end;
end;
var
pSieve : tpPrimes;
T1,T0: TDateTime;
p1,p2 :tDigits10;
pr,pr_before,
i,cnt,prLimit :nativeInt;
Begin
T0 := now;
InitAndGetPrimes;
T1 := now;
Writeln('time for sieving ',FormatDateTime('NN:SS.ZZZ',T1-T0));
pSieve := @PrimeSieve[0];
pr_before := 2;
p1 := Convert2Digits10(pr_before);
pr := pr_before;
prLimit := 100*1000;
cnt := 0;
repeat
repeat
inc(pr);
until pSieve[pr];// pr[limit+1] set to true
if pr > limit then
BREAK;
p2 := Add2Digits10(pr-pr_before,p1);
i := 9;
repeat
if p2.d10_UsedDgts[i] <> p1.d10_UsedDgts[i] then
BREAK;
dec(i);
until i < 0;
if i < 0 then
begin
inc(cnt);
IF cnt <= 30 then
OutIn(cnt,pr_before,pr);
// Out_Digits10(p1); write(pr:10,' ');Out_Digits10(p2);
end;
p1 := p2;
pr_before := pr;
if pr >=prLimit then
prlimit:= OutByPot10(cnt,prlimit);
until false;
OutByPot10(cnt,prlimit);
end.- @TIO.RUN:
time for sieving 00:00.333
[ 1,913| 1,931][18,379|18,397][19,013|19,031][25,013|25,031][34,613|34,631]
[35,617|35,671][35,879|35,897][36,979|36,997][37,379|37,397][37,813|37,831]
[40,013|40,031][40,213|40,231][40,639|40,693][45,613|45,631][48,091|48,109]
[49,279|49,297][51,613|51,631][55,313|55,331][56,179|56,197][56,713|56,731]
[58,613|58,631][63,079|63,097][63,179|63,197][64,091|64,109][65,479|65,497]
[66,413|66,431][74,779|74,797][75,913|75,931][76,213|76,231][76,579|76,597]
40 Ormiston pairs before 100,000
382 Ormiston pairs before 1,000,000
3,722 Ormiston pairs before 10,000,000
34,901 Ormiston pairs before 100,000,000
Real time: 0.986 s User time: 0.866 s Sys. time: 0.106 s CPU share: 98.64 %
@home 1E10
time for sieving 00:17.926
...
34,901 Ormiston pairs before 100,000,000
326,926 Ormiston pairs before 1,000,000,000
3,037,903 Ormiston pairs before 10,000,000,000
real 0m30,479s user 0m28,191s sys 0m2,279s
jq
Preliminaries
# Assuming . > 2, return an array, $a, of length .+1 such that
# $a[$i] is $i if $i is prime, and null otherwise.
def primeSieve:
# erase(i) sets .[i*j] to false for integral j > 1
def erase(i):
if .[i] then
reduce range(2; (1 + length) / i) as $j (.; .[i * $j] = null)
else .
end;
(. + 1) as $n
| (($n|sqrt) / 2) as $s
| [null, null, range(2; $n)]
| reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i)) ;The Task
def digits: tostring | explode;
def ormiston_pairs($limit):
($limit | primeSieve | map(select(.))) as $primes
| range(0; $primes|length-1) as $i
| $primes[$i] as $p1
| $primes[$i+1] as $p2
| select( ($p2|digits|sort) == ($p1|digits|sort) )
| [$p1, $p2] ;
def task($limit):
reduce ormiston_pairs($limit) as $pair (
{count:0, orm30: [], counts: [], j: 1e5};
if .count < 30 then .orm30 += [$pair] else . end
| if $pair[0] >= .j
then .counts += [.count]
| .j *= 10
else .
end
| .count += 1 )
| .counts += [.count]
| ("First 30 Ormiston pairs:", (.orm30 | map(tostring) | _nwise(3) | join(" "))),
"",
foreach range(0; .counts|length) as $i (.j = 1e5;
.emit = "\(.counts[$i]) Ormiston pairs before \(.j)"
| .j *= 10;
select(.emit).emit) );
task(1e7) # ten million- Output:
First 30 Ormiston pairs: [1913,1931] [18379,18397] [19013,19031] [25013,25031] [34613,34631] [35617,35671] [35879,35897] [36979,36997] [37379,37397] [37813,37831] [40013,40031] [40213,40231] [40639,40693] [45613,45631] [48091,48109] [49279,49297] [51613,51631] [55313,55331] [56179,56197] [56713,56731] [58613,58631] [63079,63097] [63179,63197] [64091,64109] [65479,65497] [66413,66431] [74779,74797] [75913,75931] [76213,76231] [76579,76597] 40 Ormiston pairs before 100000 382 Ormiston pairs before 1000000 3722 Ormiston pairs before 10000000
Phix
with javascript_semantics atom t0 = time(), t1 = time()+1 constant limit = iff(platform()=JS?1e8:1e9), -- (keep JS<10s) primes = get_primes_le(limit) sequence orm30 = {}, counts = {} integer count = 0, nc = 1e5 for i=1 to length(primes)-1 do integer p1 = primes[i], p2 = primes[i+1] if remainder(p2-p1,18)=0 and sort(sprint(p1))=sort(sprint(p2)) then if count<30 then orm30 &= {sprintf("[%5d %5d]",{p1, p2})} end if if p1>=nc then counts &= count nc *= 10 end if count += 1 end if if time()>t1 then progress("%d/%d\r",{i,length(primes)}) t1 = time()+1 end if end for progress("") counts &= count printf(1,"First 30 Ormiston pairs:\n%s\n",join_by(orm30,1,3)) for i,c in counts do printf(1,"%,d Ormiston pairs before %,d\n", {c, power(10,i+4)}) end for ?elapsed(time()-t0)
- Output:
First 30 Ormiston pairs: [ 1913 1931] [18379 18397] [19013 19031] [25013 25031] [34613 34631] [35617 35671] [35879 35897] [36979 36997] [37379 37397] [37813 37831] [40013 40031] [40213 40231] [40639 40693] [45613 45631] [48091 48109] [49279 49297] [51613 51631] [55313 55331] [56179 56197] [56713 56731] [58613 58631] [63079 63097] [63179 63197] [64091 64109] [65479 65497] [66413 66431] [74779 74797] [75913 75931] [76213 76231] [76579 76597] 40 Ormiston pairs before 100,000 382 Ormiston pairs before 1,000,000 3,722 Ormiston pairs before 10,000,000 34,901 Ormiston pairs before 100,000,000 326,926 Ormiston pairs before 1,000,000,000 "39.8s"
slower but higher limits
-- -- demo\rosetta\Ormiston_pairs.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_pairs(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 atom p1 = 2, n = 3, nc = 1e5, low = 0, t1 = time()+1 sequence isprime = repeat(true,segment_size+1), primes = {}, multiples = {}, orm30 = repeat(0,30) 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 remainder(n-p1,18)=0 and sort(sprint(p1))=sort(sprint(n)) then if p1>=nc then string e = elapsed_short(time()-t0) progress("%,d Ormiston pairs before %,d (%s)\n", {count, nc, e}) nc *= 10 end if count += 1 if count<=30 then orm30[count] = sprintf("[%5d %5d]",{p1, n}) if count=30 then printf(1,"First 30 Ormiston pairs:\n%s\n",join_by(orm30,1,3)) end if end if end if p1 = n end if n += 2 end while low += segment_size end while string e = elapsed_short(time()-t0) progress("%,d Ormiston pairs before %,d (%s)\n", {count, nc, e}) end procedure ormiston_pairs(iff(platform()=JS?1e8:1e9))
- Output:
With limit upped to 1e10
First 30 Ormiston pairs: [ 1913 1931] [18379 18397] [19013 19031] [25013 25031] [34613 34631] [35617 35671] [35879 35897] [36979 36997] [37379 37397] [37813 37831] [40013 40031] [40213 40231] [40639 40693] [45613 45631] [48091 48109] [49279 49297] [51613 51631] [55313 55331] [56179 56197] [56713 56731] [58613 58631] [63079 63097] [63179 63197] [64091 64109] [65479 65497] [66413 66431] [74779 74797] [75913 75931] [76213 76231] [76579 76597] 40 Ormiston pairs before 100,000 (0s) 382 Ormiston pairs before 1,000,000 (0s) 3,722 Ormiston pairs before 10,000,000 (0s) 34,901 Ormiston pairs before 100,000,000 (5s) 326,926 Ormiston pairs before 1,000,000,000 (55s) 3,037,903 Ormiston pairs before 10,000,000,000 (21:57)
Note that running this under pwa/p2js with a limit of 1e9 would get you a blank screen for 1min 25s, hence I've limited it to 1e8 (8s)
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
""" rosettacode.org task Ormiston_pairs """
from sympy import primerange
PRIMES1M = list(primerange(1, 1_000_000))
ASBASE10SORT = [str(sorted(list(str(i)))) for i in PRIMES1M]
ORMISTONS = [(PRIMES1M[i - 1], PRIMES1M[i]) for i in range(1, len(PRIMES1M))
if ASBASE10SORT[i - 1] == ASBASE10SORT[i]]
print('First 30 Ormiston pairs:')
for (i, o) in enumerate(ORMISTONS):
if i < 30:
print(f'({o[0] : 6} {o[1] : 6} )',
end='\n' if (i + 1) % 5 == 0 else ' ')
else:
break
print(len(ORMISTONS), 'is the count of Ormiston pairs up to one million.')
- Output:
First 30 Ormiston pairs: ( 1913 1931 ) ( 18379 18397 ) ( 19013 19031 ) ( 25013 25031 ) ( 34613 34631 ) ( 35617 35671 ) ( 35879 35897 ) ( 36979 36997 ) ( 37379 37397 ) ( 37813 37831 ) ( 40013 40031 ) ( 40213 40231 ) ( 40639 40693 ) ( 45613 45631 ) ( 48091 48109 ) ( 49279 49297 ) ( 51613 51631 ) ( 55313 55331 ) ( 56179 56197 ) ( 56713 56731 ) ( 58613 58631 ) ( 63079 63097 ) ( 63179 63197 ) ( 64091 64109 ) ( 65479 65497 ) ( 66413 66431 ) ( 74779 74797 ) ( 75913 75931 ) ( 76213 76231 ) ( 76579 76597 ) 382 is the count of Ormiston pairs up to one million.
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, @primes[$^key+1].key) if $^value.value eq @primes[$^key+1].value };
say "First thirty Ormiston pairs:";
say @Ormistons[^30].batch(3)».map( { "({.[0].fmt: "%5d"}, {.[1].fmt: "%5d"})" } ).join: "\n";
say '';
say +@Ormistons.&before( *[1] > $_ ) ~ " Ormiston pairs before " ~ .Int.&cardinal for 1e5, 1e6, 1e7;
- Output:
First thirty Ormiston pairs: ( 1913, 1931) (18379, 18397) (19013, 19031) (25013, 25031) (34613, 34631) (35617, 35671) (35879, 35897) (36979, 36997) (37379, 37397) (37813, 37831) (40013, 40031) (40213, 40231) (40639, 40693) (45613, 45631) (48091, 48109) (49279, 49297) (51613, 51631) (55313, 55331) (56179, 56197) (56713, 56731) (58613, 58631) (63079, 63097) (63179, 63197) (64091, 64109) (65479, 65497) (66413, 66431) (74779, 74797) (75913, 75931) (76213, 76231) (76579, 76597) 40 Ormiston pairs before one hundred thousand 382 Ormiston pairs before one million 3722 Ormiston pairs before ten million
Rust
// [dependencies]
// primal = "0.3"
fn get_digits(mut n: usize) -> [usize; 10] {
let mut digits = [0; 10];
while n > 0 {
digits[n % 10] += 1;
n /= 10;
}
digits
}
fn ormiston_pairs() -> impl std::iter::Iterator<Item = (usize, usize)> {
let mut digits = [0; 10];
let mut prime = 0;
let mut primes = primal::Primes::all();
std::iter::from_fn(move || {
for p in primes.by_ref() {
let prime0 = prime;
prime = p;
let digits0 = digits;
digits = get_digits(prime);
if digits == digits0 {
return Some((prime0, prime));
}
}
None
})
}
fn main() {
let mut count = 0;
let mut op = ormiston_pairs();
println!("First 30 Ormiston pairs:");
for (p1, p2) in op.by_ref() {
count += 1;
let c = if count % 3 == 0 { '\n' } else { ' ' };
print!("({:5}, {:5}){}", p1, p2, c);
if count == 30 {
break;
}
}
println!();
let mut limit = 1000000;
for (p1, _) in op.by_ref() {
if p1 > limit {
println!("Number of Ormiston pairs < {}: {}", limit, count);
limit *= 10;
if limit == 10000000000 {
break;
}
}
count += 1;
}
}
- Output:
First 30 Ormiston pairs: ( 1913, 1931) (18379, 18397) (19013, 19031) (25013, 25031) (34613, 34631) (35617, 35671) (35879, 35897) (36979, 36997) (37379, 37397) (37813, 37831) (40013, 40031) (40213, 40231) (40639, 40693) (45613, 45631) (48091, 48109) (49279, 49297) (51613, 51631) (55313, 55331) (56179, 56197) (56713, 56731) (58613, 58631) (63079, 63097) (63179, 63197) (64091, 64109) (65479, 65497) (66413, 66431) (74779, 74797) (75913, 75931) (76213, 76231) (76579, 76597) Number of Ormiston pairs < 1000000: 382 Number of Ormiston pairs < 10000000: 3722 Number of Ormiston pairs < 100000000: 34901 Number of Ormiston pairs < 1000000000: 326926
Wren
import "./math" for Int
import "./fmt" for Fmt
var limit = 1e9
var primes = Int.primeSieve(limit)
var orm30 = []
var j = 1e5
var count = 0
var counts = []
for (i in 0...primes.count-1) {
var p1 = primes[i]
var p2 = primes[i+1]
if ((p2 - p1) % 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) {
if (count < 30) orm30.add([p1, p2])
if (p1 >= j) {
counts.add(count)
j = j * 10
}
count = count + 1
}
}
counts.add(count)
System.print("First 30 Ormiston pairs:")
Fmt.tprint("[$,6d] ", orm30, 3)
System.print()
j = 1e5
for (i in 0...counts.count) {
Fmt.print("$,d Ormiston pairs before $,d", counts[i], j)
j = j * 10
}- Output:
First 30 Ormiston pairs: [ 1,913 1,931] [18,379 18,397] [19,013 19,031] [25,013 25,031] [34,613 34,631] [35,617 35,671] [35,879 35,897] [36,979 36,997] [37,379 37,397] [37,813 37,831] [40,013 40,031] [40,213 40,231] [40,639 40,693] [45,613 45,631] [48,091 48,109] [49,279 49,297] [51,613 51,631] [55,313 55,331] [56,179 56,197] [56,713 56,731] [58,613 58,631] [63,079 63,097] [63,179 63,197] [64,091 64,109] [65,479 65,497] [66,413 66,431] [74,779 74,797] [75,913 75,931] [76,213 76,231] [76,579 76,597] 40 Ormiston pairs before 100,000 382 Ormiston pairs before 1,000,000 3,722 Ormiston pairs before 10,000,000 34,901 Ormiston pairs before 100,000,000 326,926 Ormiston pairs before 1,000,000,000
XPL0
func IsPrime(N); \Return 'true' if N is prime
int N, I;
[if N <= 2 then return N = 2;
if (N&1) = 0 then \even >2\ return false;
for I:= 3 to sqrt(N) do
[if rem(N/I) = 0 then return false;
I:= I+1;
];
return true;
];
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;
];
int Cnt, N, N0, Sig, Sig0;
[Cnt:= 0; N0:= 0; Sig0:= 0; N:= 3;
Format(6, 0);
loop [if IsPrime(N) then
[Sig:= GetSig(N);
if Sig = Sig0 then
[Cnt:= Cnt+1;
if Cnt <= 30 then
[RlOut(0, float(N0)); RlOut(0, float(N));
if rem(Cnt/3) = 0 then CrLf(0) else Text(0, " ");
];
];
Sig0:= Sig;
N0:= N;
];
if N = 1_000_000-1 then
[Text(0, "^m^jOrmiston pairs up to one million: ");
IntOut(0, Cnt);
];
if N = 10_000_000-1 then
[Text(0, "^m^jOrmiston pairs up to ten million: ");
IntOut(0, Cnt);
quit;
];
N:= N+2;
];
]- Output:
1913 1931 18379 18397 19013 19031 25013 25031 34613 34631 35617 35671 35879 35897 36979 36997 37379 37397 37813 37831 40013 40031 40213 40231 40639 40693 45613 45631 48091 48109 49279 49297 51613 51631 55313 55331 56179 56197 56713 56731 58613 58631 63079 63097 63179 63197 64091 64109 65479 65497 66413 66431 74779 74797 75913 75931 76213 76231 76579 76597 Ormiston pairs up to one million: 382 Ormiston pairs up to ten million: 3722