106 lines
3.8 KiB
Text
106 lines
3.8 KiB
Text
--
|
|
-- 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))
|