Ormiston triples: Difference between revisions

Line 1,955:

=={{header|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).~~

-- 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

with javascript_semantics

atom t0 = time()

atom t0 = time()

procedure ormiston_triplets(atom limit)

procedure ormiston_triplets(atom limit)

// Generate primes using the segmented sieve of Eratosthenes.

// Generate primes using the segmented sieve of Eratosthenes.

// credit: https://gist.github.com/kimwalisch/3dc39786fab8d5b34fee

// credit: https://gist.github.com/kimwalisch/3dc39786fab8d5b34fee

integer segment_size = floor(sqrt(limit)),

integer segment_size = floor(sqrt(limit)),

count = 0, i = 3, s = 3, triplen = 1

count = 0, i = 3, s = 3, triplen = 1

atom p1 = 2, p2, n = 3, nc = min(1e9,limit), low = 0, t1 = time()+1

atom p1 = 2, p2, n = 3, nc = min(1e9,limit), low = 0, t1 = time()+1

sequence isprime = repeat(true,segment_size+1),

sequence isprime = repeat(true,segment_size+1),

primes = {},

primes = {},

multiples = {},

multiples = {},

orm25 = repeat(0,25)

orm25 = repeat(0,25)

while low<=limit do

while low<=limit do

sequence sieve = repeat(true,segment_size+1)

sequence sieve = repeat(true,segment_size+1)

if time()>t1 then

if time()>t1 then

progress("Processing %,d/%,d (%3.2f%%)\r",{low,limit,(low/limit)*100})

progress("Processing %,d/%,d (%3.2f%%)\r",{low,limit,(low/limit)*100})

t1 = time()+1

t1 = time()+1

end if

end if

// current segment = [low, high]

// current segment = [low, high]

atom high = min(low+segment_size,limit)

atom high = min(low+segment_size,limit)

// generate sieving primes using simple sieve of Eratosthenes

// generate sieving primes using simple sieve of Eratosthenes

while i*i<=min(high,segment_size) do

while i*i<=min(high,segment_size) do

if isprime[i+1] then

if isprime[i+1] then

for j=i*i to segment_size by i do

for j=i*i to segment_size by i do

isprime[j+1] = false

~~<span~~ ~~style="color:~~ ~~#008080;">end~~ ~~<span~~ ~~style~~=~~"color:~~ ~~#008080;">for~~

isprime[j+1] = false

~~<span~~ ~~style="color:~~ ~~#008080;">~~end~~ if~~

end for

end if

i += 2

i += 2

end while

end while

// initialize sieving primes for segmented sieve

// initialize sieving primes for segmented sieve

while s*s<=high do

while s*s<=high do

if isprime[s+1] then

if isprime[s+1] then

primes &= s

primes &= s

multiples &= s*s-low

multiples &= s*s-low

end if

end if

s += 2

s += 2

end while

end while

// sieve the current segment

// sieve the current segment

for mi,j in multiples do

for mi,j in multiples do

integer k = primes[mi]*2

integer k = primes[mi]*2

while j<segment_size do

while j<segment_size do

sieve[j+1] = false

sieve[j+1] = false

j += k

j += k

end while

end while

multiples[mi] = j - segment_size

multiples[mi] = j - segment_size

end for

end for

while n<=high do

while n<=high do

if sieve[n-low+1] then // n is a prime

if sieve[n-low+1] then // n is a prime

if triplen=1 then

if triplen=1 then

if remainder(n-p1,18)=0

if remainder(n-p1,18)=0

and sort(sprint(p1))=sort(sprint(n)) then

and sort(sprint(p1))=sort(sprint(n)) then

p2 = n

p2 = n

triplen = 2

~~<span~~ ~~style~~=~~"color:~~ ~~#008080;">else~~

triplen = 2

else

p1 = n

~~<span~~ ~~style="color:~~ ~~#008080;">end~~ ~~<span~~ ~~style~~=~~"color:~~ ~~#008080;">if~~

p1 = n

end if

elsif triplen=2

elsif triplen=2

and remainder(n-p2,18)=0

and remainder(n-p2,18)=0

and sort(sprint(p2))=sort(sprint(n)) then

~~-- triplet~~ ~~found!~~

and sort(sprint(p2))=sort(sprint(n)) then

-- triplet found!

if p1>=nc then

if p1>=nc then

string e = elapsed_short(time()-t0)

string e = elapsed_short(time()-t0)

progress("%,d Ormiston triplets before %,d (%s)\n", {count, nc, e})

~~<span~~ ~~style="color:~~ ~~#000000;">nc~~ ~~<span~~ ~~style=~~"~~color:~~ ~~#0000FF;">*=~~ ~~<span~~ ~~style="color: #000000;">10~~

progress("%,d Ormiston triplets before %,d (%s)\n", {count, nc, e})

~~<span~~ ~~style="color:~~ ~~#008080;">end~~ ~~<span~~ ~~style~~=~~"color:~~ ~~#008080;">if~~

nc *= 10

end if

count += 1

count += 1

if count<=25 then

if count<=25 then

orm25[count] = sprintf("%d",{p1})

orm25[count] = sprintf("%d",{p1})

if count=25 then

if count=25 then

printf(1,"Smallest members of first 25 Ormiston triplets:\n%s\n",join_by(orm25,1,5))

~~<span~~ ~~style="color:~~ ~~#008080;~~"~~>end~~ ~~<span~~ ~~style="color:~~ ~~#008080;">if~~

printf(1,"Smallest members of first 25 Ormiston triplets:\n%s\n",join_by(orm25,1,5))

~~<span~~ ~~style="color:~~ ~~#008080;">~~end~~ ~~if~~~~

end if

-- overlapping (and leave triplen set to 2):

-- overlapping (and leave triplen set to 2):

p1 = p2

p1 = p2

p2 = n

p2 = n

-- (for disjoint-only just set triplen to 0)

~~<span~~ ~~style="color:~~ ~~#008080;">else~~

-- (for disjoint-only just set triplen to 0)

else

p1 = n

p1 = n

triplen = 1

~~<span~~ ~~style="color:~~ ~~#008080;">end~~ ~~<span~~ ~~style~~=~~"color:~~ ~~#008080;">if~~

triplen = 1

~~<span~~ ~~style="color:~~ ~~#008080;">~~end~~ ~~if~~~~

end if

n += 2

n += 2

end while

end while

low += segment_size

low += segment_size

end while

end while

string e = elapsed_short(time()-t0)

string e = elapsed_short(time()-t0)

progress("%,d Ormiston triplets before %,d (%s)\n", {count, nc, e})

progress("%,d Ormiston triplets before %,d (%s)\n", {count, nc, e})

end procedure

end procedure

ormiston_triplets(iff(platform()=JS?1e8:1e9))

ormiston_triplets(iff(platform()=JS?1e8:1e9))

</syntaxhighlight>

With limit upped to 1e10

Line 2,081:

Line 2,082:

=== 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("1.0.3")

requires(WINDOWS)

requires(WINDOWS)

requires(64,true)

requires(64,true)

include builtins/primesieve.e

include builtins/primesieve.e

atom t0 = time(), t1 = time()+1,

atom t0 = time(), t1 = time()+1,

p = primesieve_next_prime(),

p = primesieve_next_prime(),

p1 = p, p2, count = 0, nc = 1e9

p1 = p, p2, count = 0, nc = 1e9

sequence orm25 = repeat(0,25)

sequence orm25 = repeat(0,25)

procedure showt()

procedure showt()

if count=25 then

if count=25 then

progress("Smallest members of first 25 Ormiston triplets:\n")

progress("Smallest members of first 25 Ormiston triplets:\n")

progress("%s\n",join_by(orm25,1,5))

progress("%s\n",join_by(orm25,1,5))

else

else

string e = elapsed_short(time()-t0)

string e = elapsed_short(time()-t0)

progress("%,6d Ormiston triplets before %,d (%s)\n",{count,nc,e})

progress("%,6d Ormiston triplets before %,d (%s)\n",{count,nc,e})

nc *= 10

nc *= 10

end if

end if

end procedure

end procedure

while p<1e10 do

while p<1e10 do

p2 = p1

p2 = p1

p1 = p

p1 = p

p = primesieve_next_prime()

p = primesieve_next_prime()

if remainder(p-p1,18)=0

if remainder(p-p1,18)=0

and remainder(p1-p2,18)=0 then

and remainder(p1-p2,18)=0 then

string s = sort(sprint(p))

string s = sort(sprint(p))

if sort(sprint(p1))=s

if sort(sprint(p1))=s

and sort(sprint(p2))=s then

and sort(sprint(p2))=s then

if p>=nc then showt() end if

if p>=nc then showt() end if

count += 1

count += 1

if count<=25 then

if count<=25 then

orm25[count] = sprintf("%d",{p2})

orm25[count] = sprintf("%d",{p2})

if count=25 then showt() end if

if count=25 then showt() end if

end if

end if

end if

end if

elsif time()>t1 then

elsif time()>t1 then

progress("%,d\r",{p})

progress("%,d\r",{p})

t1 = time()+1

t1 = time()+1

end if

end if

end while

end while

showt()

showt()

</syntaxhighlight>

<pre>