Euclid-Mullin sequence: Difference between revisions
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
m (→{{header|Wren}}: Minor tidy) |
||
(36 intermediate revisions by 13 users not shown) | |||
Line 1: | Line 1: | ||
{{ |
{{task}} |
||
;Definition |
;Definition |
||
Line 17: | Line 17: | ||
[https://oeis.org/A000945 OEIS sequence A000945] |
[https://oeis.org/A000945 OEIS sequence A000945] |
||
<br><br> |
<br><br> |
||
=={{header|ALGOL 68}}== |
|||
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}} |
|||
Uses ALGOL 68G's LONG LONG INT which has programmer specifiable precission, the default is sufficient for this task. |
|||
<br> |
|||
Although the first 16 elements will all fit in 64 bits, the product exceeds 64 bits after the ninth element. |
|||
<syntaxhighlight lang="algol68"> |
|||
BEGIN # find elements of the Euclid-Mullin sequence: starting from 2, # |
|||
# the next element is the smallest prime factor of 1 + the product # |
|||
# of the previous elements # |
|||
print( ( " 2" ) ); |
|||
LONG LONG INT product := 2; |
|||
FROM 2 TO 16 DO |
|||
LONG LONG INT next := product + 1; |
|||
# find the first prime factor of next # |
|||
LONG LONG INT p := 3; |
|||
BOOL found := FALSE; |
|||
WHILE p * p <= next AND NOT ( found := next MOD p = 0 ) DO |
|||
p +:= 2 |
|||
OD; |
|||
IF found THEN next := p FI; |
|||
print( ( " ", whole( next, 0 ) ) ); |
|||
product *:= next |
|||
OD |
|||
END |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003 |
|||
</pre> |
|||
=={{header|ALGOL W}}== |
|||
{{Trans|ALGOL 68|but only showing the first 8 elements as Algol W integers are limited to 32-bit.}} |
|||
<syntaxhighlight lang="algolw"> |
|||
begin % find elements of the Euclid-Mullin sequence: starting from 2, % |
|||
% the next element is the smallest prime factor of 1 + the product % |
|||
% of the previous elements % |
|||
integer product; |
|||
write( "2" ); |
|||
product := 2; |
|||
for i := 2 until 8 do begin |
|||
integer nextV, p; |
|||
logical found; |
|||
nextV := product + 1; |
|||
% find the first prime factor of nextV % |
|||
p := 3; |
|||
found := false; |
|||
while p * p <= nextV and not found do begin |
|||
found := nextV rem p = 0; |
|||
if not found then p := p + 2 |
|||
end while_p_squared_le_nextV_and_not_found ; |
|||
if found then nextV := p; |
|||
writeon( i_w := 1, s_w := 0, " ", nextV ); |
|||
product := product * nextV |
|||
end for_i |
|||
end. |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 |
|||
</pre> |
|||
=={{header|AWK}}== |
=={{header|AWK}}== |
||
Line 50: | Line 111: | ||
2 3 7 43 13 53 5 6221671 |
2 3 7 43 13 53 5 6221671 |
||
</pre> |
</pre> |
||
Alternative version. |
|||
{{Trans|ALGOL 68}} |
|||
<syntaxhighlight lang="awk"> |
|||
# find elements of the Euclid-Mullin sequence: starting from 2, |
|||
# the next element is the smallest prime factor of 1 + the product |
|||
# of the previous elements |
|||
BEGIN { |
|||
printf( "2" ); |
|||
product = 2; |
|||
for( i = 2; i <= 8; i ++ ) |
|||
{ |
|||
nextV = product + 1; |
|||
# find the first prime factor of nextV |
|||
p = 3; |
|||
found = 0; |
|||
while( p * p <= nextV && ! ( found = nextV % p == 0 ) ) |
|||
{ |
|||
p += 2; |
|||
} |
|||
if( found ) |
|||
{ |
|||
nextV = p; |
|||
} |
|||
printf( " %d", nextV ); |
|||
product *= nextV |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 |
|||
</pre> |
|||
=={{header|BASIC}}== |
|||
==={{header|Craft Basic}}=== |
|||
<syntaxhighlight lang="basic">define size = 16, em = 0 |
|||
dim list[size] |
|||
let list[0] = 2 |
|||
print 2, " ", |
|||
for i = 1 to 15 |
|||
let k = 3 |
|||
do |
|||
let em = 1 |
|||
for j = 0 to i - 1 |
|||
let em = (em * list[j]) % k |
|||
next j |
|||
let em = (em + 1) % k |
|||
if em = 0 then |
|||
let list[i] = k |
|||
print list[i], " ", |
|||
break |
|||
endif |
|||
let k = k + 2 |
|||
wait |
|||
loop |
|||
next i</syntaxhighlight> |
|||
==={{header|FreeBASIC}}=== |
|||
Naive and takes forever to find the largest term, but does get there in the end. |
|||
<syntaxhighlight lang="freebasic"> |
|||
dim as ulongint E(0 to 15), k |
|||
dim as integer i, em |
|||
E(0) = 2 : print 2 |
|||
for i=1 to 15 |
|||
k=3 |
|||
do |
|||
em = 1 |
|||
for j as uinteger = 0 to i-1 |
|||
em = (em*E(j)) mod k |
|||
next j |
|||
em = (em + 1) mod k |
|||
if em = 0 then |
|||
E(i)=k |
|||
print E(i) |
|||
exit do |
|||
end if |
|||
k = k + 2 |
|||
loop |
|||
next i</syntaxhighlight> |
|||
=={{header|EasyLang}}== |
|||
{{trans|AWK}} |
|||
<syntaxhighlight> |
|||
limit = 8 |
|||
arr[] = [ 2 ] |
|||
write 2 & " " |
|||
for i = 2 to limit |
|||
k = 3 |
|||
repeat |
|||
em = 1 |
|||
for j = 1 to i - 1 |
|||
em = em * arr[j] mod k |
|||
. |
|||
em = (em + 1) mod k |
|||
until em = 0 |
|||
k += 2 |
|||
. |
|||
arr[] &= k |
|||
write k & " " |
|||
. |
|||
</syntaxhighlight> |
|||
=={{header|F_Sharp|F#}}== |
=={{header|F_Sharp|F#}}== |
||
<syntaxhighlight lang="fsharp"> |
<syntaxhighlight lang="fsharp"> |
||
Line 97: | Line 278: | ||
{{out}}<pre> 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre> |
{{out}}<pre> 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre> |
||
=={{header| |
=={{header|Go}}== |
||
{{trans|Wren}} |
|||
Naive and takes forever to find the largest term, but does get there in the end. |
|||
{{libheader|GMP(Go wrapper)}} |
|||
<syntaxhighlight lang="freebasic"> |
|||
This runs in about 54 seconds which, puzzlingly, is a good bit slower than Wren even though both are using GMP and the Pollard's rho algorithm. I have no idea why. |
|||
dim as ulongint E(0 to 15), k |
|||
<syntaxhighlight lang="go">package main |
|||
dim as integer i, em |
|||
E(0) = 2 : print 2 |
|||
for i=1 to 15 |
|||
k=3 |
|||
do |
|||
em = 1 |
|||
for j as uinteger = 0 to i-1 |
|||
em = (em*E(j)) mod k |
|||
next j |
|||
em = (em + 1) mod k |
|||
if em = 0 then |
|||
E(i)=k |
|||
print E(i) |
|||
exit do |
|||
end if |
|||
k = k + 2 |
|||
loop |
|||
next i</syntaxhighlight> |
|||
import ( |
|||
"fmt" |
|||
big "github.com/ncw/gmp" |
|||
"log" |
|||
) |
|||
var ( |
|||
zero = big.NewInt(0) |
|||
one = big.NewInt(1) |
|||
two = big.NewInt(2) |
|||
three = big.NewInt(3) |
|||
four = big.NewInt(4) |
|||
five = big.NewInt(5) |
|||
six = big.NewInt(6) |
|||
ten = big.NewInt(10) |
|||
k100 = big.NewInt(100000) |
|||
) |
|||
func pollardRho(n, c *big.Int) *big.Int { |
|||
g := func(x, y *big.Int) *big.Int { |
|||
x2 := new(big.Int) |
|||
x2.Mul(x, x) |
|||
x2.Add(x2, c) |
|||
return x2.Mod(x2, y) |
|||
} |
|||
x, y, z := big.NewInt(2), big.NewInt(2), big.NewInt(1) |
|||
d := new(big.Int) |
|||
count := 0 |
|||
for { |
|||
x = g(x, n) |
|||
y = g(g(y, n), n) |
|||
d.Sub(x, y) |
|||
d.Abs(d) |
|||
d.Mod(d, n) |
|||
z.Mul(z, d) |
|||
count++ |
|||
if count == 100 { |
|||
d.GCD(nil, nil, z, n) |
|||
if d.Cmp(one) != 0 { |
|||
break |
|||
} |
|||
z.Set(one) |
|||
count = 0 |
|||
} |
|||
} |
|||
if d.Cmp(n) == 0 { |
|||
return zero |
|||
} |
|||
return d |
|||
} |
|||
func smallestPrimeFactorWheel(n, max *big.Int) *big.Int { |
|||
if n.ProbablyPrime(15) { |
|||
return n |
|||
} |
|||
z := new(big.Int) |
|||
if z.Rem(n, two).Cmp(zero) == 0 { |
|||
return two |
|||
} |
|||
if z.Rem(n, three).Cmp(zero) == 0 { |
|||
return three |
|||
} |
|||
if z.Rem(n, five).Cmp(zero) == 0 { |
|||
return five |
|||
} |
|||
k := big.NewInt(7) |
|||
i := 0 |
|||
inc := []*big.Int{four, two, four, two, four, six, two, six} |
|||
for z.Mul(k, k).Cmp(n) <= 0 { |
|||
if z.Rem(n, k).Cmp(zero) == 0 { |
|||
return k |
|||
} |
|||
k.Add(k, inc[i]) |
|||
if k.Cmp(max) > 0 { |
|||
break |
|||
} |
|||
i = (i + 1) % 8 |
|||
} |
|||
return nil |
|||
} |
|||
func smallestPrimeFactor(n *big.Int) *big.Int { |
|||
s := smallestPrimeFactorWheel(n, k100) |
|||
if s != nil { |
|||
return s |
|||
} |
|||
c := big.NewInt(1) |
|||
s = new(big.Int).Set(n) |
|||
for { |
|||
d := pollardRho(n, c) |
|||
if d.Cmp(zero) == 0 { |
|||
if c.Cmp(ten) == 0 { |
|||
log.Fatal("Pollard Rho doesn't appear to be working.") |
|||
} |
|||
c.Add(c, one) |
|||
} else { |
|||
// get the smallest prime factor of 'd' |
|||
factor := smallestPrimeFactorWheel(d, d) |
|||
// check whether n/d has a smaller prime factor |
|||
s = smallestPrimeFactorWheel(n.Quo(n, d), factor) |
|||
if s != nil { |
|||
if s.Cmp(factor) < 0 { |
|||
return s |
|||
} else { |
|||
return factor |
|||
} |
|||
} else { |
|||
return factor |
|||
} |
|||
} |
|||
} |
|||
} |
|||
func main() { |
|||
k := 19 |
|||
fmt.Println("First", k, "terms of the Euclid–Mullin sequence:") |
|||
fmt.Println(2) |
|||
prod := big.NewInt(2) |
|||
z := new(big.Int) |
|||
count := 1 |
|||
for count < k { |
|||
z.Add(prod, one) |
|||
t := smallestPrimeFactor(z) |
|||
fmt.Println(t) |
|||
prod.Mul(prod, t) |
|||
count++ |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
First 19 terms of the Euclid–Mullin sequence: |
|||
2 |
|||
3 |
|||
7 |
|||
43 |
|||
13 |
|||
53 |
|||
5 |
|||
6221671 |
|||
38709183810571 |
|||
139 |
|||
2801 |
|||
11 |
|||
17 |
|||
5471 |
|||
52662739 |
|||
23003 |
|||
30693651606209 |
|||
37 |
|||
1741 |
|||
</pre> |
|||
=={{header|J}}== |
|||
<syntaxhighlight lang="j"> 2x (, <./@(>:&.(*/)))@[&_~ 15 |
|||
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</syntaxhighlight> |
|||
=={{header|Java}}== |
|||
<syntaxhighlight lang="java"> |
|||
import java.math.BigInteger; |
|||
import java.util.ArrayList; |
|||
import java.util.BitSet; |
|||
import java.util.List; |
|||
import java.util.concurrent.ThreadLocalRandom; |
|||
public final class EuclidMullinSequence { |
|||
public static void main(String[] aArgs) { |
|||
primes = listPrimesUpTo(1_000_000); |
|||
System.out.println("The first 27 terms of the Euclid-Mullin sequence:"); |
|||
System.out.print(2 + " "); |
|||
for ( int i = 1; i < 27; i++ ) { |
|||
System.out.print(String.format("%s%s", nextEuclidMullin(), ( i == 14 || i == 27 ) ? "\n" : " ")); |
|||
} |
|||
} |
|||
private static BigInteger nextEuclidMullin() { |
|||
BigInteger smallestPrime = smallestPrimeFactor(product.add(BigInteger.ONE)); |
|||
product = product.multiply(smallestPrime); |
|||
return smallestPrime; |
|||
} |
|||
private static BigInteger smallestPrimeFactor(BigInteger aNumber) { |
|||
if ( aNumber.isProbablePrime(CERTAINTY_LEVEL) ) { |
|||
return aNumber; |
|||
} |
|||
for ( BigInteger prime : primes ) { |
|||
if ( aNumber.mod(prime).signum() == 0 ) { |
|||
return prime; |
|||
} |
|||
} |
|||
BigInteger factor = pollardsRho(aNumber); |
|||
return smallestPrimeFactor(factor); |
|||
} |
|||
private static BigInteger pollardsRho(BigInteger aN) { |
|||
if ( aN.equals(BigInteger.ONE) ) { |
|||
return BigInteger.ONE; |
|||
} |
|||
if ( aN.mod(BigInteger.TWO).signum() == 0 ) { |
|||
return BigInteger.TWO; |
|||
} |
|||
final BigInteger core = new BigInteger(aN.bitLength(), random); |
|||
BigInteger x = new BigInteger(aN.bitLength(), random); |
|||
BigInteger xx = x; |
|||
BigInteger divisor = null; |
|||
do { |
|||
x = x.multiply(x).mod(aN).add(core).mod(aN); |
|||
xx = xx.multiply(xx).mod(aN).add(core).mod(aN); |
|||
xx = xx.multiply(xx).mod(aN).add(core).mod(aN); |
|||
divisor = x.subtract(xx).gcd(aN); |
|||
} while ( divisor.equals(BigInteger.ONE) ); |
|||
return divisor; |
|||
} |
|||
private static List<BigInteger> listPrimesUpTo(int aLimit) { |
|||
BitSet sieve = new BitSet(aLimit + 1); |
|||
sieve.set(2, aLimit + 1); |
|||
for ( int i = 2; i * i <= aLimit; i = sieve.nextSetBit(i + 1) ) { |
|||
for ( int j = i * i; j <= aLimit; j = j + i ) { |
|||
sieve.clear(j); |
|||
} |
|||
} |
|||
List<BigInteger> result = new ArrayList<BigInteger>(sieve.cardinality()); |
|||
for ( int i = 2; i >= 0; i = sieve.nextSetBit(i + 1) ) { |
|||
result.add(BigInteger.valueOf(i)); |
|||
} |
|||
return result; |
|||
} |
|||
private static List<BigInteger> primes; |
|||
private static BigInteger product = BigInteger.TWO; |
|||
private static ThreadLocalRandom random = ThreadLocalRandom.current(); |
|||
private static final int CERTAINTY_LEVEL = 20; |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
The first 27 terms of the Euclid-Mullin sequence: |
|||
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 |
|||
23003 30693651606209 37 1741 1313797957 887 71 7127 109 23 97 159227 |
|||
</pre> |
|||
=={{header|jq}}== |
|||
'''Works with jq, gojq and jaq, that is, the C, Go and Rust implementations of jq.''' (*) |
|||
'''Adapted from [[#Algol_68|Algol 68]]''' |
|||
(*) The precision of jq and jaq is insufficient to compute more than the first nine numbers |
|||
in the sequence (i.e., beyond 38709183810571). gojq's memory consumption |
|||
becomes excessive after the first 16 numbers in the sequence are produced. |
|||
<syntaxhighlight lang=jq> |
|||
# Output: the Euclid-Mullins sequence, beginning with 2 |
|||
def euclid_mullins: |
|||
foreach range(1; infinite|floor) as $i ( { product: 1 }; |
|||
.next = .product + 1 |
|||
# find the first prime factor of .next |
|||
| .p = 3 |
|||
| .found = false |
|||
| until( .p * .p > .next or .found; |
|||
.found = ((.next % .p) == 0) |
|||
| if .found then . else .p += 2 end) |
|||
| if .found then .next = .p else . end |
|||
| .product *= .next) |
|||
| .next ; |
|||
# Produce 16 terms |
|||
limit(16; euclid_mullins) |
|||
</syntaxhighlight> |
|||
{{output}} |
|||
gojq supports unbounded-precision integer arithmetic |
|||
and was accordingly used to produce the following output |
|||
in accordance with the basic task requirements |
|||
Beyond this number, gojq's memory consumption becomes excessive. |
|||
Invocation: $ gojq -n -f euclid-mullin-sequence.algol.jq |
|||
<pre> |
|||
2 |
|||
3 |
|||
7 |
|||
43 |
|||
13 |
|||
53 |
|||
5 |
|||
6221671 |
|||
38709183810571 |
|||
139 |
|||
2801 |
|||
11 |
|||
17 |
|||
5471 |
|||
52662739 |
|||
23003 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
<syntaxhighlight lang="julia">using Primes |
<syntaxhighlight lang="julia">using Primes |
||
Line 150: | Line 622: | ||
<pre>{2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003, 30693651606209, 37, 1741, 1313797957, 887}</pre> |
<pre>{2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003, 30693651606209, 37, 1741, 1313797957, 887}</pre> |
||
Others may be found by adjusting the range of the Do loop but it will take a while. |
Others may be found by adjusting the range of the Do loop but it will take a while. |
||
=={{header|Lua}}== |
|||
{{Trans|ALGOL W}} |
|||
<syntaxhighlight lang="lua"> |
|||
-- find elements of the Euclid-Mullin sequence: starting from 2, |
|||
-- the next element is the smallest prime factor of 1 + the product |
|||
-- of the previous elements |
|||
do |
|||
io.write( "2" ) |
|||
local product = 2 |
|||
for i = 2, 8 do |
|||
local nextV = product + 1 |
|||
-- find the first prime factor of nextV |
|||
local p = 3 |
|||
local found = false |
|||
while p * p <= nextV and not found do |
|||
found = nextV % p == 0 |
|||
if not found then p = p + 2 end |
|||
end |
|||
if found then nextV = p end |
|||
io.write( " ", nextV ) |
|||
product = product * nextV |
|||
end |
|||
end |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 |
|||
</pre> |
|||
Alternative using Pollard's Rho algorithm. Uses the iterative gcd function from the [[Greatest common divisor]] task.<br> |
|||
{{Trans|Ruby|for the Pollard's Rho algorithm}}. |
|||
Note that, as discussed on the Talk page, Pollard's Rho algorithm won't necessarily find the lowest factor, however it does for the first 16 elements.<br> |
|||
As with the other Lua sample, only 8 elements are found due to the size of some of the subsequent ones. |
|||
<syntaxhighlight lang="lua"> |
|||
function gcd(a,b) |
|||
while b~=0 do |
|||
a,b=b,a%b |
|||
end |
|||
return math.abs(a) |
|||
end |
|||
function pollard_rho(n) |
|||
local x, y, d = 2, 2, 1 |
|||
local g = function(x) return (x*x+1) % n end |
|||
while d == 1 do |
|||
x = g(x) |
|||
y = g(g(y)) |
|||
d = gcd(math.abs(x-y),n) |
|||
end |
|||
if d == n then return d end |
|||
return math.min(d, math.floor( n/d ) ) |
|||
end |
|||
local ar, product = {2}, 2 |
|||
repeat |
|||
ar[ #ar + 1 ] = pollard_rho( product + 1 ) |
|||
product = product * ar[ #ar ] |
|||
until #ar >= 8 |
|||
print( table.concat(ar, " ") ) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 |
|||
</pre> |
|||
=={{header|Maxima}}== |
|||
<syntaxhighlight lang="maxima"> |
|||
euclid_mullin(n):=if n=1 then 2 else ifactors(1+product(euclid_mullin(i),i,1,n-1))[1][1]$ |
|||
/* Test case */ |
|||
makelist(euclid_mullin(k),k,16); |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[2,3,7,43,13,53,5,6221671,38709183810571,139,2801,11,17,5471,52662739,23003] |
|||
</pre> |
|||
=={{header|MiniScript}}== |
|||
{{Trans|Lua}} |
|||
<syntaxhighlight lang="miniscript"> |
|||
// find elements of the Euclid-Mullin sequence: starting from 2, |
|||
// the next element is the smallest prime factor of 1 + the product |
|||
// of the previous elements |
|||
seq = [2] |
|||
product = 2 |
|||
for i in range( 2, 8 ) |
|||
nextV = product + 1 |
|||
// find the first prime factor of nextV |
|||
p = 3 |
|||
found = false |
|||
while p * p <= nextV and not found |
|||
found = nextV % p == 0 |
|||
if not found then p = p + 2 |
|||
end while |
|||
if found then nextV = p |
|||
seq.push( nextV ) |
|||
product = product * nextV |
|||
end for |
|||
print seq.join( " ") |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 |
|||
</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|Wren}} |
|||
{{libheader|integers}} |
|||
<syntaxhighlight lang="Nim">import integers |
|||
let |
|||
Zero = newInteger() |
|||
One = newInteger(1) |
|||
Two = newInteger(2) |
|||
Three = newInteger(3) |
|||
Five = newInteger(5) |
|||
Ten = newInteger(10) |
|||
Max = newInteger(100000) |
|||
None = newInteger(-1) |
|||
proc pollardRho(n, c: Integer): Integer = |
|||
template g(x: Integer): Integer = (x * x + c) mod n |
|||
var |
|||
x = newInteger(2) |
|||
y = newInteger(2) |
|||
z = newInteger(1) |
|||
d = Max + 1 |
|||
count = 0 |
|||
while true: |
|||
x = g(x) |
|||
y = g(g(y)) |
|||
d = abs(x - y) mod n |
|||
z *= d |
|||
inc count |
|||
if count == 100: |
|||
d = gcd(z, n) |
|||
if d != One: break |
|||
z = newInteger(1) |
|||
count = 0 |
|||
result = if d == n: Zero else: d |
|||
template isEven(n: Integer): bool = isZero(n and 1) |
|||
proc smallestPrimeFactorWheel(n: Integer): Integer = |
|||
if n.isPrime(5): return n |
|||
if n.isEven: return Two |
|||
if isZero(n mod 3): return Three |
|||
if isZero(n mod 5): return Five |
|||
var k = newInteger(7) |
|||
var i = 0 |
|||
const Inc = [4, 2, 4, 2, 4, 6, 2, 6] |
|||
while k * k <= n: |
|||
if isZero(n mod k): return k |
|||
k += Inc[i] |
|||
if k > Max: return None |
|||
i = (i + 1) mod 8 |
|||
proc smallestPrimeFactor(n: Integer): Integer = |
|||
var n = n |
|||
result = smallestPrimeFactorWheel(n) |
|||
if result != None: return |
|||
var c = One |
|||
result = newInteger(n) |
|||
while n > Max: |
|||
var d = pollardRho(n, c) |
|||
if d.isZero: |
|||
if c == Ten: |
|||
quit "Pollard Rho doesn't appear to be working.", QuitFailure |
|||
inc c |
|||
else: |
|||
# Can't be sure PR will find the smallest prime factor first. |
|||
result = min(result, d) |
|||
n = n div d |
|||
if n.isPrime(2): |
|||
return min(result, n) |
|||
proc main() = |
|||
var k = 19 |
|||
echo "First ", k, " terms of the Euclid–Mullin sequence:" |
|||
echo 2 |
|||
var prod = newInteger(2) |
|||
var count = 1 |
|||
while count < k: |
|||
let t = smallestPrimeFactor(prod + One) |
|||
echo t |
|||
prod *= t |
|||
inc count |
|||
main() |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
The first 16 numbers are displayed instantly, but it took about 9 seconds on my (moderately powerful) laptop to get the next three. I gave up for the 20th number. |
|||
<pre>First 19 terms of the Euclid–Mullin sequence: |
|||
2 |
|||
3 |
|||
7 |
|||
43 |
|||
13 |
|||
53 |
|||
5 |
|||
6221671 |
|||
38709183810571 |
|||
139 |
|||
2801 |
|||
11 |
|||
17 |
|||
5471 |
|||
52662739 |
|||
23003 |
|||
30693651606209 |
|||
37 |
|||
1741 |
|||
</pre> |
|||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
Line 226: | Line 913: | ||
{{out}} |
{{out}} |
||
<pre>First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre> |
<pre>First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre> |
||
=={{header|Ring}}== |
|||
{{Trans|Lua}} |
|||
<syntaxhighlight lang="ring"> |
|||
// find elements of the Euclid-Mullin sequence: starting from 2, |
|||
// the next element is the smallest prime factor of 1 + the product |
|||
// of the previous elements |
|||
see "2" |
|||
product = 2 |
|||
for i = 2 to 8 |
|||
nextV = product + 1 |
|||
// find the first prime factor of nextV |
|||
p = 3 |
|||
found = false |
|||
while p * p <= nextV and not found |
|||
found = ( nextV % p ) = 0 |
|||
if not found p = p + 2 ok |
|||
end |
|||
if found nextV = p ok |
|||
see " " + nextV |
|||
product = product * nextV |
|||
next |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2 3 7 43 13 53 5 6221671 |
|||
</pre> |
|||
=={{header|RPL}}== |
|||
{{works with|Halcyon Calc|4.2.7}} |
|||
{| class="wikitable" |
|||
! RPL code |
|||
! Comment |
|||
|- |
|||
| |
|||
≪ |
|||
'''IF''' # 1d DUP2 AND ≠ '''THEN''' DROP # 2d |
|||
'''ELSE IF''' DUP 3 DUP2 / * == '''THEN''' DROP # 3d |
|||
'''ELSE''' DUP B→R √ → divm |
|||
≪ 4 5 divm '''FOR''' n |
|||
'''IF''' OVER n DUP2 / * == |
|||
'''THEN''' SWAP DROP n R→B SWAP divm 'n' STO '''END''' |
|||
6 SWAP - DUP '''STEP''' DROP |
|||
≫ '''END END''' |
|||
≫ ‘'''bDIV1'''’ STO |
|||
≪ |
|||
DUP SIZE 1 1 ROT '''FOR''' j OVER j GET * '''NEXT''' |
|||
1 + '''bDIV1''' + |
|||
≫ ''''NXTEM'''' STO |
|||
| |
|||
'''bDIV1''' ''( #m -- #first_divisor )'' |
|||
is #2 a divisor ? |
|||
is #3 a divisor ? |
|||
otherwise get sqrt(m) |
|||
d = 4 ; for n = 5 to sqrt(m) |
|||
if n divides m |
|||
replace m by n and prepare loop exit |
|||
d = 6 - d ; n += d |
|||
'''NXTEM''' ''( { #EM(1) .. #EM(n) } -- { #EM(1) .. #EM(n+1) } )'' |
|||
get EM(1)*..*EM(n) |
|||
get least prime factor of 1+EM(1)*..*EM(n) and add to list |
|||
|} |
|||
{{in}} |
|||
<pre> |
|||
≪ { # 2 } WHILE DUP SIZE ≤ 16 REPEAT NXTEM END ≫ EVAL |
|||
</pre> |
|||
The emulator's watchdog timer prevents checking the primality of EM(9) = # 38709183810571d. Even if this device stayed idle, EM(10) could not be calculated, since the product of all earlier elements is more than 64 bits long. |
|||
{{out}} |
|||
<pre> |
|||
1: { # 2d # 3d # 7d # 43d # 13d # 53d # 5d # 6221671d } |
|||
</pre> |
|||
=={{header|Ruby}}== |
|||
<syntaxhighlight lang="ruby">def pollard_rho(n) |
|||
x, y, d = 2, 2, 1 |
|||
g = proc{|x|(x*x+1) % n} |
|||
while d == 1 do |
|||
x = g[x] |
|||
y = g[g[y]] |
|||
d = (x-y).abs.gcd(n) |
|||
end |
|||
return d if d == n |
|||
[d, n/d].compact.min |
|||
end |
|||
ar = [2] |
|||
ar << pollard_rho(ar.inject(&:*)+1) until ar.size >= 16 |
|||
puts ar.join(", ") |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003 |
|||
</pre> |
|||
=={{header|SETL}}== |
|||
<syntaxhighlight lang="setl">program euclid_mullin; |
|||
print(2); |
|||
product := 2; |
|||
loop for i in [2..16] do |
|||
next := smallest_factor(product + 1); |
|||
product *:= next; |
|||
print(next); |
|||
end loop; |
|||
op smallest_factor(n); |
|||
if even n then return 2; end if; |
|||
d := 3; |
|||
loop while d*d <= n do |
|||
if n mod d=0 then return d; end if; |
|||
d +:= 2; |
|||
end loop; |
|||
return n; |
|||
end op; |
|||
end program;</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 |
|||
3 |
|||
7 |
|||
43 |
|||
13 |
|||
53 |
|||
5 |
|||
6221671 |
|||
38709183810571 |
|||
139 |
|||
2801 |
|||
11 |
|||
17 |
|||
5471 |
|||
52662739 |
|||
23003</pre> |
|||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
Line 244: | Line 1,064: | ||
===Wren-cli=== |
===Wren-cli=== |
||
This uses the [https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm Pollard Rho algorithm] to try and speed up the factorization of the 15th element but overall time still slow at around 32 seconds. |
This uses the [https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm Pollard Rho algorithm] to try and speed up the factorization of the 15th element but overall time still slow at around 32 seconds. |
||
<syntaxhighlight lang=" |
<syntaxhighlight lang="wren">import "./big" for BigInt |
||
var zero = BigInt.zero |
var zero = BigInt.zero |
||
Line 250: | Line 1,070: | ||
var two = BigInt.two |
var two = BigInt.two |
||
var ten = BigInt.ten |
var ten = BigInt.ten |
||
var |
var k100 = BigInt.new(100000) |
||
var |
var smallestPrimeFactorWheel = Fn.new { |n, max| |
||
var g = Fn.new { |x, y| (x*x + c) % n } |
|||
var x = two |
|||
var y = two |
|||
var z = one |
|||
var d = max + one |
|||
var count = 0 |
|||
while (true) { |
|||
x = g.call(x, n) |
|||
y = g.call(g.call(y, n), n) |
|||
d = (x - y).abs % n |
|||
z = z * d |
|||
count = count + 1 |
|||
if (count == 100) { |
|||
d = BigInt.gcd(z, n) |
|||
if (d != one) break |
|||
z = one |
|||
count = 0 |
|||
} |
|||
} |
|||
if (d == n) return zero |
|||
return d |
|||
} |
|||
var smallestPrimeFactorWheel = Fn.new { |n| |
|||
if (n.isProbablePrime(5)) return n |
if (n.isProbablePrime(5)) return n |
||
if (n % 2 == zero) return BigInt.two |
if (n % 2 == zero) return BigInt.two |
||
Line 293: | Line 1,089: | ||
var smallestPrimeFactor = Fn.new { |n| |
var smallestPrimeFactor = Fn.new { |n| |
||
var s = smallestPrimeFactorWheel.call(n) |
var s = smallestPrimeFactorWheel.call(n, k100) |
||
if (s) return s |
if (s) return s |
||
var c = one |
var c = one |
||
while (true) { |
|||
var d = BigInt.pollardRho(n, 2, c) |
|||
while (n > max) { |
|||
var d = pollardRho.call(n, c) |
|||
if (d == 0) { |
if (d == 0) { |
||
if (c == ten) Fiber.abort("Pollard Rho doesn't appear to be working.") |
if (c == ten) Fiber.abort("Pollard Rho doesn't appear to be working.") |
||
c = c + one |
c = c + one |
||
} else { |
} else { |
||
// |
// get the smallest prime factor of 'd' |
||
var factor = smallestPrimeFactorWheel.call(d, d) |
|||
// check whether n/d has a smaller prime factor |
|||
s = smallestPrimeFactorWheel.call(n/d, factor) |
|||
return s ? BigInt.min(s, factor) : factor |
|||
} |
} |
||
} |
} |
||
return s |
|||
} |
} |
||
Line 322: | Line 1,117: | ||
prod = prod * t |
prod = prod * t |
||
count = count + 1 |
count = count + 1 |
||
} |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 345: | Line 1,140: | ||
</pre> |
</pre> |
||
<br> |
<br> |
||
===Embedded=== |
===Embedded=== |
||
{{libheader|Wren-gmp}} |
{{libheader|Wren-gmp}} |
||
This finds the first 16 in 0.11 seconds |
This finds the first 16 in 0.11 seconds but the next 11 takes 6 minutes 17 seconds. |
||
<syntaxhighlight lang="ecmascript">/* euclid_mullin_gmp.wren */ |
|||
If we could assume that Pollard's Rho will always find the smallest prime factor (or a multiple thereof) first, then this would bring the runtime down to 44 seconds and still produce the correct answers for this particular task. However, in general it is not safe to assume that - see discussion in Talk Page. |
|||
<syntaxhighlight lang="wren">/* Euclid_mullin_sequence_2.wren */ |
|||
import "./gmp" for Mpz |
import "./gmp" for Mpz |
||
var |
var k100 = Mpz.from(100000) |
||
var |
var smallestPrimeFactorTrial = Fn.new { |n, max| |
||
if (n.probPrime(15) > 0) return n |
if (n.probPrime(15) > 0) return n |
||
var k = Mpz.one |
|||
if (n.isDivisibleUi(3)) return Mpz.three |
|||
if (n.isDivisibleUi(5)) return Mpz.five |
|||
var k = Mpz.from(7) |
|||
var i = 0 |
|||
var inc = [4, 2, 4, 2, 4, 6, 2, 6] |
|||
while (k * k <= n) { |
while (k * k <= n) { |
||
k.nextPrime |
|||
k.add(inc[i]) |
|||
if (k > max) return null |
if (k > max) return null |
||
if (n.isDivisible(k)) return k |
|||
} |
} |
||
} |
} |
||
var smallestPrimeFactor = Fn.new { |n| |
var smallestPrimeFactor = Fn.new { |n| |
||
var s = |
var s = smallestPrimeFactorTrial.call(n, k100) |
||
if (s) return s |
if (s) return s |
||
var c = Mpz.one |
var c = Mpz.one |
||
while (true) { |
|||
while (n > max) { |
|||
var d = Mpz.pollardRho(n, 2, c) |
var d = Mpz.pollardRho(n, 2, c) |
||
if (d.isZero) { |
if (d.isZero) { |
||
Line 381: | Line 1,172: | ||
c.inc |
c.inc |
||
} else { |
} else { |
||
// |
// get the smallest prime factor of 'd' |
||
var factor = smallestPrimeFactorTrial.call(d, d) |
|||
n |
// check whether n/d has a smaller prime factor |
||
s = smallestPrimeFactorTrial.call(n/d, factor) |
|||
return s ? Mpz.min(s, factor) : factor |
|||
} |
} |
||
} |
} |
||
return s |
|||
} |
} |
||
var k = |
var k = 27 |
||
System.print("First %(k) terms of the Euclid–Mullin sequence:") |
System.print("First %(k) terms of the Euclid–Mullin sequence:") |
||
System.print(2) |
System.print(2) |
||
Line 403: | Line 1,194: | ||
{{out}} |
{{out}} |
||
As Wren-cli plus |
As Wren-cli plus eleven more: |
||
<pre> |
<pre> |
||
30693651606209 |
30693651606209 |
||
37 |
37 |
||
1741 |
1741 |
||
1313797957 |
|||
887 |
|||
71 |
|||
7127 |
|||
109 |
|||
23 |
|||
97 |
|||
159227 |
|||
</pre> |
</pre> |
Latest revision as of 11:08, 30 November 2023
You are encouraged to solve this task according to the task description, using any language you may know.
- Definition
The Euclid–Mullin sequence is an infinite sequence of distinct prime numbers, in which each element is the least prime factor of one plus the product of all earlier elements.
The first element is usually assumed to be 2. So the second element is : (2) + 1 = 3 and the third element is : (2 x 3) + 1 = 7 as this is prime.
Although intermingled with smaller elements, the sequence can produce very large elements quite quickly and only the first 51 have been computed at the time of writing.
- Task
Compute and show here the first 16 elements of the sequence or, if your language does not support arbitrary precision arithmetic, as many as you can.
- Stretch goal
Compute the next 11 elements of the sequence.
- Reference
ALGOL 68
Uses ALGOL 68G's LONG LONG INT which has programmer specifiable precission, the default is sufficient for this task.
Although the first 16 elements will all fit in 64 bits, the product exceeds 64 bits after the ninth element.
BEGIN # find elements of the Euclid-Mullin sequence: starting from 2, #
# the next element is the smallest prime factor of 1 + the product #
# of the previous elements #
print( ( " 2" ) );
LONG LONG INT product := 2;
FROM 2 TO 16 DO
LONG LONG INT next := product + 1;
# find the first prime factor of next #
LONG LONG INT p := 3;
BOOL found := FALSE;
WHILE p * p <= next AND NOT ( found := next MOD p = 0 ) DO
p +:= 2
OD;
IF found THEN next := p FI;
print( ( " ", whole( next, 0 ) ) );
product *:= next
OD
END
- Output:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
ALGOL W
begin % find elements of the Euclid-Mullin sequence: starting from 2, %
% the next element is the smallest prime factor of 1 + the product %
% of the previous elements %
integer product;
write( "2" );
product := 2;
for i := 2 until 8 do begin
integer nextV, p;
logical found;
nextV := product + 1;
% find the first prime factor of nextV %
p := 3;
found := false;
while p * p <= nextV and not found do begin
found := nextV rem p = 0;
if not found then p := p + 2
end while_p_squared_le_nextV_and_not_found ;
if found then nextV := p;
writeon( i_w := 1, s_w := 0, " ", nextV );
product := product * nextV
end for_i
end.
- Output:
2 3 7 43 13 53 5 6221671
AWK
# syntax: GAWK -f EUCLID-MULLIN_SEQUENCE.AWK
# converted from FreeBASIC
BEGIN {
limit = 7 # we'll stop here
arr[0] = 2
printf("%s ",arr[0])
for (i=1; i<=limit; i++) {
k = 3
while (1) {
em = 1
for (j=0; j<=i-1; j++) {
em = (em * arr[j]) % k
}
em = (em + 1) % k
if (em == 0) {
arr[i] = k
printf("%s ",arr[i])
break
}
k += 2
}
}
printf("\n")
exit(0)
}
- Output:
2 3 7 43 13 53 5 6221671
Alternative version.
# find elements of the Euclid-Mullin sequence: starting from 2,
# the next element is the smallest prime factor of 1 + the product
# of the previous elements
BEGIN {
printf( "2" );
product = 2;
for( i = 2; i <= 8; i ++ )
{
nextV = product + 1;
# find the first prime factor of nextV
p = 3;
found = 0;
while( p * p <= nextV && ! ( found = nextV % p == 0 ) )
{
p += 2;
}
if( found )
{
nextV = p;
}
printf( " %d", nextV );
product *= nextV
}
}
- Output:
2 3 7 43 13 53 5 6221671
BASIC
Craft Basic
define size = 16, em = 0
dim list[size]
let list[0] = 2
print 2, " ",
for i = 1 to 15
let k = 3
do
let em = 1
for j = 0 to i - 1
let em = (em * list[j]) % k
next j
let em = (em + 1) % k
if em = 0 then
let list[i] = k
print list[i], " ",
break
endif
let k = k + 2
wait
loop
next i
FreeBASIC
Naive and takes forever to find the largest term, but does get there in the end.
dim as ulongint E(0 to 15), k
dim as integer i, em
E(0) = 2 : print 2
for i=1 to 15
k=3
do
em = 1
for j as uinteger = 0 to i-1
em = (em*E(j)) mod k
next j
em = (em + 1) mod k
if em = 0 then
E(i)=k
print E(i)
exit do
end if
k = k + 2
loop
next i
EasyLang
limit = 8
arr[] = [ 2 ]
write 2 & " "
for i = 2 to limit
k = 3
repeat
em = 1
for j = 1 to i - 1
em = em * arr[j] mod k
.
em = (em + 1) mod k
until em = 0
k += 2
.
arr[] &= k
write k & " "
.
F#
//Euclid-Mullin sequence. Nigel Galloway: October 29th., 2021
let(|Prime|_|)(n,g)=if Open.Numeric.Primes.MillerRabin.IsProbablePrime &g then Some(n*g,n*g+1I) else None
let n=Seq.unfold(fun(n,g)->match n,g with Prime n->Some(g,n) |_->let g=Open.Numeric.Primes.Extensions.PrimeExtensions.PrimeFactors g|>Seq.item 1 in Some(g,(n*g,n*g+1I)))(1I,2I)
n|>Seq.take 16|>Seq.iter(printfn "%A")
- Output:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Fermat
Func Firstfac(n) =
j := 3;
up := Sqrt(n);
while j <= up do
if Divides(j,n) then Return(j) fi;
j:=j+2;
od;
Return(n).;
Array eu[16];
eu[1]:=2;
!(eu[1],' ');
for i=2 to 16 do
eu[i]:=Firstfac(1+Prod<k=1,i-1>[eu[k]]);
!(eu[i],' ');
od;
- Output:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Go
This runs in about 54 seconds which, puzzlingly, is a good bit slower than Wren even though both are using GMP and the Pollard's rho algorithm. I have no idea why.
package main
import (
"fmt"
big "github.com/ncw/gmp"
"log"
)
var (
zero = big.NewInt(0)
one = big.NewInt(1)
two = big.NewInt(2)
three = big.NewInt(3)
four = big.NewInt(4)
five = big.NewInt(5)
six = big.NewInt(6)
ten = big.NewInt(10)
k100 = big.NewInt(100000)
)
func pollardRho(n, c *big.Int) *big.Int {
g := func(x, y *big.Int) *big.Int {
x2 := new(big.Int)
x2.Mul(x, x)
x2.Add(x2, c)
return x2.Mod(x2, y)
}
x, y, z := big.NewInt(2), big.NewInt(2), big.NewInt(1)
d := new(big.Int)
count := 0
for {
x = g(x, n)
y = g(g(y, n), n)
d.Sub(x, y)
d.Abs(d)
d.Mod(d, n)
z.Mul(z, d)
count++
if count == 100 {
d.GCD(nil, nil, z, n)
if d.Cmp(one) != 0 {
break
}
z.Set(one)
count = 0
}
}
if d.Cmp(n) == 0 {
return zero
}
return d
}
func smallestPrimeFactorWheel(n, max *big.Int) *big.Int {
if n.ProbablyPrime(15) {
return n
}
z := new(big.Int)
if z.Rem(n, two).Cmp(zero) == 0 {
return two
}
if z.Rem(n, three).Cmp(zero) == 0 {
return three
}
if z.Rem(n, five).Cmp(zero) == 0 {
return five
}
k := big.NewInt(7)
i := 0
inc := []*big.Int{four, two, four, two, four, six, two, six}
for z.Mul(k, k).Cmp(n) <= 0 {
if z.Rem(n, k).Cmp(zero) == 0 {
return k
}
k.Add(k, inc[i])
if k.Cmp(max) > 0 {
break
}
i = (i + 1) % 8
}
return nil
}
func smallestPrimeFactor(n *big.Int) *big.Int {
s := smallestPrimeFactorWheel(n, k100)
if s != nil {
return s
}
c := big.NewInt(1)
s = new(big.Int).Set(n)
for {
d := pollardRho(n, c)
if d.Cmp(zero) == 0 {
if c.Cmp(ten) == 0 {
log.Fatal("Pollard Rho doesn't appear to be working.")
}
c.Add(c, one)
} else {
// get the smallest prime factor of 'd'
factor := smallestPrimeFactorWheel(d, d)
// check whether n/d has a smaller prime factor
s = smallestPrimeFactorWheel(n.Quo(n, d), factor)
if s != nil {
if s.Cmp(factor) < 0 {
return s
} else {
return factor
}
} else {
return factor
}
}
}
}
func main() {
k := 19
fmt.Println("First", k, "terms of the Euclid–Mullin sequence:")
fmt.Println(2)
prod := big.NewInt(2)
z := new(big.Int)
count := 1
for count < k {
z.Add(prod, one)
t := smallestPrimeFactor(z)
fmt.Println(t)
prod.Mul(prod, t)
count++
}
}
- Output:
First 19 terms of the Euclid–Mullin sequence: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003 30693651606209 37 1741
J
2x (, <./@(>:&.(*/)))@[&_~ 15
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Java
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;
import java.util.concurrent.ThreadLocalRandom;
public final class EuclidMullinSequence {
public static void main(String[] aArgs) {
primes = listPrimesUpTo(1_000_000);
System.out.println("The first 27 terms of the Euclid-Mullin sequence:");
System.out.print(2 + " ");
for ( int i = 1; i < 27; i++ ) {
System.out.print(String.format("%s%s", nextEuclidMullin(), ( i == 14 || i == 27 ) ? "\n" : " "));
}
}
private static BigInteger nextEuclidMullin() {
BigInteger smallestPrime = smallestPrimeFactor(product.add(BigInteger.ONE));
product = product.multiply(smallestPrime);
return smallestPrime;
}
private static BigInteger smallestPrimeFactor(BigInteger aNumber) {
if ( aNumber.isProbablePrime(CERTAINTY_LEVEL) ) {
return aNumber;
}
for ( BigInteger prime : primes ) {
if ( aNumber.mod(prime).signum() == 0 ) {
return prime;
}
}
BigInteger factor = pollardsRho(aNumber);
return smallestPrimeFactor(factor);
}
private static BigInteger pollardsRho(BigInteger aN) {
if ( aN.equals(BigInteger.ONE) ) {
return BigInteger.ONE;
}
if ( aN.mod(BigInteger.TWO).signum() == 0 ) {
return BigInteger.TWO;
}
final BigInteger core = new BigInteger(aN.bitLength(), random);
BigInteger x = new BigInteger(aN.bitLength(), random);
BigInteger xx = x;
BigInteger divisor = null;
do {
x = x.multiply(x).mod(aN).add(core).mod(aN);
xx = xx.multiply(xx).mod(aN).add(core).mod(aN);
xx = xx.multiply(xx).mod(aN).add(core).mod(aN);
divisor = x.subtract(xx).gcd(aN);
} while ( divisor.equals(BigInteger.ONE) );
return divisor;
}
private static List<BigInteger> listPrimesUpTo(int aLimit) {
BitSet sieve = new BitSet(aLimit + 1);
sieve.set(2, aLimit + 1);
for ( int i = 2; i * i <= aLimit; i = sieve.nextSetBit(i + 1) ) {
for ( int j = i * i; j <= aLimit; j = j + i ) {
sieve.clear(j);
}
}
List<BigInteger> result = new ArrayList<BigInteger>(sieve.cardinality());
for ( int i = 2; i >= 0; i = sieve.nextSetBit(i + 1) ) {
result.add(BigInteger.valueOf(i));
}
return result;
}
private static List<BigInteger> primes;
private static BigInteger product = BigInteger.TWO;
private static ThreadLocalRandom random = ThreadLocalRandom.current();
private static final int CERTAINTY_LEVEL = 20;
}
- Output:
The first 27 terms of the Euclid-Mullin sequence: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003 30693651606209 37 1741 1313797957 887 71 7127 109 23 97 159227
jq
Works with jq, gojq and jaq, that is, the C, Go and Rust implementations of jq. (*)
Adapted from Algol 68
(*) The precision of jq and jaq is insufficient to compute more than the first nine numbers in the sequence (i.e., beyond 38709183810571). gojq's memory consumption becomes excessive after the first 16 numbers in the sequence are produced.
# Output: the Euclid-Mullins sequence, beginning with 2
def euclid_mullins:
foreach range(1; infinite|floor) as $i ( { product: 1 };
.next = .product + 1
# find the first prime factor of .next
| .p = 3
| .found = false
| until( .p * .p > .next or .found;
.found = ((.next % .p) == 0)
| if .found then . else .p += 2 end)
| if .found then .next = .p else . end
| .product *= .next)
| .next ;
# Produce 16 terms
limit(16; euclid_mullins)
- Output:
gojq supports unbounded-precision integer arithmetic and was accordingly used to produce the following output in accordance with the basic task requirements Beyond this number, gojq's memory consumption becomes excessive.
Invocation: $ gojq -n -f euclid-mullin-sequence.algol.jq
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Julia
using Primes
struct EuclidMullin end
Base.length(em::EuclidMullin) = 1000 # not expected to get to 1000
Base.eltype(em::EuclidMullin) = BigInt
Base.iterate(em::EuclidMullin, t=big"1") = (p = first(first(factor(t + 1).pe)); (p, t * p))
println("First 16 Euclid-Mullin numbers: ", join(Iterators.take(EuclidMullin(), 16), ", "))
- Output:
First 16 Euclid-Mullin numbers: 2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
Mathematica/Wolfram Language
list = {2};
Do[
prod = Times @@ list;
prod++;
new = Min[FactorInteger[prod][[All, 1]]];
AppendTo[list, new]
,
{21 - 1}
];
list
- Output:
The first 21 numbers of the sequence:
{2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003, 30693651606209, 37, 1741, 1313797957, 887}
Others may be found by adjusting the range of the Do loop but it will take a while.
Lua
-- find elements of the Euclid-Mullin sequence: starting from 2,
-- the next element is the smallest prime factor of 1 + the product
-- of the previous elements
do
io.write( "2" )
local product = 2
for i = 2, 8 do
local nextV = product + 1
-- find the first prime factor of nextV
local p = 3
local found = false
while p * p <= nextV and not found do
found = nextV % p == 0
if not found then p = p + 2 end
end
if found then nextV = p end
io.write( " ", nextV )
product = product * nextV
end
end
- Output:
2 3 7 43 13 53 5 6221671
Alternative using Pollard's Rho algorithm. Uses the iterative gcd function from the Greatest common divisor task.
.
Note that, as discussed on the Talk page, Pollard's Rho algorithm won't necessarily find the lowest factor, however it does for the first 16 elements.
As with the other Lua sample, only 8 elements are found due to the size of some of the subsequent ones.
function gcd(a,b)
while b~=0 do
a,b=b,a%b
end
return math.abs(a)
end
function pollard_rho(n)
local x, y, d = 2, 2, 1
local g = function(x) return (x*x+1) % n end
while d == 1 do
x = g(x)
y = g(g(y))
d = gcd(math.abs(x-y),n)
end
if d == n then return d end
return math.min(d, math.floor( n/d ) )
end
local ar, product = {2}, 2
repeat
ar[ #ar + 1 ] = pollard_rho( product + 1 )
product = product * ar[ #ar ]
until #ar >= 8
print( table.concat(ar, " ") )
- Output:
2 3 7 43 13 53 5 6221671
Maxima
euclid_mullin(n):=if n=1 then 2 else ifactors(1+product(euclid_mullin(i),i,1,n-1))[1][1]$
/* Test case */
makelist(euclid_mullin(k),k,16);
- Output:
[2,3,7,43,13,53,5,6221671,38709183810571,139,2801,11,17,5471,52662739,23003]
MiniScript
// find elements of the Euclid-Mullin sequence: starting from 2,
// the next element is the smallest prime factor of 1 + the product
// of the previous elements
seq = [2]
product = 2
for i in range( 2, 8 )
nextV = product + 1
// find the first prime factor of nextV
p = 3
found = false
while p * p <= nextV and not found
found = nextV % p == 0
if not found then p = p + 2
end while
if found then nextV = p
seq.push( nextV )
product = product * nextV
end for
print seq.join( " ")
- Output:
2 3 7 43 13 53 5 6221671
Nim
import integers
let
Zero = newInteger()
One = newInteger(1)
Two = newInteger(2)
Three = newInteger(3)
Five = newInteger(5)
Ten = newInteger(10)
Max = newInteger(100000)
None = newInteger(-1)
proc pollardRho(n, c: Integer): Integer =
template g(x: Integer): Integer = (x * x + c) mod n
var
x = newInteger(2)
y = newInteger(2)
z = newInteger(1)
d = Max + 1
count = 0
while true:
x = g(x)
y = g(g(y))
d = abs(x - y) mod n
z *= d
inc count
if count == 100:
d = gcd(z, n)
if d != One: break
z = newInteger(1)
count = 0
result = if d == n: Zero else: d
template isEven(n: Integer): bool = isZero(n and 1)
proc smallestPrimeFactorWheel(n: Integer): Integer =
if n.isPrime(5): return n
if n.isEven: return Two
if isZero(n mod 3): return Three
if isZero(n mod 5): return Five
var k = newInteger(7)
var i = 0
const Inc = [4, 2, 4, 2, 4, 6, 2, 6]
while k * k <= n:
if isZero(n mod k): return k
k += Inc[i]
if k > Max: return None
i = (i + 1) mod 8
proc smallestPrimeFactor(n: Integer): Integer =
var n = n
result = smallestPrimeFactorWheel(n)
if result != None: return
var c = One
result = newInteger(n)
while n > Max:
var d = pollardRho(n, c)
if d.isZero:
if c == Ten:
quit "Pollard Rho doesn't appear to be working.", QuitFailure
inc c
else:
# Can't be sure PR will find the smallest prime factor first.
result = min(result, d)
n = n div d
if n.isPrime(2):
return min(result, n)
proc main() =
var k = 19
echo "First ", k, " terms of the Euclid–Mullin sequence:"
echo 2
var prod = newInteger(2)
var count = 1
while count < k:
let t = smallestPrimeFactor(prod + One)
echo t
prod *= t
inc count
main()
- Output:
The first 16 numbers are displayed instantly, but it took about 9 seconds on my (moderately powerful) laptop to get the next three. I gave up for the 20th number.
First 19 terms of the Euclid–Mullin sequence: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003 30693651606209 37 1741
PARI/GP
E=vector(16)
E[1]=2
for(i=2,16,E[i]=factor(prod(n=1,i-1,E[n])+1)[1,1])
print(E)
- Output:
[2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003]
Perl
use strict;
use warnings;
use feature 'say';
use ntheory <factor vecprod vecmin>;
my @Euclid_Mullin = 2;
push @Euclid_Mullin, vecmin factor (1 + vecprod @Euclid_Mullin) for 2..16+11;
say "First sixteen: @Euclid_Mullin[ 0..15]";
say "Next eleven: @Euclid_Mullin[16..26]";
- Output:
First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003 Next eleven: 30693651606209 37 1741 1313797957 887 71 7127 109 23 97 159227
Phix
with javascript_semantics requires("1.0.1") -- (added mpz_set_v()) include mpfr.e sequence res = {} mpz {total,tmp} = mpz_inits(2,1) while length(res)<16 do mpz_add_si(tmp,total,1) mpz_set_v(tmp,mpz_pollard_rho(tmp)[1][1]) res = append(res,mpz_get_str(tmp)) mpz_mul(total,total,tmp) end while printf(1,"The first 16 Euclid-Mulin numbers: %s\n",{join(res)})
- Output:
The first 16 Euclid-Mulin numbers: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
While the first 16 are pretty fast, mpz_pollard_rho("723023114226131400979589798874734076807875188379971") took 3 minutes, and yielded the next element as 30693651606209, but beyond that I gave up.
Python
""" Rosetta code task: Euclid-Mullin_sequence """
from primePy import primes
def euclid_mullin():
""" generate Euclid-Mullin sequence """
total = 1
while True:
next_iter = primes.factor(total + 1)
total *= next_iter
yield next_iter
GEN = euclid_mullin()
print('First 16 Euclid-Mullin numbers:', ', '.join(str(next(GEN)) for _ in range(16)))
- Output:
First 16 Euclid-Mullin numbers: 2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
Raku
use Prime::Factor;
my @Euclid-Mullin = 2, { state $i = 1; (1 + [×] @Euclid-Mullin[^$i++]).&prime-factors.min } … *;
put 'First sixteen: ', @Euclid-Mullin[^16];
- Output:
First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Ring
// find elements of the Euclid-Mullin sequence: starting from 2,
// the next element is the smallest prime factor of 1 + the product
// of the previous elements
see "2"
product = 2
for i = 2 to 8
nextV = product + 1
// find the first prime factor of nextV
p = 3
found = false
while p * p <= nextV and not found
found = ( nextV % p ) = 0
if not found p = p + 2 ok
end
if found nextV = p ok
see " " + nextV
product = product * nextV
next
- Output:
2 3 7 43 13 53 5 6221671
RPL
RPL code | Comment |
---|---|
≪ IF # 1d DUP2 AND ≠ THEN DROP # 2d ELSE IF DUP 3 DUP2 / * == THEN DROP # 3d ELSE DUP B→R √ → divm ≪ 4 5 divm FOR n IF OVER n DUP2 / * == THEN SWAP DROP n R→B SWAP divm 'n' STO END 6 SWAP - DUP STEP DROP ≫ END END ≫ ‘bDIV1’ STO ≪ DUP SIZE 1 1 ROT FOR j OVER j GET * NEXT 1 + bDIV1 + ≫ 'NXTEM' STO |
bDIV1 ( #m -- #first_divisor ) is #2 a divisor ? is #3 a divisor ? otherwise get sqrt(m) d = 4 ; for n = 5 to sqrt(m) if n divides m replace m by n and prepare loop exit d = 6 - d ; n += d NXTEM ( { #EM(1) .. #EM(n) } -- { #EM(1) .. #EM(n+1) } ) get EM(1)*..*EM(n) get least prime factor of 1+EM(1)*..*EM(n) and add to list |
- Input:
≪ { # 2 } WHILE DUP SIZE ≤ 16 REPEAT NXTEM END ≫ EVAL
The emulator's watchdog timer prevents checking the primality of EM(9) = # 38709183810571d. Even if this device stayed idle, EM(10) could not be calculated, since the product of all earlier elements is more than 64 bits long.
- Output:
1: { # 2d # 3d # 7d # 43d # 13d # 53d # 5d # 6221671d }
Ruby
def pollard_rho(n)
x, y, d = 2, 2, 1
g = proc{|x|(x*x+1) % n}
while d == 1 do
x = g[x]
y = g[g[y]]
d = (x-y).abs.gcd(n)
end
return d if d == n
[d, n/d].compact.min
end
ar = [2]
ar << pollard_rho(ar.inject(&:*)+1) until ar.size >= 16
puts ar.join(", ")
- Output:
2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
SETL
program euclid_mullin;
print(2);
product := 2;
loop for i in [2..16] do
next := smallest_factor(product + 1);
product *:= next;
print(next);
end loop;
op smallest_factor(n);
if even n then return 2; end if;
d := 3;
loop while d*d <= n do
if n mod d=0 then return d; end if;
d +:= 2;
end loop;
return n;
end op;
end program;
- Output:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Sidef
func f(n) is cached {
return 2 if (n == 1)
lpf(1 + prod(1..^n, {|k| f(k) }))
}
say f.map(1..16)
say f.map(17..27)
- Output:
[2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003] [30693651606209, 37, 1741, 1313797957, 887, 71, 7127, 109, 23, 97, 159227]
Wren
Wren-cli
This uses the Pollard Rho algorithm to try and speed up the factorization of the 15th element but overall time still slow at around 32 seconds.
import "./big" for BigInt
var zero = BigInt.zero
var one = BigInt.one
var two = BigInt.two
var ten = BigInt.ten
var k100 = BigInt.new(100000)
var smallestPrimeFactorWheel = Fn.new { |n, max|
if (n.isProbablePrime(5)) return n
if (n % 2 == zero) return BigInt.two
if (n % 3 == zero) return BigInt.three
if (n % 5 == zero) return BigInt.five
var k = BigInt.new(7)
var i = 0
var inc = [4, 2, 4, 2, 4, 6, 2, 6]
while (k * k <= n) {
if (n % k == zero) return k
k = k + inc[i]
if (k > max) return null
i = (i + 1) % 8
}
}
var smallestPrimeFactor = Fn.new { |n|
var s = smallestPrimeFactorWheel.call(n, k100)
if (s) return s
var c = one
while (true) {
var d = BigInt.pollardRho(n, 2, c)
if (d == 0) {
if (c == ten) Fiber.abort("Pollard Rho doesn't appear to be working.")
c = c + one
} else {
// get the smallest prime factor of 'd'
var factor = smallestPrimeFactorWheel.call(d, d)
// check whether n/d has a smaller prime factor
s = smallestPrimeFactorWheel.call(n/d, factor)
return s ? BigInt.min(s, factor) : factor
}
}
}
var k = 16
System.print("First %(k) terms of the Euclid–Mullin sequence:")
System.print(2)
var prod = BigInt.two
var count = 1
while (count < k) {
var t = smallestPrimeFactor.call(prod + one)
System.print(t)
prod = prod * t
count = count + 1
}
- Output:
First 16 terms of the Euclid–Mullin sequence: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Embedded
This finds the first 16 in 0.11 seconds but the next 11 takes 6 minutes 17 seconds.
If we could assume that Pollard's Rho will always find the smallest prime factor (or a multiple thereof) first, then this would bring the runtime down to 44 seconds and still produce the correct answers for this particular task. However, in general it is not safe to assume that - see discussion in Talk Page.
/* Euclid_mullin_sequence_2.wren */
import "./gmp" for Mpz
var k100 = Mpz.from(100000)
var smallestPrimeFactorTrial = Fn.new { |n, max|
if (n.probPrime(15) > 0) return n
var k = Mpz.one
while (k * k <= n) {
k.nextPrime
if (k > max) return null
if (n.isDivisible(k)) return k
}
}
var smallestPrimeFactor = Fn.new { |n|
var s = smallestPrimeFactorTrial.call(n, k100)
if (s) return s
var c = Mpz.one
while (true) {
var d = Mpz.pollardRho(n, 2, c)
if (d.isZero) {
if (c == 100) Fiber.abort("Pollard Rho doesn't appear to be working.")
c.inc
} else {
// get the smallest prime factor of 'd'
var factor = smallestPrimeFactorTrial.call(d, d)
// check whether n/d has a smaller prime factor
s = smallestPrimeFactorTrial.call(n/d, factor)
return s ? Mpz.min(s, factor) : factor
}
}
}
var k = 27
System.print("First %(k) terms of the Euclid–Mullin sequence:")
System.print(2)
var prod = Mpz.two
var count = 1
while (count < k) {
var t = smallestPrimeFactor.call(prod + Mpz.one)
System.print(t)
prod.mul(t)
count = count + 1
}
- Output:
As Wren-cli plus eleven more:
30693651606209 37 1741 1313797957 887 71 7127 109 23 97 159227