Euclid-Mullin sequence: Difference between revisions

Content added Content deleted

Inline

Revision as of 10:19, 8 April 2023

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

OEIS sequence A000945

ALGOL 68

Works with: ALGOL 68G version 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.
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

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

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

print "done."

end

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

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

Go

Translation of: Wren

Library: GMP(Go wrapper)

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)
    max   = 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 *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)
    if s != nil {
        return s
    }
    c := big.NewInt(1)
    s = new(big.Int).Set(n)
    for n.Cmp(max) > 0 {
        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 {
            // can't be sure PR will find the smallest prime factor first
            if d.Cmp(s) < 0 {
                s.Set(d)
            }
            n.Quo(n, d)
            if n.ProbablyPrime(5) {
                if n.Cmp(s) < 0 {
                    return n
                }
                return s
            }
        }
    }
    return s
}

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

Java

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.

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

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

RPL

Works with: Halcyon Calc version 4.2.7

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 }

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

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 = 19
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 three more:

30693651606209
37
1741

@@ Line 217: / Line 217: @@
     max   = big.NewInt(100000)
 )
-=={{header|Java}}==
 func pollardRho(n, c *big.Int) *big.Int {
@@ Line 353: / Line 351: @@
 </pre>
+=={{header|Java}}==
 =={{header|Julia}}==