Jump to content

Euclid-Mullin sequence

From Rosetta Code
Task
Euclid-Mullin sequence
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

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

ALGOL W

Translation of: ALGOL 68 – but only showing the first 8 elements as Algol W integers are limited to 32-bit.
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.

Translation of: ALGOL 68
# 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

C++ 

#include <cstdint>
#include <iostream>

uint64_t product = 2;

uint64_t smallest_prime_factor(const uint64_t& number) {
	if ( number % 3 == 0 ) { return 3; }
	if ( number % 5 == 0 ) { return 5; }

	for ( uint64_t divisor = 7; divisor * divisor <= number; divisor += 2 ) {
		if ( number % divisor == 0 ) { return divisor; }
	}

	return number;
}

uint64_t next_euclid_mullin() {
	const uint64_t smallest_prime = smallest_prime_factor(product + 1);
	product *= smallest_prime;
	return smallest_prime;
}

int main() {
	std::cout << "The first 9 terms of the Euclid-Mullin sequence:" << "\n";
	std::cout << 2 << "  ";
	for ( uint32_t i = 1; i < 9; ++i ) {
		std::cout << next_euclid_mullin() << "  ";
	}
	std::cout << "\n";
}
Output:
The first 9 terms of the Euclid-Mullin sequence:
2  3  7  43  13  53  5  6221671  38709183810571

EasyLang

Translation of: AWK
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

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)
    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 += 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

Translation of: ALGOL W
-- 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.

Translation of: 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.
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

Translation of: 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
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

Translation of: Wren
Library: integers
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]

PascalABC.NET

Translation of: C++
function smallest_prime_factor(number: biginteger): biginteger;
begin
  var divisor := 3bi;
  repeat
    if (number mod divisor) = 0 then 
    begin
      result := divisor; 
      exit
    end;
    divisor += 2;
  until divisor * divisor > number;
  result := number;
end;

function euclid_mullin(): sequence of biginteger;
begin
  var product := 2bi;
  yield product;
  while true do
  begin
    var smallest_prime := smallest_prime_factor(product + 1);
    product *= smallest_prime;
    yield smallest_prime;
  end;
end;

begin
  euclid_mullin.take(16).println;
end.
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

REXX

Modules: How to use
Modules: Source code

The most direct way to solve this task turned out to use a procedure delivering the first n numbers, Euclids(x) in module Sequences. There procedure Ffactor() is used, first prime factor. The stretch task is out of REXX' reach. 

include Settings

say 'EUCLID-MULLIN SEQUENCE - 4 Mar 2025'
say version
say
numeric digits 100
call Euclids 16
call Task 16
say Format(Time('e'),,3) 'seconds'
exit

Task:
arg xx
say 'The first' xx 'Euclid-Mullin numbers are:'
do i = 1 to xx
   call charout ,eucl.prime.i' '
end
say
return

include Functions
include Sequences
include Numbers
include Abend
Output:
EUCLID-MULLIN SEQUENCE - 4 Mar 2025
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022

The first 16 Euclid-Mullin numbers are:
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
34.837 seconds

Ring

Translation of: 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
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

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 }

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

Library: Wren-gmp

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
Retrieved from "https://rosettacode.org/wiki/Euclid-Mullin_sequence?oldid=377993"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information