Euclid-Mullin sequence
- 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
F#
<lang fsharp> //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") </lang>
- Output:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Fermat
<lang 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;</lang>
- Output:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
FreeBASIC
Naive and takes forever to find the largest term, but does get there in the end. <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</lang>
Julia
<lang 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), ", "))
</lang>
- Output:
First 16 Euclid-Mullin numbers: 2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
PARI/GP
<lang parigp>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)</lang>
- Output:
[2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003]
Perl
<lang 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]";</lang>
- 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
Python
<lang python>""" Rosetta code task: Euclid-Mullin_sequence """
from primePy import primes
def euler_mullin():
""" generate Euler-Mullin sequence """ total = 1 while True: next_iter = primes.factor(total + 1) total *= next_iter yield next_iter
GEN = euler_mullin() print('First 16 Euler-Mullin numbers:', ', '.join(str(next(GEN)) for _ in range(16)))
</lang>
- Output:
First 16 Euler-Mullin numbers: 2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
Raku
<lang perl6>use Prime::Factor;
my @Euclid-Mullin = 2, { state $i = 1; (1 + [×] @Euclid-Mullin[^$i++]).&prime-factors.min } … *;
put 'First sixteen: ', @Euclid-Mullin[^16];</lang>
- Output:
First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Wren
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. <lang ecmascript>import "./big" for BigInt
var zero = BigInt.zero var one = BigInt.one var two = BigInt.two var ten = BigInt.ten var max = BigInt.new(100000)
var pollardRho = Fn.new { |n, c|
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 % 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) if (s) return s var c = one s = n while (n > max) { var d = pollardRho.call(n, c) if (d == 0) { if (c == ten) Fiber.abort("Pollard Rho doesn't appear to be working.") c = c + one } else { // can't be sure PR will find the smallest prime factor first s = BigInt.min(s, d) n = n / d if (n.isProbablePrime(2)) return BigInt.min(s, n) } } return s
}
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
} </lang>
- 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