I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Euclid-Mullin sequence

Euclid-Mullin sequence 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.
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.

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

## AWK

` # syntax: GAWK -f EUCLID-MULLIN_SEQUENCE.AWK# converted from FreeBASICBEGIN {    limit = 7 # we'll stop here    arr = 2    printf("%s ",arr)    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
```

## F#

` //Euclid-Mullin sequence. Nigel Galloway: October 29th., 2021let(|Prime|_|)(n,g)=if Open.Numeric.Primes.MillerRabin.IsProbablePrime &g then Some(n*g,n*g+1I) else Nonelet 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;eu:=2;!(eu,' ');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`

## FreeBASIC

Naive and takes forever to find the largest term, but does get there in the end.

` dim as ulongint E(0 to 15), kdim as integer i, emE(0) = 2 : print 2for 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    loopnext i`

## Julia

`using Primes struct EuclidMullin end Base.length(em::EuclidMullin) = 1000  # not expected to get to 1000Base.eltype(em::EuclidMullin) = BigIntBase.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
```

## PARI/GP

`E=vector(16)E=2for(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

Library: ntheory
`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
include mpfr.e

sequence res = {}
mpz {total,tmp} = mpz_inits(2,1)
while length(res)<16 do
mpz_set_v(tmp,mpz_pollard_rho(tmp))
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`

## 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.zerovar one  = BigInt.onevar two  = BigInt.twovar ten  = BigInt.tenvar 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 = 16System.print("First %(k) terms of the Euclid–Mullin sequence:")System.print(2)var prod = BigInt.twovar count = 1while (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

Library: Wren-gmp

This finds the first 16 in 0.11 seconds and the next 3 in around 39 seconds. I gave up after that as it would take too long for the Pollard's Rho algorithm to find any more.

`/* euclid_mullin_gmp.wren */ import "./gmp" for Mpz var max = Mpz.from(100000) var smallestPrimeFactorWheel = Fn.new { |n|    if (n.probPrime(15) > 0) return n    if (n.isEven) return Mpz.two    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) {        if (n.isDivisible(k)) return k        k.add(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 = Mpz.one    s = n.copy()    while (n > max) {        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 {            // can't be sure PR will find the smallest prime factor first            s.min(d)            n.div(d)            if (n.probPrime(5) > 0) return Mpz.min(s, n)        }    }    return s} var k = 19System.print("First %(k) terms of the Euclid–Mullin sequence:")System.print(2)var prod = Mpz.twovar count = 1while (count < k) {    var t = smallestPrimeFactor.call(prod + Mpz.one)    System.print(t)    prod.mul(t)    count = count + 1}`
Output:

As Wren-cli plus three more:

```30693651606209
37
1741
```