Euclid-Mullin sequence: Difference between revisions
m
→{{header|Wren}}: Minor tidy
m (→{{header|Wren}}: Minor tidy) |
|||
| (28 intermediate revisions by 12 users not shown) | |||
Line 46:
<pre>
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003
</pre>
=={{header|ALGOL W}}==
{{Trans|ALGOL 68|but only showing the first 8 elements as Algol W integers are limited to 32-bit.}}
<syntaxhighlight lang="algolw">
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.
</syntaxhighlight>
{{out}}
<pre>
2 3 7 43 13 53 5 6221671
</pre>
Line 81 ⟶ 112:
</pre>
Alternative version.
{{Trans|ALGOL 68}}
<syntaxhighlight lang="awk">
# 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
}
}
</syntaxhighlight>
{{out}}
<pre>
2 3 7 43 13 53 5 6221671
</pre>
=={{header|BASIC}}==
==={{header|Craft Basic}}===
<syntaxhighlight lang="basic">define size = 16, em = 0
dim list[size]
let list[0] = 2
print 2, " ",
for i = 1 to 15
Line 98 ⟶ 164:
for j = 0 to i - 1
let em = (
next j
let em = (
if em = 0 then
let list[i] = k
print list[i], " ",
break
Line 118 ⟶ 184:
loop
next i</syntaxhighlight>
==={{header|FreeBASIC}}===
Naive and takes forever to find the largest term, but does get there in the end.
<syntaxhighlight 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</syntaxhighlight>
=={{header|EasyLang}}==
{{trans|AWK}}
<syntaxhighlight>
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 & " "
.
</syntaxhighlight>
=={{header|F_Sharp|F#}}==
Line 170 ⟶ 277:
od;</syntaxhighlight>
{{out}}<pre> 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre>
=={{header|Go}}==
Line 215 ⟶ 299:
six = big.NewInt(6)
ten = big.NewInt(10)
)
Line 251 ⟶ 335:
}
func smallestPrimeFactorWheel(n, max *big.Int) *big.Int {
if n.ProbablyPrime(15) {
return n
Line 282 ⟶ 366:
func smallestPrimeFactor(n *big.Int) *big.Int {
s := smallestPrimeFactorWheel(n, k100)
if s != nil {
return s
Line 288 ⟶ 372:
c := big.NewInt(1)
s = new(big.Int).Set(n)
for
d := pollardRho(n, c)
if d.Cmp(zero) == 0 {
Line 296 ⟶ 380:
c.Add(c, one)
} else {
//
// check whether
s = smallestPrimeFactorWheel(n.Quo(n, d), factor)
if
} else
return factor
}
} else
return factor
}
}
}
}
Line 351 ⟶ 436:
1741
</pre>
=={{header|J}}==
<syntaxhighlight lang="j"> 2x (, <./@(>:&.(*/)))@[&_~ 15
2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</syntaxhighlight>
=={{header|Java}}==
<syntaxhighlight lang="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 = 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;
}
</syntaxhighlight>
{{ out }}
<pre>
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
</pre>
=={{header|jq}}==
'''Works with jq, gojq and jaq, that is, the C, Go and Rust implementations of jq.''' (*)
'''Adapted from [[#Algol_68|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.
<syntaxhighlight lang=jq>
# 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)
</syntaxhighlight>
{{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
<pre>
2
3
7
43
13
53
5
6221671
38709183810571
139
2801
11
17
5471
52662739
23003
</pre>
=={{header|Julia}}==
<syntaxhighlight lang="julia">using Primes
Line 384 ⟶ 622:
<pre>{2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003, 30693651606209, 37, 1741, 1313797957, 887}</pre>
Others may be found by adjusting the range of the Do loop but it will take a while.
=={{header|Lua}}==
{{Trans|ALGOL W}}
<syntaxhighlight lang="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
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
</syntaxhighlight>
{{out}}
<pre>
2 3 7 43 13 53 5 6221671
</pre>
Alternative using Pollard's Rho algorithm. Uses the iterative gcd function from the [[Greatest common divisor]] task.<br>
{{Trans|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.<br>
As with the other Lua sample, only 8 elements are found due to the size of some of the subsequent ones.
<syntaxhighlight lang="lua">
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, " ") )
</syntaxhighlight>
{{out}}
<pre>
2 3 7 43 13 53 5 6221671
</pre>
=={{header|Maxima}}==
<syntaxhighlight lang="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);
</syntaxhighlight>
{{out}}
<pre>
[2,3,7,43,13,53,5,6221671,38709183810571,139,2801,11,17,5471,52662739,23003]
</pre>
=={{header|MiniScript}}==
{{Trans|Lua}}
<syntaxhighlight lang="miniscript">
// 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( " ")
</syntaxhighlight>
{{out}}
<pre>
2 3 7 43 13 53 5 6221671
</pre>
=={{header|Nim}}==
{{trans|Wren}}
{{libheader|integers}}
<syntaxhighlight lang="Nim">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()
</syntaxhighlight>
{{out}}
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.
<pre>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
</pre>
=={{header|PARI/GP}}==
Line 460 ⟶ 913:
{{out}}
<pre>First sixteen: 2 3 7 43 13 53 5 6221671 38709183810571 139 2801 11 17 5471 52662739 23003</pre>
=={{header|Ring}}==
{{Trans|Lua}}
<syntaxhighlight lang="ring">
// 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
</syntaxhighlight>
{{out}}
<pre>
2 3 7 43 13 53 5 6221671
</pre>
=={{header|RPL}}==
Line 509 ⟶ 989:
1: { # 2d # 3d # 7d # 43d # 13d # 53d # 5d # 6221671d }
</pre>
=={{header|Ruby}}==
<syntaxhighlight lang="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(", ")
</syntaxhighlight>
{{out}}
<pre>2, 3, 7, 43, 13, 53, 5, 6221671, 38709183810571, 139, 2801, 11, 17, 5471, 52662739, 23003
</pre>
=={{header|SETL}}==
<syntaxhighlight lang="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;</syntaxhighlight>
{{out}}
<pre>2
3
7
43
13
53
5
6221671
38709183810571
139
2801
11
17
5471
52662739
23003</pre>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">func f(n) is cached {
Line 526 ⟶ 1,064:
===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.
<syntaxhighlight lang="
var zero = BigInt.zero
Line 532 ⟶ 1,070:
var two = BigInt.two
var ten = BigInt.ten
var
var smallestPrimeFactorWheel = Fn.new { |n, max|
if (n.isProbablePrime(5)) return n
if (n % 2 == zero) return BigInt.two
Line 575 ⟶ 1,089:
var smallestPrimeFactor = Fn.new { |n|
var s = smallestPrimeFactorWheel.call(n, k100)
if (s) return s
var c = one
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 {
//
return s ? BigInt.min(s, factor) : factor
}
}
}
Line 604 ⟶ 1,117:
prod = prod * t
count = count + 1
}
{{out}}
Line 627 ⟶ 1,140:
</pre>
<br>
===Embedded===
{{libheader|Wren-gmp}}
This finds the first 16 in 0.11 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.
<syntaxhighlight lang="wren">/* Euclid_mullin_sequence_2.wren */
import "./gmp" for Mpz
var
var
if (n.probPrime(15) > 0) return n
while (k * k <= n) {
if (k > max) return null
}
}
var smallestPrimeFactor = Fn.new { |n|
var s =
if (s) return s
var c = Mpz.one
var d = Mpz.pollardRho(n, 2, c)
if (d.isZero) {
Line 663 ⟶ 1,172:
c.inc
} else {
//
// check whether n
return s ? Mpz.min(s, factor) : factor
}
}
}
var k =
System.print("First %(k) terms of the Euclid–Mullin sequence:")
System.print(2)
Line 685 ⟶ 1,194:
{{out}}
As Wren-cli plus
<pre>
30693651606209
37
1741
1313797957
887
71
7127
109
23
97
159227
</pre>
| |||