Euclid-Mullin sequence: Difference between revisions
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
|||
Line 19:
=={{header|AWK}}==
<syntaxhighlight lang="awk">
# syntax: GAWK -f EUCLID-MULLIN_SEQUENCE.AWK
# converted from FreeBASIC
Line 45:
exit(0)
}
</syntaxhighlight>
{{out}}
<pre>
Line 51:
</pre>
=={{header|F_Sharp|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")
</syntaxhighlight>
{{out}}
<pre>
Line 78:
=={{header|Fermat}}==
<
j := 3;
up := Sqrt(n);
Line 94:
eu[i]:=Firstfac(1+Prod<k=1,i-1>[eu[k]]);
!(eu[i],' ');
od;</
{{out}}<pre> 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre>
=={{header|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
Line 118:
k = k + 2
loop
next i</
=={{header|Julia}}==
<
struct EuclidMullin end
Line 130:
println("First 16 Euclid-Mullin numbers: ", join(Iterators.take(EuclidMullin(), 16), ", "))
</
<pre>
First 16 Euclid-Mullin numbers: 2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
Line 136:
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
Do[
prod = Times @@ list;
Line 145:
{21 - 1}
];
list</
{{out}}
The first 21 numbers of the sequence:
Line 152:
=={{header|PARI/GP}}==
<
E[1]=2
for(i=2,16,E[i]=factor(prod(n=1,i-1,E[n])+1)[1,1])
print(E)</
{{out}}<pre>[2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003]</pre>
=={{header|Perl}}==
{{libheader|ntheory}}
<
use warnings;
use feature 'say';
Line 169:
say "First sixteen: @Euclid_Mullin[ 0..15]";
say "Next eleven: @Euclid_Mullin[16..26]";</
{{out}}
<pre>First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
Line 175:
=={{header|Phix}}==
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.1"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (added mpz_set_v())</span>
Line 189:
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The first 16 Euclid-Mulin numbers: %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)})</span>
<!--</
{{out}}
<pre>
Line 198:
=={{header|Python}}==
<
from primePy import primes
Line 212:
GEN = euclid_mullin()
print('First 16 Euclid-Mullin numbers:', ', '.join(str(next(GEN)) for _ in range(16)))
</
<pre>
First 16 Euclid-Mullin numbers: 2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
Line 219:
=={{header|Raku}}==
<syntaxhighlight lang="raku"
my @Euclid-Mullin = 2, { state $i = 1; (1 + [×] @Euclid-Mullin[^$i++]).&prime-factors.min } … *;
put 'First sixteen: ', @Euclid-Mullin[^16];</
{{out}}
<pre>First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre>
=={{header|Sidef}}==
<
return 2 if (n == 1)
lpf(1 + prod(1..^n, {|k| f(k) }))
Line 234:
say f.map(1..16)
say f.map(17..27)</
{{out}}
<pre>
Line 244:
===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.
<
var zero = BigInt.zero
Line 322:
prod = prod * t
count = count + 1
} </
{{out}}
Line 348:
{{libheader|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.
<
import "./gmp" for Mpz
Line 400:
prod.mul(t)
count = count + 1
}</
{{out}}
|
Revision as of 10:30, 27 August 2022
- 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
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
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
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
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
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
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