Factors of a Mersenne number: Difference between revisions
(→{{header|Julia}}: Corrected the solution in two important ways and optimised it a bit) |
(Added Easylang) |
||
(51 intermediate revisions by 19 users not shown) | |||
Line 57: | Line 57: | ||
;See also: |
;See also: |
||
* [https://www.youtube.com/watch?v=SNwvJ7psoow Computers in 1948: 2 |
* [https://www.youtube.com/watch?v=SNwvJ7psoow Computers in 1948: 2<sup>127</sup> - 1] <br> (Note: This video is no longer available because the YouTube account associated with this video has been terminated.) |
||
<br><br> |
<br><br> |
||
=={{header|11l}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="11l">F is_prime(a) |
|||
I a == 2 {R 1B} |
|||
I a < 2 | a % 2 == 0 {R 0B} |
|||
L(i) (3 .. Int(sqrt(a))).step(2) |
|||
I a % i == 0 |
|||
R 0B |
|||
R 1B |
|||
F m_factor(p) |
|||
V max_k = 16384 I/ p |
|||
L(k) 0 .< max_k |
|||
V q = 2 * p * k + 1 |
|||
I !is_prime(q) |
|||
L.continue |
|||
E I q % 8 != 1 & q % 8 != 7 |
|||
L.continue |
|||
E I pow(2, p, q) == 1 |
|||
R q |
|||
R 0 |
|||
V exponent = Int(input(‘Enter exponent of Mersenne number: ’)) |
|||
I !is_prime(exponent) |
|||
print(‘Exponent is not prime: #.’.format(exponent)) |
|||
E |
|||
V factor = m_factor(exponent) |
|||
I factor == 0 |
|||
print(‘No factor found for M#.’.format(exponent)) |
|||
E |
|||
print(‘M#. has a factor: #.’.format(exponent, factor))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Enter exponent of Mersenne number: 929 |
|||
M929 has a factor: 13007 |
|||
</pre> |
|||
=={{header|8086 Assembly}}== |
|||
<syntaxhighlight lang="asm">P: equ 929 ; P for 2^P-1 |
|||
cpu 8086 |
|||
bits 16 |
|||
org 100h |
|||
section .text |
|||
mov ax,P ; Is P prime? |
|||
call prime |
|||
mov dx,notprm |
|||
jc msg ; If not, say so and stop. |
|||
xor bp,bp ; Let BP hold k |
|||
test_k: inc bp ; k += 1 |
|||
mov ax,P ; Calculate 2kP + 1 |
|||
mul bp ; AX = kP |
|||
shl ax,1 ; AX = 2kP |
|||
inc ax ; AX = 2kP + 1 |
|||
mov dx,ovfl ; If AX overflows (16 bits), say so and stop |
|||
jc msg |
|||
mov bx,ax ; What is 2^P mod (2kP+1)? |
|||
mov cx,P |
|||
call modpow |
|||
dec ax ; If it is 1, we're done |
|||
jnz test_k ; If not, increment K and try again |
|||
mov dx,factor ; If so, we found a factor |
|||
call msg |
|||
prbx: mov ax,10 ; The factor is still in BX |
|||
xchg bx,ax ; Put factor in AX and divisor (10) in BX |
|||
mov di,number ; Generate ASCII representation of number |
|||
digit: xor dx,dx |
|||
div bx ; Divide current number by 10, |
|||
add dl,'0' ; add '0' to remainder, |
|||
dec di ; move pointer back, |
|||
mov [di],dl ; store digit, |
|||
test ax,ax ; and if there are more digits, |
|||
jnz digit ; find the next digit. |
|||
mov dx,di ; Finally, print the number. |
|||
jmp msg |
|||
;;; Calculate 2^CX mod BX |
|||
;;; Output: AX |
|||
;;; Destroyed: CX, DX |
|||
modpow: shl cx,1 ; Shift CX left until top bit in high bit |
|||
jnc modpow ; Keep shifting while carry zero |
|||
rcr cx,1 ; Put the top bit back into CX |
|||
mov ax,1 ; Start with square = 1 |
|||
.step: mul ax ; Square (result is 32-bit, goes in DX:AX) |
|||
shl cx,1 ; Grab a bit from CX |
|||
jnc .nodbl ; If zero, don't multiply by two |
|||
shl ax,1 ; Shift DX:AX left (mul by two) |
|||
rcl dx,1 |
|||
.nodbl: div bx ; Divide DX:AX by BX (putting modulus in DX) |
|||
mov ax,dx ; Continue with modulus |
|||
jcxz .done ; When CX reaches 0, we're done |
|||
jmp .step ; Otherwise, do the next step |
|||
.done: ret |
|||
;;; Is AX prime? |
|||
;;; Output: carry clear if prime, set if not prime. |
|||
;;; Destroyed: AX, BX, CX, DX, SI, DI, BP |
|||
prime: mov cx,[prcnt] ; See if AX is already in the list of primes |
|||
mov di,primes |
|||
repne scasw ; If so, return |
|||
je modpow.done ; Reuse the RET just above here (carry clear) |
|||
mov bp,ax ; Move AX out of the way |
|||
mov bx,[di-2] ; Start generating new primes |
|||
.sieve: inc bx ; BX = last prime + 2 |
|||
inc bx |
|||
cmp bp,bx ; If BX higher than number to test, |
|||
jb modpow.done ; stop, number is not prime. (carry set) |
|||
mov cx,[prcnt] ; CX = amount of current primes |
|||
mov si,primes ; SI = start of primes |
|||
.try: mov ax,bx ; BX divisible by current prime? |
|||
xor dx,dx |
|||
div word [si] |
|||
test dx,dx ; If so, BX is not prime. |
|||
jz .sieve |
|||
inc si |
|||
inc si |
|||
loop .try ; Otherwise, try next prime. |
|||
mov ax,bx ; If we get here, BX _is_ prime |
|||
stosw ; So add it to the list |
|||
inc word [prcnt] ; We have one more prime |
|||
cmp ax,bp ; Is it the prime we are looking for? |
|||
jne .sieve ; If not, try next prime |
|||
ret |
|||
;;; Print message in DX |
|||
msg: mov ah,9 |
|||
int 21h |
|||
ret |
|||
section .data |
|||
db "*****" ; Placeholder for number |
|||
number: db "$" |
|||
notprm: db "P is not prime.$" |
|||
ovfl: db "Range of factor exceeded (max 16 bits)." |
|||
factor: db "Found factor: $" |
|||
prcnt: dw 2 ; Amount of primes currently in list |
|||
primes: dw 2, 3 ; List of primes to be extended</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Found factor: 13007</pre> |
|||
=={{header|360 Assembly}}== |
=={{header|360 Assembly}}== |
||
Line 64: | Line 204: | ||
Use of bitwise operations |
Use of bitwise operations |
||
(TM (Test under Mask), SLA (Shift Left Arithmetic),SRA (Shift Right Arithmetic)). |
(TM (Test under Mask), SLA (Shift Left Arithmetic),SRA (Shift Right Arithmetic)). |
||
<lang>* Factors of a Mersenne number 11/09/2015 |
<syntaxhighlight lang="text">* Factors of a Mersenne number 11/09/2015 |
||
MERSENNE CSECT |
MERSENNE CSECT |
||
USING MERSENNE,R15 |
USING MERSENNE,R15 |
||
Line 153: | Line 293: | ||
PG DS CL24 buffer |
PG DS CL24 buffer |
||
YREGS |
YREGS |
||
END MERSENNE</ |
END MERSENNE</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 161: | Line 301: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
mersenne.adb: |
mersenne.adb: |
||
< |
<syntaxhighlight lang="ada">with Ada.Text_IO; |
||
-- reuse Is_Prime from [[Primality by Trial Division]] |
-- reuse Is_Prime from [[Primality by Trial Division]] |
||
with Is_Prime; |
with Is_Prime; |
||
Line 228: | Line 368: | ||
Ada.Text_IO.Put_Line ("is not a Mersenne number"); |
Ada.Text_IO.Put_Line ("is not a Mersenne number"); |
||
end; |
end; |
||
end Mersenne;</ |
end Mersenne;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 240: | Line 380: | ||
<!-- {{works with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386}} |
<!-- {{works with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386}} |
||
Compiles, but I couldn't maxint not in library, works with manually entered maxint, bits width. Leaving some issue with newline --> |
Compiles, but I couldn't maxint not in library, works with manually entered maxint, bits width. Leaving some issue with newline --> |
||
< |
<syntaxhighlight lang="algol68">MODE ISPRIMEINT = INT; |
||
PR READ "prelude/is_prime.a68" PR; |
PR READ "prelude/is_prime.a68" PR; |
||
Line 290: | Line 430: | ||
FI |
FI |
||
END</ |
END</syntaxhighlight> |
||
Example: |
Example: |
||
<pre> |
<pre> |
||
Line 296: | Line 436: | ||
M +929 has a factor: +13007 |
M +929 has a factor: +13007 |
||
</pre> |
</pre> |
||
=={{header|Arturo}}== |
|||
<syntaxhighlight lang="rebol">mersenneFactors: function [q][ |
|||
if not? prime? q -> print "number not prime!" |
|||
r: new q |
|||
while -> r > 0 |
|||
-> shl 'r 1 |
|||
d: new 1 + 2 * q |
|||
while [true][ |
|||
i: new 1 |
|||
p: new r |
|||
while [p <> 0][ |
|||
i: new (i * i) % d |
|||
if p < 0 -> 'i * 2 |
|||
if i > d -> 'i - d |
|||
shl 'p 1 |
|||
] |
|||
if? i <> 1 -> 'd + 2 * q |
|||
else -> break |
|||
] |
|||
print ["2 ^" q "- 1 = 0 ( mod" d ")"] |
|||
] |
|||
mersenneFactors 929</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 ^ 929 - 1 = 0 ( mod 13007 )</pre> |
|||
=={{header|AutoHotkey}}== |
=={{header|AutoHotkey}}== |
||
ahk [http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=144 discussion] |
ahk [http://www.autohotkey.com/forum/viewtopic.php?t=44657&postdays=0&postorder=asc&start=144 discussion] |
||
< |
<syntaxhighlight lang="autohotkey">MsgBox % MFact(27) ;-1: 27 is not prime |
||
MsgBox % MFact(2) ; 0 |
MsgBox % MFact(2) ; 0 |
||
MsgBox % MFact(3) ; 0 |
MsgBox % MFact(3) ; 0 |
||
Line 354: | Line 523: | ||
y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1 |
y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1 |
||
Return y |
Return y |
||
}</ |
}</syntaxhighlight> |
||
=={{header|BBC BASIC}}== |
=={{header|BBC BASIC}}== |
||
< |
<syntaxhighlight lang="bbcbasic"> PRINT "A factor of M929 is "; FNmersenne_factor(929) |
||
PRINT "A factor of M937 is "; FNmersenne_factor(937) |
PRINT "A factor of M937 is "; FNmersenne_factor(937) |
||
END |
END |
||
Line 394: | Line 563: | ||
ENDWHILE |
ENDWHILE |
||
= Y% |
= Y% |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre>A factor of M929 is 13007 |
<pre>A factor of M929 is 13007 |
||
Line 400: | Line 569: | ||
=={{header|Bracmat}}== |
=={{header|Bracmat}}== |
||
< |
<syntaxhighlight lang="bracmat">( ( modPow |
||
= square P divisor highbit log 2pow |
= square P divisor highbit log 2pow |
||
. !arg:(?P.?divisor) |
. !arg:(?P.?divisor) |
||
Line 461: | Line 630: | ||
| out$"no divisors found" |
| out$"no divisors found" |
||
) |
) |
||
);</ |
);</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>found some divisors of 2^!P-1 : 13007 and 348890248924938259750454781163390930305120269538278042934009621348894657205785 |
<pre>found some divisors of 2^!P-1 : 13007 and 348890248924938259750454781163390930305120269538278042934009621348894657205785 |
||
Line 469: | Line 638: | ||
=={{header|C}}== |
=={{header|C}}== |
||
< |
<syntaxhighlight lang="c">int isPrime(int n){ |
||
if (n%2==0) return n==2; |
if (n%2==0) return n==2; |
||
if (n%3==0) return n==3; |
if (n%3==0) return n==3; |
||
Line 492: | Line 661: | ||
else break; |
else break; |
||
} while(1); |
} while(1); |
||
printf("2^%d - 1 = 0 (mod %d)\n", q, d);}</ |
printf("2^%d - 1 = 0 (mod %d)\n", q, d);}</syntaxhighlight> |
||
=={{header|C sharp}}== |
=={{header|C sharp|C#}}== |
||
< |
<syntaxhighlight lang="csharp">using System; |
||
namespace prog |
namespace prog |
||
Line 540: | Line 709: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
=={{header|C++}}== |
|||
<syntaxhighlight lang="cpp">#include <iostream> |
|||
#include <cstdint> |
|||
typedef uint64_t integer; |
|||
integer bit_count(integer n) { |
|||
integer count = 0; |
|||
for (; n > 0; count++) |
|||
n >>= 1; |
|||
return count; |
|||
} |
|||
integer mod_pow(integer p, integer n) { |
|||
integer square = 1; |
|||
for (integer bits = bit_count(p); bits > 0; square %= n) { |
|||
square *= square; |
|||
if (p & (1 << --bits)) |
|||
square <<= 1; |
|||
} |
|||
return square; |
|||
} |
|||
bool is_prime(integer n) { |
|||
if (n < 2) |
|||
return false; |
|||
if (n % 2 == 0) |
|||
return n == 2; |
|||
for (integer p = 3; p * p <= n; p += 2) |
|||
if (n % p == 0) |
|||
return false; |
|||
return true; |
|||
} |
|||
integer find_mersenne_factor(integer p) { |
|||
for (integer k = 0, q = 1;;) { |
|||
q = 2 * ++k * p + 1; |
|||
if ((q % 8 == 1 || q % 8 == 7) && mod_pow(p, q) == 1 && is_prime(q)) |
|||
return q; |
|||
} |
|||
return 0; |
|||
} |
|||
int main() { |
|||
std::cout << find_mersenne_factor(929) << std::endl; |
|||
return 0; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
13007 |
|||
</pre> |
|||
=={{header|Clojure}}== |
=={{header|Clojure}}== |
||
{{trans|Python}} |
{{trans|Python}} |
||
< |
<syntaxhighlight lang="lisp">(ns mersennenumber |
||
(:gen-class)) |
(:gen-class)) |
||
Line 603: | Line 825: | ||
:let [s (-main p)]] |
:let [s (-main p)]] |
||
(println s)) |
(println s)) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Output}} |
{{Output}} |
||
<pre> |
<pre> |
||
Line 631: | Line 853: | ||
{{trans|Ruby}} |
{{trans|Ruby}} |
||
< |
<syntaxhighlight lang="coffeescript">mersenneFactor = (p) -> |
||
limit = Math.sqrt(Math.pow(2,p) - 1) |
limit = Math.sqrt(Math.pow(2,p) - 1) |
||
k = 1 |
k = 1 |
||
Line 658: | Line 880: | ||
checkMersenne(prime) for prime in ["2","3","4","5","7","11","13","17","19","23","29","31","37","41","43","47","53","929"] |
checkMersenne(prime) for prime in ["2","3","4","5","7","11","13","17","19","23","29","31","37","41","43","47","53","929"] |
||
</syntaxhighlight> |
|||
</lang> |
|||
<pre>M2 = 2^2-1 is prime |
<pre>M2 = 2^2-1 is prime |
||
Line 682: | Line 904: | ||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
||
{{trans|Maxima}} |
{{trans|Maxima}} |
||
< |
<syntaxhighlight lang="lisp">(defun mersenne-fac (p &aux (m (1- (expt 2 p)))) |
||
(loop for k from 1 |
(loop for k from 1 |
||
for n = (1+ (* 2 k p)) |
for n = (1+ (* 2 k p)) |
||
Line 688: | Line 910: | ||
finally (return n))) |
finally (return n))) |
||
(print (mersenne-fac 929))</ |
(print (mersenne-fac 929))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 697: | Line 919: | ||
We can use a primality test from the [[Primality by Trial Division#Common_Lisp|Primality by Trial Division]] task. |
We can use a primality test from the [[Primality by Trial Division#Common_Lisp|Primality by Trial Division]] task. |
||
< |
<syntaxhighlight lang="lisp">(defun primep (n) |
||
"Is N prime?" |
"Is N prime?" |
||
(and (> n 1) |
(and (> n 1) |
||
(or (= n 2) (oddp n)) |
(or (= n 2) (oddp n)) |
||
(loop for i from 3 to (isqrt n) by 2 |
(loop for i from 3 to (isqrt n) by 2 |
||
never (zerop (rem n i)))))</ |
never (zerop (rem n i)))))</syntaxhighlight> |
||
Specific to this task, we define modulo-power and mersenne-prime-p. |
Specific to this task, we define modulo-power and mersenne-prime-p. |
||
< |
<syntaxhighlight lang="lisp">(defun modulo-power (base power modulus) |
||
(loop with square = 1 |
(loop with square = 1 |
||
for bit across (format nil "~b" power) |
for bit across (format nil "~b" power) |
||
Line 724: | Line 946: | ||
(primep q) |
(primep q) |
||
(= 1 (modulo-power 2 power q))) |
(= 1 (modulo-power 2 power q))) |
||
(return (values nil q)))))</ |
(return (values nil q)))))</syntaxhighlight> |
||
We can run the same tests from the [[#Ruby|Ruby]] entry. |
We can run the same tests from the [[#Ruby|Ruby]] entry. |
||
Line 751: | Line 973: | ||
M53 = 2**53-1 is composite with factor 6361. |
M53 = 2**53-1 is composite with factor 6361. |
||
M929 = 2**929-1 is composite with factor 13007.</pre> |
M929 = 2**929-1 is composite with factor 13007.</pre> |
||
=={{header|Crystal}}== |
|||
{{trans|Ruby}} |
|||
<syntaxhighlight lang="ruby">require "big" |
|||
def prime?(n) # P3 Prime Generator primality test |
|||
return n | 1 == 3 if n < 5 # n: 0,1,4|false, 2,3|true |
|||
return false if n.gcd(6) != 1 # for n a P3 prime candidate (pc) |
|||
pc1, pc2 = -1, 1 # use P3's prime candidates sequence |
|||
until (pc1 += 6) > Math.sqrt(n).to_i # pcs are only 1/3 of all integers |
|||
return false if n % pc1 == 0 || n % (pc2 += 6) == 0 # if n is composite |
|||
end |
|||
true |
|||
end |
|||
# Compute b**e mod m |
|||
def powmod(b, e, m) |
|||
r, b = 1.to_big_i, b.to_big_i |
|||
while e > 0 |
|||
r = (r * b) % m if e.odd? |
|||
b = (b * b) % m |
|||
e >>= 1 |
|||
end |
|||
r |
|||
end |
|||
def mersenne_factor(p) |
|||
mers_num = 2.to_big_i ** p - 1 |
|||
kp2 = p2 = 2.to_big_i * p |
|||
while (kp2 - 1) ** 2 < mers_num |
|||
q = kp2 + 1 # return q if it's a factor |
|||
return q if [1, 7].includes?(q % 8) && prime?(q) && (powmod(2, p, q) == 1) |
|||
kp2 += p2 |
|||
end |
|||
true # could also set to `0` value to check for |
|||
end |
|||
def check_mersenne(p) |
|||
print "M#{p} = 2**#{p}-1 is " |
|||
f = mersenne_factor(p) |
|||
(puts "prime"; return) if f.is_a?(Bool) # or f == 0 |
|||
puts "composite with factor #{f}" |
|||
end |
|||
(2..53).each { |p| check_mersenne(p) if prime?(p) } |
|||
check_mersenne 929</syntaxhighlight> |
|||
{{out}} |
|||
<pre>M2 = 2**2-1 is prime |
|||
M3 = 2**3-1 is prime |
|||
M5 = 2**5-1 is prime |
|||
M7 = 2**7-1 is prime |
|||
M11 = 2**11-1 is composite with factor 23 |
|||
M13 = 2**13-1 is prime |
|||
M17 = 2**17-1 is prime |
|||
M19 = 2**19-1 is prime |
|||
M23 = 2**23-1 is composite with factor 47 |
|||
M29 = 2**29-1 is composite with factor 233 |
|||
M31 = 2**31-1 is prime |
|||
M37 = 2**37-1 is composite with factor 223 |
|||
M41 = 2**41-1 is composite with factor 13367 |
|||
M43 = 2**43-1 is composite with factor 431 |
|||
M47 = 2**47-1 is composite with factor 2351 |
|||
M53 = 2**53-1 is composite with factor 6361 |
|||
M929 = 2**929-1 is composite with factor 13007</pre> |
|||
=={{header|D}}== |
=={{header|D}}== |
||
< |
<syntaxhighlight lang="d">import std.stdio, std.math, std.traits; |
||
ulong mersenneFactor(in ulong p) pure nothrow @nogc { |
ulong mersenneFactor(in ulong p) pure nothrow @nogc { |
||
Line 789: | Line 1,076: | ||
void main() { |
void main() { |
||
writefln("Factor of M929: %d", 929.mersenneFactor); |
writefln("Factor of M929: %d", 929.mersenneFactor); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>Factor of M929: 13007</pre> |
<pre>Factor of M929: 13007</pre> |
||
=={{header|Delphi}}== |
|||
See [https://rosettacode.org/wiki/Factors_of_a_Mersenne_number#Pascal Pascal]. |
|||
=={{header|EasyLang}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight> |
|||
fastfunc isprim num . |
|||
i = 2 |
|||
while i <= sqrt num |
|||
if num mod i = 0 |
|||
return 0 |
|||
. |
|||
i += 1 |
|||
. |
|||
return 1 |
|||
. |
|||
func bit_count n . |
|||
while n > 0 |
|||
n = bitshift n -1 |
|||
cnt += 1 |
|||
. |
|||
return cnt |
|||
. |
|||
func mod_pow p n . |
|||
square = 1 |
|||
bits = bit_count p |
|||
while bits > 0 |
|||
square *= square |
|||
bits -= 1 |
|||
if bitand p bitshift 1 bits > 0 |
|||
square = bitshift square 1 |
|||
. |
|||
square = square mod n |
|||
. |
|||
return square |
|||
. |
|||
func mersenne_factor p . |
|||
while 1 = 1 |
|||
k += 1 |
|||
q = 2 * k * p + 1 |
|||
if (q mod 8 = 1 or q mod 8 = 7) and mod_pow p q = 1 and isprim q = 1 |
|||
return q |
|||
. |
|||
. |
|||
. |
|||
print mersenne_factor 929 |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
13007 |
|||
</pre> |
|||
=={{header|EchoLisp}}== |
=={{header|EchoLisp}}== |
||
< |
<syntaxhighlight lang="scheme"> |
||
;; M = 2^P - 1 , P prime |
;; M = 2^P - 1 , P prime |
||
;; look for a prime divisor q such as : q < √ M, q = 1 or 7 modulo 8, q = 1 + 2kP |
;; look for a prime divisor q such as : q < √ M, q = 1 or 7 modulo 8, q = 1 + 2kP |
||
Line 818: | Line 1,157: | ||
→ #t |
→ #t |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Elixir}}== |
=={{header|Elixir}}== |
||
{{trans|Ruby}} |
{{trans|Ruby}} |
||
< |
<syntaxhighlight lang="elixir">defmodule Mersenne do |
||
def mersenne_factor(p) do |
def mersenne_factor(p) do |
||
limit = :math.sqrt(:math.pow(2, p) - 1) |
limit = :math.sqrt(:math.pow(2, p) - 1) |
||
Line 857: | Line 1,196: | ||
[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,929] |
[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,929] |
||
|> Enum.each(fn p -> Mersenne.check_mersenne(p) end)</ |
|> Enum.each(fn p -> Mersenne.check_mersenne(p) end)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 884: | Line 1,223: | ||
The modpow function is not my original. This is a translation of python, more or less. |
The modpow function is not my original. This is a translation of python, more or less. |
||
< |
<syntaxhighlight lang="erlang"> |
||
-module(mersene2). |
-module(mersene2). |
||
-export([prime/1,modpow/3,mf/1]). |
-export([prime/1,modpow/3,mf/1]). |
||
Line 927: | Line 1,266: | ||
false -> divisors(N, C-1) |
false -> divisors(N, C-1) |
||
end. |
end. |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 942: | Line 1,281: | ||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
< |
<syntaxhighlight lang="factor">USING: combinators.short-circuit interpolate io kernel locals |
||
math math.bits math.functions math.primes sequences ; |
math math.bits math.functions math.primes sequences ; |
||
IN: rosetta-code.mersenne-factors |
IN: rosetta-code.mersenne-factors |
||
Line 973: | Line 1,312: | ||
[ [I No factor found for M${}.I] ] if* nl ; |
[ [I No factor found for M${}.I] ] if* nl ; |
||
929 test-mersenne</ |
929 test-mersenne</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 980: | Line 1,319: | ||
=={{header|Forth}}== |
=={{header|Forth}}== |
||
< |
<syntaxhighlight lang="forth">: prime? ( odd -- ? ) |
||
3 |
3 |
||
begin 2dup dup * >= |
begin 2dup dup * >= |
||
Line 1,012: | Line 1,351: | ||
929 factor-mersenne . \ 13007 |
929 factor-mersenne . \ 13007 |
||
4423 factor-mersenne . \ 0</ |
4423 factor-mersenne . \ 0</syntaxhighlight> |
||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
{{works with|Fortran|90 and later}} |
{{works with|Fortran|90 and later}} |
||
< |
<syntaxhighlight lang="fortran">PROGRAM EXAMPLE |
||
IMPLICIT NONE |
IMPLICIT NONE |
||
INTEGER :: exponent, factor |
INTEGER :: exponent, factor |
||
Line 1,067: | Line 1,406: | ||
Mfactor = 0 |
Mfactor = 0 |
||
END FUNCTION |
END FUNCTION |
||
END PROGRAM EXAMPLE</ |
END PROGRAM EXAMPLE</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
M929 has a factor: 13007 |
M929 has a factor: 13007 |
||
Line 1,073: | Line 1,412: | ||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="freebasic">' FB 1.05.0 Win64 |
||
Function isPrime(n As Integer) As Boolean |
Function isPrime(n As Integer) As Boolean |
||
Line 1,117: | Line 1,456: | ||
Print |
Print |
||
Print "Press any key to quit" |
Print "Press any key to quit" |
||
Sleep</ |
Sleep</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,137: | Line 1,476: | ||
2^97 - 1 = 0 (mod 11447) |
2^97 - 1 = 0 (mod 11447) |
||
2^929 - 1 = 0 (mod 13007) |
2^929 - 1 = 0 (mod 13007) |
||
</pre> |
|||
=={{header|Frink}}== |
|||
Frink has built-in routines for iterating through prime numbers and modular exponentiation. The following program will find all of the factors of the number given enough runtime. |
|||
<syntaxhighlight lang="frink">for p = primes[] |
|||
if modPow[2, 929, p] - 1 == 0 |
|||
println[p]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
13007 |
|||
</pre> |
</pre> |
||
=={{header|GAP}}== |
=={{header|GAP}}== |
||
< |
<syntaxhighlight lang="gap">MersenneSmallFactor := function(n) |
||
local k, m, d; |
local k, m, d; |
||
if IsPrime(n) then |
if IsPrime(n) then |
||
Line 1,171: | Line 1,520: | ||
FactorsInt(2^101-1); |
FactorsInt(2^101-1); |
||
# [ 7432339208719, 341117531003194129 ]</ |
# [ 7432339208719, 341117531003194129 ]</syntaxhighlight> |
||
=={{header|Go}}== |
=={{header|Go}}== |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 1,254: | Line 1,603: | ||
} |
} |
||
return int32(pow) |
return int32(pow) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,266: | Line 1,615: | ||
Using David Amos module Primes [https://web.archive.org/web/20121130222921/http://www.polyomino.f2s.com/david/haskell/codeindex.html] for prime number testing: |
Using David Amos module Primes [https://web.archive.org/web/20121130222921/http://www.polyomino.f2s.com/david/haskell/codeindex.html] for prime number testing: |
||
< |
<syntaxhighlight lang="haskell">import Data.List |
||
import HFM.Primes (isPrime) |
import HFM.Primes (isPrime) |
||
import Control.Monad |
import Control.Monad |
||
Line 1,278: | Line 1,627: | ||
map (succ.(2*m*)). enumFromTo 1 $ m `div` 2 |
map (succ.(2*m*)). enumFromTo 1 $ m `div` 2 |
||
bs = int2bin m |
bs = int2bin m |
||
pm n b = 2^b*n*n</ |
pm n b = 2^b*n*n</syntaxhighlight> |
||
< |
<syntaxhighlight lang="haskell">*Main> trialfac 929 |
||
[13007]</ |
[13007]</syntaxhighlight> |
||
=={{header|Icon}} and {{header|Unicon}}== |
=={{header|Icon}} and {{header|Unicon}}== |
||
Line 1,287: | Line 1,636: | ||
The following works in both languages: |
The following works in both languages: |
||
< |
<syntaxhighlight lang="unicon">procedure main(A) |
||
p := integer(A[1]) | 929 |
p := integer(A[1]) | 929 |
||
write("M",p," has a factor ",mfactor(p)) |
write("M",p," has a factor ",mfactor(p)) |
||
Line 1,318: | Line 1,667: | ||
} |
} |
||
return |
return |
||
end</ |
end</syntaxhighlight> |
||
Sample runs: |
Sample runs: |
||
Line 1,331: | Line 1,680: | ||
=={{header|J}}== |
=={{header|J}}== |
||
< |
<syntaxhighlight lang="j">trialfac=: 3 : 0 |
||
qs=. (#~8&(1=|+.7=|))(#~1&p:)1+(*(1x+i.@<:@<.)&.-:)y |
qs=. (#~8&(1=|+.7=|))(#~1&p:)1+(*(1x+i.@<:@<.)&.-:)y |
||
qs#~1=qs&|@(2&^@[**:@])/ 1,~ |.#: y |
qs#~1=qs&|@(2&^@[**:@])/ 1,~ |.#: y |
||
)</ |
)</syntaxhighlight> |
||
{{out|Examples}} |
{{out|Examples}} |
||
< |
<syntaxhighlight lang="j">trialfac 929 |
||
13007</ |
13007</syntaxhighlight> |
||
<lang |
<syntaxhighlight lang="j">trialfac 44497</syntaxhighlight> |
||
Empty output --> No factors found. |
Empty output --> No factors found. |
||
=={{header|Java}}== |
=={{header|Java}}== |
||
< |
<syntaxhighlight lang="java"> |
||
import java.math.BigInteger; |
import java.math.BigInteger; |
||
Line 1,412: | Line 1,761: | ||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 1,439: | Line 1,788: | ||
=={{header|JavaScript}}== |
=={{header|JavaScript}}== |
||
< |
<syntaxhighlight lang="javascript">function mersenne_factor(p){ |
||
var limit, k, q |
var limit, k, q |
||
limit = Math.sqrt(Math.pow(2,p) - 1) |
limit = Math.sqrt(Math.pow(2,p) - 1) |
||
Line 1,479: | Line 1,828: | ||
f = mersenne_factor(p) |
f = mersenne_factor(p) |
||
console.log(f == null ? "prime" : "composite with factor "+f) |
console.log(f == null ? "prime" : "composite with factor "+f) |
||
}</ |
}</syntaxhighlight> |
||
<pre> |
<pre> |
||
Line 1,491: | Line 1,840: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
< |
<syntaxhighlight lang="julia"># v0.6 |
||
using Primes |
using Primes |
||
Line 1,511: | Line 1,860: | ||
if mf != -1 println("M$i = ", mf, " × ", (big(2) ^ i - 1) ÷ mf) |
if mf != -1 println("M$i = ", mf, " × ", (big(2) ^ i - 1) ÷ mf) |
||
else println("M$i is prime") end |
else println("M$i is prime") end |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,529: | Line 1,878: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="scala">// version 1.0.6 |
||
fun isPrime(n: Int): Boolean { |
fun isPrime(n: Int): Boolean { |
||
Line 1,575: | Line 1,924: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,598: | Line 1,947: | ||
=={{header|Lingo}}== |
=={{header|Lingo}}== |
||
< |
<syntaxhighlight lang="lingo">on modPow (b, e, m) |
||
bits = getBits(e) |
bits = getBits(e) |
||
sq = 1 |
sq = 1 |
||
Line 1,620: | Line 1,969: | ||
end repeat |
end repeat |
||
return bits |
return bits |
||
end</ |
end</syntaxhighlight> |
||
< |
<syntaxhighlight lang="lingo">repeat with i = 2 to the maxInteger |
||
if modPow(2, 929, i)=1 then |
if modPow(2, 929, i)=1 then |
||
put "M929 has a factor: " & i |
put "M929 has a factor: " & i |
||
exit repeat |
exit repeat |
||
end if |
end if |
||
end repeat</ |
end repeat</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,634: | Line 1,983: | ||
</pre> |
</pre> |
||
=={{header|Mathematica}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
Believe it or not, this type of test runs faster in Mathematica than the squaring version described above. |
Believe it or not, this type of test runs faster in Mathematica than the squaring version described above. |
||
<syntaxhighlight lang="mathematica">For[i = 2, i < Prime[1000000], i = NextPrime[i], |
|||
<lang mathematica> |
|||
For[i = 2, i < Prime[1000000], i = NextPrime[i], |
|||
If[Mod[2^44497, i] == 1, |
If[Mod[2^44497, i] == 1, |
||
Print["divisible by "<>i]]]; Print["prime test passed; call Lucas and Lehmer"]</ |
Print["divisible by "<>i]]]; Print["prime test passed; call Lucas and Lehmer"]</syntaxhighlight> |
||
=={{header|Maxima}}== |
=={{header|Maxima}}== |
||
< |
<syntaxhighlight lang="maxima">mersenne_fac(p) := block([m: 2^p - 1, k: 1], |
||
while mod(m, 2 * k * p + 1) # 0 do k: k + 1, |
while mod(m, 2 * k * p + 1) # 0 do k: k + 1, |
||
2 * k * p + 1 |
2 * k * p + 1 |
||
Line 1,649: | Line 1,997: | ||
mersenne_fac(929); |
mersenne_fac(929); |
||
/* 13007 */</ |
/* 13007 */</syntaxhighlight> |
||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
{{trans|C}} |
{{trans|C}} |
||
< |
<syntaxhighlight lang="nim">import math |
||
proc isPrime(a: int): bool = |
proc isPrime(a: int): bool = |
||
Line 1,678: | Line 2,026: | ||
if i != 1: d += 2 * q |
if i != 1: d += 2 * q |
||
else: break |
else: break |
||
echo "2^",q," - 1 = 0 (mod ",d,")"</ |
echo "2^",q," - 1 = 0 (mod ",d,")"</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>2^929 - 1 = 0 (mod 13007)</pre> |
<pre>2^929 - 1 = 0 (mod 13007)</pre> |
||
Line 1,687: | Line 2,035: | ||
(GNU Octave has a <code>isprime</code> built-in test) |
(GNU Octave has a <code>isprime</code> built-in test) |
||
< |
<syntaxhighlight lang="octave">% test a bit; lsb is 1 (like built-in bit* ops) |
||
function b = bittst(n, p) |
function b = bittst(n, p) |
||
b = bitand(n, 2^(p-1)) > 0; |
b = bitand(n, 2^(p-1)) > 0; |
||
Line 1,720: | Line 2,068: | ||
endfunction |
endfunction |
||
printf("%d\n", Mfactor(929));</ |
printf("%d\n", Mfactor(929));</syntaxhighlight> |
||
=={{header|PARI/GP}}== |
=={{header|PARI/GP}}== |
||
This version takes about 15 microseconds to find a factor of 2<sup>929</sup> − 1. |
This version takes about 15 microseconds to find a factor of 2<sup>929</sup> − 1. |
||
< |
<syntaxhighlight lang="parigp">factorMersenne(p)={ |
||
forstep(q=2*p+1,sqrt(2)<<(p\2),2*p, |
forstep(q=2*p+1,sqrt(2)<<(p\2),2*p, |
||
[1,0,0,0,0,0,1][q%8] && Mod(2, q)^p==1 && return(q) |
[1,0,0,0,0,0,1][q%8] && Mod(2, q)^p==1 && return(q) |
||
Line 1,730: | Line 2,078: | ||
1<<p-1 |
1<<p-1 |
||
}; |
}; |
||
factorMersenne(929)</ |
factorMersenne(929)</syntaxhighlight> |
||
This implementation seems to be broken: |
This implementation seems to be broken: |
||
< |
<syntaxhighlight lang="parigp">TM(p) = local(status=1, i=1, len=0, S=0);{ |
||
printp("Test TM \t..."); |
printp("Test TM \t..."); |
||
S=2*p+1; |
S=2*p+1; |
||
Line 1,752: | Line 2,100: | ||
); |
); |
||
return(status); |
return(status); |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
{{trans|Fortran}} |
{{trans|Fortran}} |
||
< |
<syntaxhighlight lang="pascal">program FactorsMersenneNumber(input, output); |
||
function isPrime(n: longint): boolean; |
function isPrime(n: longint): boolean; |
||
Line 1,839: | Line 2,187: | ||
else |
else |
||
writeln('M', exponent, ' (2**', exponent, ' - 1) has the factor: ', factor); |
writeln('M', exponent, ' (2**', exponent, ' - 1) has the factor: ', factor); |
||
end.</ |
end.</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,848: | Line 2,196: | ||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use utf8; |
use utf8; |
||
Line 1,908: | Line 2,256: | ||
print $f? "M$m = $x = $q × @{[$x / $q]}\n" |
print $f? "M$m = $x = $q × @{[$x / $q]}\n" |
||
: "M$m = $x is prime\n"; |
: "M$m = $x is prime\n"; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>M2 = 3 is prime |
<pre>M2 = 3 is prime |
||
Line 1,923: | Line 2,271: | ||
Following the task introduction, this uses GMP's modular exponentiation and simple probable prime test for the small numbers, then looks for small factors before doing a Lucas-Lehmer test. For ranges above about p=2000, looking for small factors this way saves time (the amount of testing should be adjusted based on the input size and platform -- this example just uses a fixed amount). Note as well that the Lucas-Lehmer test shown here is ignoring the large speedup we can get by optimizing the modulo operation, but that's a different task. |
Following the task introduction, this uses GMP's modular exponentiation and simple probable prime test for the small numbers, then looks for small factors before doing a Lucas-Lehmer test. For ranges above about p=2000, looking for small factors this way saves time (the amount of testing should be adjusted based on the input size and platform -- this example just uses a fixed amount). Note as well that the Lucas-Lehmer test shown here is ignoring the large speedup we can get by optimizing the modulo operation, but that's a different task. |
||
< |
<syntaxhighlight lang="perl">use Math::GMP; |
||
# Use GMP's simple probable prime test. |
# Use GMP's simple probable prime test. |
||
Line 1,952: | Line 2,300: | ||
print "M$p is ", is_mersenne_prime($p) ? "prime" : "composite", "\n"; |
print "M$p is ", is_mersenne_prime($p) ? "prime" : "composite", "\n"; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,982: | Line 2,330: | ||
M929 = 13007 * 348890248924[.....]64184273 |
M929 = 13007 * 348890248924[.....]64184273 |
||
</pre> |
</pre> |
||
=={{header|Perl 6}}== |
|||
{{Works with|rakudo|2015.12}} |
|||
<lang perl6>my @primes = 2, 3, -> $n is copy { |
|||
repeat { $n += 2 } until $n %% none do for @primes -> $p { |
|||
last if $p > sqrt($n); |
|||
$p; |
|||
} |
|||
$n; |
|||
} ... *; |
|||
multi factors(1) { 1 } |
|||
multi factors(Int $remainder is copy) { |
|||
gather for @primes -> $factor { |
|||
if $factor * $factor > $remainder { |
|||
take $remainder if $remainder > 1; |
|||
last; |
|||
} |
|||
while $remainder %% $factor { |
|||
take $factor; |
|||
last if ($remainder div= $factor) === 1; |
|||
} |
|||
} |
|||
} |
|||
sub is_prime($x) { (state %){$x} //= factors($x) == 1 } |
|||
sub mtest($bits, $p) { |
|||
my @bits = $bits.base(2).comb; |
|||
loop (my $sq = 1; @bits; $sq %= $p) { |
|||
$sq *= $sq; |
|||
$sq += $sq if 1 == @bits.shift; |
|||
} |
|||
$sq == 1; |
|||
} |
|||
for flat 2 .. 60, 929 -> $m { |
|||
next unless is_prime($m); |
|||
my $f = 0; |
|||
my $x = 2**$m - 1; |
|||
my $q; |
|||
for 1..* -> $k { |
|||
$q = 2 * $k * $m + 1; |
|||
next unless $q % 8 == 1|7 or is_prime($q); |
|||
last if $q * $q > $x or $f = mtest($m, $q); |
|||
} |
|||
say $f ?? "M$m = $x\n\t= $q × { $x div $q }" |
|||
!! "M$m = $x is prime"; |
|||
}</lang> |
|||
{{out}} |
|||
<pre>M2 = 3 is prime |
|||
M3 = 7 is prime |
|||
M5 = 31 is prime |
|||
M7 = 127 is prime |
|||
M11 = 2047 |
|||
= 23 × 89 |
|||
M13 = 8191 is prime |
|||
M17 = 131071 is prime |
|||
M19 = 524287 is prime |
|||
M23 = 8388607 |
|||
= 47 × 178481 |
|||
M29 = 536870911 |
|||
= 233 × 2304167 |
|||
M31 = 2147483647 is prime |
|||
M37 = 137438953471 |
|||
= 223 × 616318177 |
|||
M41 = 2199023255551 |
|||
= 13367 × 164511353 |
|||
M43 = 8796093022207 |
|||
= 431 × 20408568497 |
|||
M47 = 140737488355327 |
|||
= 2351 × 59862819377 |
|||
M53 = 9007199254740991 |
|||
= 6361 × 1416003655831 |
|||
M59 = 576460752303423487 |
|||
= 179951 × 3203431780337 |
|||
M929 = 4538015467766671944574165338592225830478699345884382504442663144885072806275648112625635725391102144133907238129251016389326737199538896813326509341743147661691195191795226666084858428449394948944821764472508048114220424520501343042471615418544488778723282182172070046459244838911 |
|||
= 13007 × 348890248924938259750454781163390930305120269538278042934009621348894657205785201247454118966026150852149399410259938217062100192168747352450719561908445272675574320888385228421992652298715687625495638077382028762529439880103124705348782610789919949159935587158612289264184273</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Translation/Amalgamation of BBC BASIC, D, and Go |
Translation/Amalgamation of BBC BASIC, D, and Go |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>function is_prime(integer n) |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
if n<2 then return 0 end if |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">modpow</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
if n=2 then return 1 end if |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> |
|||
if remainder(n,2)=0 then return 0 end if |
|||
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> |
|||
for i=3 to floor(sqrt(n)) by 2 do |
|||
<span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span> |
|||
if remainder(n,i)=0 then |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">i</span> <span style="color: #008080;">do</span> |
|||
return 0 |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">and_bits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> |
|||
end if |
|||
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
return 1 |
|||
<span style="color: #000000;">z</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
function modpow(atom x, atom n, atom m) |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">y</span> |
|||
atom i = n, |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
y = 1, |
|||
z = x |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">mersenne_factor</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> |
|||
while i do |
|||
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #7060A8;">is_prime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
if and_bits(i,1) then |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">limit</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">))-</span><span style="color: #000000;">1</span> |
|||
y = mod(y*z,m) |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> |
|||
end if |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
|||
z = mod(z*z,m) |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">p</span><span style="color: #0000FF;">*</span><span style="color: #000000;">k</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span> |
|||
i = floor(i/2) |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">q</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">limit</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
end while |
|||
<span style="color: #008080;">if</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">8</span><span style="color: #0000FF;">),{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span><span style="color: #0000FF;">})</span> |
|||
return y |
|||
<span style="color: #008080;">and</span> <span style="color: #7060A8;">is_prime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">and</span> <span style="color: #000000;">modpow</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">q</span> |
|||
function mersenne_factor(integer p) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
if not is_prime(p) then return -1 end if |
|||
<span style="color: #000000;">k</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
|||
atom limit = sqrt(power(2,p))-1 |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
integer k = 1 |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">0</span> |
|||
while 1 do |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
atom q = 2*p*k + 1 |
|||
if q>=limit then exit end if |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">11</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">23</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">29</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">37</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">41</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">43</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">47</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">53</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">59</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">67</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">71</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">73</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">79</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">83</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">97</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">929</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">937</span><span style="color: #0000FF;">}</span> |
|||
if find(mod(q,8),{1,7}) |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
and is_prime(q) |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">ti</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> |
|||
and modpow(2,p,q)=1 then |
|||
<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;">"A factor of M%d is %d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">ti</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mersenne_factor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ti</span><span style="color: #0000FF;">)})</span> |
|||
return q |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
end if |
|||
<!--</syntaxhighlight>--> |
|||
k += 1 |
|||
end while |
|||
return 0 |
|||
end function |
|||
sequence tests = {11, 23, 29, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 97, 929, 937} |
|||
for i=1 to length(tests) do |
|||
integer ti = tests[i] |
|||
printf(1,"A factor of M%d is %d\n",{ti,mersenne_factor(ti)}) |
|||
end for</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 2,136: | Line 2,396: | ||
{{trans|D}} |
{{trans|D}} |
||
Requires bcmath |
Requires bcmath |
||
< |
<syntaxhighlight lang="php">echo 'M929 has a factor: ', mersenneFactor(929), '</br>'; |
||
function mersenneFactor($p) { |
function mersenneFactor($p) { |
||
Line 2,157: | Line 2,417: | ||
} |
} |
||
return true; |
return true; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,163: | Line 2,423: | ||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
< |
<syntaxhighlight lang="picolisp">(de **Mod (X Y N) |
||
(let M 1 |
(let M 1 |
||
(loop |
(loop |
||
Line 2,193: | Line 2,453: | ||
(prime? Q) |
(prime? Q) |
||
(= 1 (**Mod 2 P Q)) ) |
(= 1 (**Mod 2 P Q)) ) |
||
Q ) ) ) )</ |
Q ) ) ) )</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>: (for P (2 3 4 5 7 11 13 17 19 23 29 31 37 41 43 47 53 929) |
<pre>: (for P (2 3 4 5 7 11 13 17 19 23 29 31 37 41 43 47 53 929) |
||
Line 2,220: | Line 2,480: | ||
M53 = 2**53-1 is composite with factor 6361 |
M53 = 2**53-1 is composite with factor 6361 |
||
M929 = 2**929-1 is composite with factor 13007</pre> |
M929 = 2**929-1 is composite with factor 13007</pre> |
||
=={{header|Prolog}}== |
|||
<syntaxhighlight lang="prolog"> |
|||
mersenne_factor(P, F) :- |
|||
prime(P), |
|||
once(( |
|||
between(1, 100_000, K), % Fail if we can't find a small factor |
|||
Q is 2*K*P + 1, |
|||
test_factor(Q, P, F))). |
|||
test_factor(Q, P, prime) :- Q*Q > (1 << P - 1), !. |
|||
test_factor(Q, P, Q) :- |
|||
R is Q /\ 7, member(R, [1, 7]), |
|||
prime(Q), |
|||
powm(2, P, Q) =:= 1. |
|||
wheel235(L) :- |
|||
W = [4, 2, 4, 2, 4, 6, 2, 6 | W], |
|||
L = [1, 2, 2 | W]. |
|||
prime(N) :- |
|||
N >= 2, |
|||
wheel235(W), |
|||
prime(N, 2, W). |
|||
prime(N, D, _) :- D*D > N, !. |
|||
prime(N, D, [A|As]) :- |
|||
N mod D =\= 0, |
|||
D2 is D + A, prime(N, D2, As). |
|||
</syntaxhighlight> |
|||
{{Out}} |
|||
<pre> |
|||
?- mersenne_factor(23, X). |
|||
X = 47. |
|||
?- mersenne_factor(5,X). |
|||
X = prime. |
|||
?- mersenne_factor(25,X). |
|||
false. |
|||
?- mersenne_factor(929,X). |
|||
X = 13007. |
|||
?- mersenne_factor(127,X). |
|||
false. |
|||
</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
||
< |
<syntaxhighlight lang="python">def is_prime(number): |
||
return True # code omitted - see Primality by Trial Division |
return True # code omitted - see Primality by Trial Division |
||
Line 2,247: | Line 2,555: | ||
print "No factor found for M%d" % exponent |
print "No factor found for M%d" % exponent |
||
else: |
else: |
||
print "M%d has a factor: %d" % (exponent, factor)</ |
print "M%d has a factor: %d" % (exponent, factor)</syntaxhighlight> |
||
{{out|Example}} |
{{out|Example}} |
||
Line 2,256: | Line 2,564: | ||
=={{header|Racket}}== |
=={{header|Racket}}== |
||
< |
<syntaxhighlight lang="racket"> |
||
#lang racket |
#lang racket |
||
Line 2,277: | Line 2,585: | ||
(mersenne-factor 929) |
(mersenne-factor 929) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
< |
<syntaxhighlight lang="racket"> |
||
13007 |
13007 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
|||
(formerly Perl 6) |
|||
<syntaxhighlight lang="raku" line>sub mtest($bits, $p) { |
|||
my @bits = $bits.base(2).comb; |
|||
loop (my $sq = 1; @bits; $sq %= $p) { |
|||
$sq ×= $sq; |
|||
$sq += $sq if 1 == @bits.shift; |
|||
} |
|||
$sq == 1; |
|||
} |
|||
for flat 2 .. 60, 929 -> $m { |
|||
next unless is-prime($m); |
|||
my $f = 0; |
|||
my $x = 2**$m - 1; |
|||
my $q; |
|||
for 1..* -> $k { |
|||
$q = 2 × $k × $m + 1; |
|||
next unless $q % 8 == 1|7 or is-prime($q); |
|||
last if $q × $q > $x or $f = mtest($m, $q); |
|||
} |
|||
say $f ?? "M$m = $x\n\t= $q × { $x div $q }" |
|||
!! "M$m = $x is prime"; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>M2 = 3 is prime |
|||
M3 = 7 is prime |
|||
M5 = 31 is prime |
|||
M7 = 127 is prime |
|||
M11 = 2047 |
|||
= 23 × 89 |
|||
M13 = 8191 is prime |
|||
M17 = 131071 is prime |
|||
M19 = 524287 is prime |
|||
M23 = 8388607 |
|||
= 47 × 178481 |
|||
M29 = 536870911 |
|||
= 233 × 2304167 |
|||
M31 = 2147483647 is prime |
|||
M37 = 137438953471 |
|||
= 223 × 616318177 |
|||
M41 = 2199023255551 |
|||
= 13367 × 164511353 |
|||
M43 = 8796093022207 |
|||
= 431 × 20408568497 |
|||
M47 = 140737488355327 |
|||
= 2351 × 59862819377 |
|||
M53 = 9007199254740991 |
|||
= 6361 × 1416003655831 |
|||
M59 = 576460752303423487 |
|||
= 179951 × 3203431780337 |
|||
M929 = 4538015467766671944574165338592225830478699345884382504442663144885072806275648112625635725391102144133907238129251016389326737199538896813326509341743147661691195191795226666084858428449394948944821764472508048114220424520501343042471615418544488778723282182172070046459244838911 |
|||
= 13007 × 348890248924938259750454781163390930305120269538278042934009621348894657205785201247454118966026150852149399410259938217062100192168747352450719561908445272675574320888385228421992652298715687625495638077382028762529439880103124705348782610789919949159935587158612289264184273</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
REXX practically has no limit (well, up to around 8 million) on the number of decimal digits (precision). |
REXX practically has no limit (well, up to around 8 million) on the number of decimal digits (precision). |
||
This REXX version automatically adjusts the '''numeric digits''' to whatever is needed. |
This REXX version automatically adjusts the '''numeric digits''' to whatever is needed. |
||
< |
<syntaxhighlight lang="rexx">/*REXX program uses exponent─and─mod operator to test possible Mersenne numbers. */ |
||
numeric digits 20 /*this will be increased if necessary. */ |
numeric digits 20 /*this will be increased if necessary. */ |
||
parse arg N spec /*obtain optional arguments from the CL*/ |
parse arg N spec /*obtain optional arguments from the CL*/ |
||
if N=='' | N=="," |
if N=='' | N=="," then N= 88 /*Not specified? Then use the default.*/ |
||
if spec=='' | spec=="," |
if spec=='' | spec=="," then spec= 920 970 /* " " " " " " */ |
||
do j=1; |
do j=1; z= j /*process a range, & then do some more.*/ |
||
if j=N |
if j==N then j= word(spec, 1) /*now, use the high range of numbers. */ |
||
if j>word(spec, 2) then leave /*done with |
if j>word(spec, 2) then leave /*done with " " " " " */ |
||
if \isPrime(z) then iterate /*if Z isn't a prime, keep plugging.*/ |
if \isPrime(z) then iterate /*if Z isn't a prime, keep plugging.*/ |
||
r=testMer(z) |
r= commas( testMer(z) ); L= length(r) /*add commas; get its new length. */ |
||
if r==0 then say right('M'z, 10) "──────── is a Mersenne prime." |
|||
r=commas(r); L=length(r) /*add commas to R; get it's new length.*/ |
|||
else say right('M'z, 50) "is composite, a factor:"right(r, max(L, 13) ) |
|||
else say right('M'z, 50) "is composite, a factor:" right(r, max(L, 11) ) |
|||
end /*j*/ |
end /*j*/ |
||
exit /*stick a fork in it, we're all done. */ |
exit /*stick a fork in it, we're all done. */ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
commas: |
commas: parse arg _; do jc=length(_)-3 to 1 by -3; _=insert(',', _, jc); end; return _ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
isPrime: procedure; parse arg x; if wordpos(x, '2 3 5 7') \== 0 then return 1 |
isPrime: procedure; parse arg x; if wordpos(x, '2 3 5 7') \== 0 then return 1 |
||
Line 2,311: | Line 2,674: | ||
end /*j*/ /* └─◄ Is j>√ x ? Then return 1*/ |
end /*j*/ /* └─◄ Is j>√ x ? Then return 1*/ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
iSqrt: procedure; parse arg x; |
iSqrt: procedure; parse arg x; #= 1; r= 0; do while #<=x; #= # * 4 |
||
end /*while*/ |
|||
do while #>1; #= # % 4; _= x-r-#; r= r % 2 |
|||
if _>=0 then do; x= _; r= r + # |
|||
end |
|||
end /*while*/ /*iSqrt ≡ integer square root.*/ |
end /*while*/ /*iSqrt ≡ integer square root.*/ |
||
return r /*───── ─ ── ─ ─ */ |
return r /*───── ─ ── ─ ─ */ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
testMer: procedure; parse arg x; p=2**x |
testMer: procedure; parse arg x; p= 2**x /* [↓] do we have enough digits?*/ |
||
$$=x2b( d2x(x) ) + 0 |
$$=x2b( d2x(x) ) + 0 |
||
if pos('E',p)\==0 then do; parse var p "E" _; numeric digits _+2; p=2**x |
if pos('E',p)\==0 then do; parse var p "E" _; numeric digits _ + 2; p= 2**x |
||
end |
|||
!.= 1; !.1= 0; !.7= 0 /*array used for a quicker test. */ |
|||
R= iSqrt(p) /*obtain integer square root of P*/ |
|||
do k=2 by 2; q= k*x + 1 /*(shortcut) compute value of Q. */ |
|||
m= q // 8 /*obtain the remainder when ÷ 8. */ |
|||
if !.m then iterate /*M must be either one or seven.*/ |
|||
parse var q '' -1 _; if _==5 then iterate /*last digit a five ? */ |
|||
if q// |
if q// 3==0 then iterate /*divisible by three? */ |
||
if q// 7==0 then iterate /* " " seven? */ |
|||
if q//11==0 then iterate /* " " eleven?*/ |
|||
/* ____ */ |
|||
if q>R then return 0 /*Is q>√2**x ? A Mersenne prime*/ |
if q>R then return 0 /*Is q>√2**x ? A Mersenne prime*/ |
||
sq=1; |
sq= 1; $= $$ /*obtain binary version from $. */ |
||
do until $==''; sq=sq*sq |
do until $==''; sq= sq*sq |
||
parse var $ _ 2 $ /*obtain 1st digit and the rest. */ |
parse var $ _ 2 $ /*obtain 1st digit and the rest. */ |
||
if _ then sq=(sq+sq) // q |
if _ then sq= (sq+sq) // q |
||
end /*until*/ |
end /*until*/ |
||
if sq==1 then return q /*Not a prime? Return a factor.*/ |
if sq==1 then return q /*Not a prime? Return a factor.*/ |
||
end /*k*/</ |
end /*k*/</syntaxhighlight> |
||
Program note: the '''iSqrt''' function computes the integer square root of a non-negative integer without using any floating point, just integers. |
Program note: the '''iSqrt''' function computes the integer square root of a non-negative integer without using any floating point, just integers. |
||
Line 2,343: | Line 2,712: | ||
M5 ──────── is a Mersenne prime. |
M5 ──────── is a Mersenne prime. |
||
M7 ──────── is a Mersenne prime. |
M7 ──────── is a Mersenne prime. |
||
M11 is composite, a factor: 23 |
M11 is composite, a factor: 23 |
||
M13 ──────── is a Mersenne prime. |
M13 ──────── is a Mersenne prime. |
||
M17 ──────── is a Mersenne prime. |
M17 ──────── is a Mersenne prime. |
||
M19 ──────── is a Mersenne prime. |
M19 ──────── is a Mersenne prime. |
||
M23 is composite, a factor: 47 |
M23 is composite, a factor: 47 |
||
M29 is composite, a factor: 233 |
M29 is composite, a factor: 233 |
||
M31 ──────── is a Mersenne prime. |
M31 ──────── is a Mersenne prime. |
||
M37 is composite, a factor: 223 |
M37 is composite, a factor: 223 |
||
M41 is composite, a factor: 13,367 |
M41 is composite, a factor: 13,367 |
||
M43 is composite, a factor: 431 |
M43 is composite, a factor: 431 |
||
M47 is composite, a factor: 2,351 |
M47 is composite, a factor: 2,351 |
||
M53 is composite, a factor: 6,361 |
M53 is composite, a factor: 6,361 |
||
M59 is composite, a factor: 179,951 |
M59 is composite, a factor: 179,951 |
||
M61 ──────── is a Mersenne prime. |
M61 ──────── is a Mersenne prime. |
||
M67 is composite, a factor: 193,707,721 |
M67 is composite, a factor: 193,707,721 |
||
M71 is composite, a factor: 228,479 |
M71 is composite, a factor: 228,479 |
||
M73 is composite, a factor: 439 |
M73 is composite, a factor: 439 |
||
M79 is composite, a factor: 2,687 |
M79 is composite, a factor: 2,687 |
||
M83 is composite, a factor: 167 |
M83 is composite, a factor: 167 |
||
M929 is composite, a factor: 13,007 |
M929 is composite, a factor: 13,007 |
||
M937 is composite, a factor: 28,111 |
M937 is composite, a factor: 28,111 |
||
M941 is composite, a factor: 7,529 |
M941 is composite, a factor: 7,529 |
||
M947 is composite, a factor: 295,130,657 |
M947 is composite, a factor: 295,130,657 |
||
M953 is composite, a factor: 343,081 |
M953 is composite, a factor: 343,081 |
||
M967 is composite, a factor: 23,209 |
M967 is composite, a factor: 23,209 |
||
</pre> |
</pre> |
||
=={{header|Ring}}== |
=={{header|Ring}}== |
||
< |
<syntaxhighlight lang="ring"> |
||
# Project : Factors of a Mersenne number |
# Project : Factors of a Mersenne number |
||
Line 2,413: | Line 2,782: | ||
end |
end |
||
return y |
return y |
||
</syntaxhighlight> |
|||
</lang> |
|||
Output: |
Output: |
||
<pre> |
<pre> |
||
Line 2,419: | Line 2,788: | ||
A factor of M937 is 28111 |
A factor of M937 is 28111 |
||
</pre> |
</pre> |
||
=={{header|RPL}}== |
|||
{{works with|HP|48}} |
|||
<code>PRIM?</code> is defined at [[Primality by trial division#RPL|Primality by trial division]] |
|||
{| class="wikitable" ≪ |
|||
! RPL code |
|||
! Comment |
|||
|- |
|||
| |
|||
≪ SWAP R→B → quotient power |
|||
≪ 2 power B→R LN 2 LN / FLOOR ^ R→B |
|||
1 |
|||
'''WHILE''' OVER B→R '''REPEAT''' |
|||
SQ |
|||
'''IF''' OVER power AND B→R '''THEN''' DUP + '''END''' |
|||
quotient MOD |
|||
SWAP SR SWAP |
|||
'''END''' SWAP DROP |
|||
≫ ≫ '<span style="color:blue">MODPOW</span>' STO |
|||
≪ 2 OVER ^ 1 - √ 0 → power max k |
|||
≪ 1 |
|||
'''WHILE''' 'k' INCR 2 * 1 + DUP max ≤ '''REPEAT''' |
|||
'''IF''' { 1 7 } OVER 8 MOD POS THEN |
|||
'''IF''' DUP <span style="color:blue">PRIM?</span> THEN |
|||
'''IF''' power OVER <span style="color:blue">MODPOW</span> 1 == '''THEN''' |
|||
SWAP max 'k' STO '''END''' |
|||
'''END END''' |
|||
DROP |
|||
'''END''' DROP |
|||
≫ '<span style="color:blue">MFACT</span>' STO |
|||
| |
|||
<span style="color:blue">MODPOW</span> ''( power quotient → remainder )'' |
|||
create top-bit mask |
|||
square = 1 |
|||
while mask is not zero |
|||
square *= square |
|||
if unmasked bit = 1 then square += square |
|||
square = square mod quotient |
|||
shift mask right |
|||
clean stack |
|||
return square |
|||
<span style="color:blue">MFACT</span> ''( N → factor ) '' |
|||
factor = 1 |
|||
while 2k+1 ≤ sqrt(M(N)) |
|||
if 2k+1 mod 8 is equal to 1 or 7 |
|||
if 2k+1 prime |
|||
is 2K+1 a factor of M(N) ? |
|||
if yes, exit loop |
|||
return factor |
|||
|} |
|||
929 <span style="color:blue">MFACT</span> |
|||
{{out}} |
|||
<pre> |
|||
1: 13007 |
|||
</pre> |
|||
Factor found in 69 minutes on a 4-bit HP-48SX calculator. |
|||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
{{works with|Ruby|1.9.3+}} |
{{works with|Ruby|1.9.3+}} |
||
< |
<syntaxhighlight lang="ruby">require 'prime' |
||
def mersenne_factor(p) |
def mersenne_factor(p) |
||
Line 2,455: | Line 2,885: | ||
Prime.each(53) { |p| check_mersenne p } |
Prime.each(53) { |p| check_mersenne p } |
||
check_mersenne 929</ |
check_mersenne 929</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,475: | Line 2,905: | ||
M53 = 2**53-1 is composite with factor 6361 |
M53 = 2**53-1 is composite with factor 6361 |
||
M929 = 2**929-1 is composite with factor 13007</pre> |
M929 = 2**929-1 is composite with factor 13007</pre> |
||
=={{header|Rust}}== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="rust">fn bit_count(mut n: usize) -> usize { |
|||
let mut count = 0; |
|||
while n > 0 { |
|||
n >>= 1; |
|||
count += 1; |
|||
} |
|||
count |
|||
} |
|||
fn mod_pow(p: usize, n: usize) -> usize { |
|||
let mut square = 1; |
|||
let mut bits = bit_count(p); |
|||
while bits > 0 { |
|||
square = square * square; |
|||
bits -= 1; |
|||
if (p & (1 << bits)) != 0 { |
|||
square <<= 1; |
|||
} |
|||
square %= n; |
|||
} |
|||
return square; |
|||
} |
|||
fn is_prime(n: usize) -> bool { |
|||
if n < 2 { |
|||
return false; |
|||
} |
|||
if n % 2 == 0 { |
|||
return n == 2; |
|||
} |
|||
if n % 3 == 0 { |
|||
return n == 3; |
|||
} |
|||
let mut p = 5; |
|||
while p * p <= n { |
|||
if n % p == 0 { |
|||
return false; |
|||
} |
|||
p += 2; |
|||
if n % p == 0 { |
|||
return false; |
|||
} |
|||
p += 4; |
|||
} |
|||
true |
|||
} |
|||
fn find_mersenne_factor(p: usize) -> usize { |
|||
let mut k = 0; |
|||
loop { |
|||
k += 1; |
|||
let q = 2 * k * p + 1; |
|||
if q % 8 == 1 || q % 8 == 7 { |
|||
if mod_pow(p, q) == 1 && is_prime(p) { |
|||
return q; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
fn main() { |
|||
println!("{}", find_mersenne_factor(929)); |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
13007 |
|||
</pre> |
|||
=={{header|Scala}}== |
=={{header|Scala}}== |
||
Line 2,480: | Line 2,981: | ||
===Full-blown version=== |
===Full-blown version=== |
||
<syntaxhighlight lang="scala"> |
|||
<lang Scala>/** Find factors of a Mersenne number |
|||
/** Find factors of a Mersenne number |
|||
* |
* |
||
* |
* The implementation finds factors for M929 and further. |
||
* |
* |
||
* |
* @example M59 = 2^059 - 1 = 576460752303423487 ( 2 msec) |
||
* |
* @example = 179951 × 3203431780337. |
||
*/ |
*/ |
||
object |
object FactorsOfAMersenneNumber extends App { |
||
val two: BigInt = 2 |
val two: BigInt = 2 |
||
def sieve(nums: Stream[Int]): Stream[Int] = |
|||
Stream.cons(nums.head, sieve((nums.tail) filter (_ % nums.head != 0))) |
|||
// An infinite stream of primes, lazy evaluation and memo-ized |
// An infinite stream of primes, lazy evaluation and memo-ized |
||
val oddPrimes = sieve( |
val oddPrimes = sieve(LazyList.from(3, 2)) |
||
def primes = sieve(2 #:: oddPrimes) |
|||
def |
def sieve(nums: LazyList[Int]): LazyList[Int] = |
||
LazyList.cons(nums.head, sieve((nums.tail) filter (_ % nums.head != 0))) |
|||
def primes: LazyList[Int] = sieve(2 #:: oddPrimes) |
|||
def factorMersenne(p: Int): Option[Long] = { |
def factorMersenne(p: Int): Option[Long] = { |
||
Line 2,507: | Line 3,008: | ||
// Build a stream of factors from (2*p+1) step-by (2*p) |
// Build a stream of factors from (2*p+1) step-by (2*p) |
||
def s(a: Long): |
def s(a: Long): LazyList[Long] = a #:: s(a + (2 * p)) // Build stream of possible factors |
||
// Limit and Filter Stream and then take the head element |
// Limit and Filter Stream and then take the head element |
||
Line 2,513: | Line 3,014: | ||
e.headOption |
e.headOption |
||
} |
} |
||
def mersenne(p: Int): BigInt = (two pow p) - 1 |
|||
// Test |
// Test |
||
(primes takeWhile (_ <= 97)) ++ List(929, 937) foreach { p => |
(primes takeWhile (_ <= 97)) ++ List(929, 937) foreach { p => { // Needs some intermediate results for nice formatting |
||
val nMersenne = mersenne(p); |
|||
{ // Needs some intermediate results for nice formatting |
|||
val lit = s"${nMersenne}" |
|||
val preAmble = f"${s"M${p}"}%4s = 2^$p%03d - 1 = ${lit}%s" |
|||
val datum = System.nanoTime |
|||
val result = factorMersenne(p) |
|||
val mSec = ((System.nanoTime - datum) / 1.0e+6).round |
|||
def decStr = { |
|||
if (lit.length > 30) f"(M has ${lit.length}%3d dec)" else "" |
|||
} |
|||
def sPrime: String = { |
|||
println(f"$preAmble${sPrime} ${f"($mSec%,1d"}%13s msec)") |
|||
if ( |
if (result.isEmpty) " is a Mersenne prime number." else " " * 28 |
||
println(f"${decStr}%-17s = ${result.get} × ${nMersenne / result.get}") |
|||
} |
} |
||
println(f"$preAmble${sPrime} ${f"($mSec%,1d"}%13s msec)") |
|||
if (result.isDefined) |
|||
println(f"${decStr}%-17s = ${result.get} × ${nMersenne / result.get}") |
|||
} |
} |
||
} |
|||
}</lang> |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre style="height:40ex;overflow:scroll"> M2 = 2^002 - 1 = 3 is a Mersenne prime number. (63 msec) |
<pre style="height:40ex;overflow:scroll"> M2 = 2^002 - 1 = 3 is a Mersenne prime number. (63 msec) |
||
Line 2,584: | Line 3,093: | ||
This works with PLT Scheme, other implementations only need to change the inclusion. |
This works with PLT Scheme, other implementations only need to change the inclusion. |
||
< |
<syntaxhighlight lang="scheme"> |
||
#lang scheme |
#lang scheme |
||
Line 2,606: | Line 3,115: | ||
(= 1 (modpow p i)))) |
(= 1 (modpow p i)))) |
||
i)) |
i)) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,618: | Line 3,127: | ||
=={{header|Seed7}}== |
=={{header|Seed7}}== |
||
< |
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i"; |
||
const func boolean: isPrime (in integer: number) is func |
const func boolean: isPrime (in integer: number) is func |
||
Line 2,681: | Line 3,190: | ||
begin |
begin |
||
writeln("Factor of M929: " <& mersenneFactor(929)); |
writeln("Factor of M929: " <& mersenneFactor(929)); |
||
end func;</ |
end func;</syntaxhighlight> |
||
Original source: [http://seed7.sourceforge.net/algorith/math.htm#isPrime isPrime], |
Original source: [http://seed7.sourceforge.net/algorith/math.htm#isPrime isPrime], |
||
Line 2,692: | Line 3,201: | ||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
< |
<syntaxhighlight lang="ruby">func mtest(b, p) { |
||
var bits = b.base(2).digits |
var bits = b.base(2).digits |
||
for (var sq = 1; bits; sq %= p) { |
for (var sq = 1; bits; sq %= p) { |
||
Line 2,712: | Line 3,221: | ||
say (f ? "M#{m} is composite with factor #{q}" |
say (f ? "M#{m} is composite with factor #{q}" |
||
: "M#{m} is prime") |
: "M#{m} is prime") |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,734: | Line 3,243: | ||
M929 is composite with factor 13007 |
M929 is composite with factor 13007 |
||
</pre> |
</pre> |
||
=={{header|Swift}}== |
|||
<syntaxhighlight lang="swift">import Foundation |
|||
extension BinaryInteger { |
|||
var isPrime: Bool { |
|||
if self == 0 || self == 1 { |
|||
return false |
|||
} else if self == 2 { |
|||
return true |
|||
} |
|||
let max = Self(ceil((Double(self).squareRoot()))) |
|||
for i in stride(from: 2, through: max, by: 1) where self % i == 0 { |
|||
return false |
|||
} |
|||
return true |
|||
} |
|||
func modPow(exp: Self, mod: Self) -> Self { |
|||
guard exp != 0 else { |
|||
return 1 |
|||
} |
|||
var res = Self(1) |
|||
var base = self % mod |
|||
var exp = exp |
|||
while true { |
|||
if exp & 1 == 1 { |
|||
res *= base |
|||
res %= mod |
|||
} |
|||
if exp == 1 { |
|||
return res |
|||
} |
|||
exp >>= 1 |
|||
base *= base |
|||
base %= mod |
|||
} |
|||
} |
|||
} |
|||
func mFactor(exp: Int) -> Int? { |
|||
for k in 0..<16384 { |
|||
let q = 2*exp*k + 1 |
|||
if !q.isPrime { |
|||
continue |
|||
} else if q % 8 != 1 && q % 8 != 7 { |
|||
continue |
|||
} else if 2.modPow(exp: exp, mod: q) == 1 { |
|||
return q |
|||
} |
|||
} |
|||
return nil |
|||
} |
|||
print(mFactor(exp: 929)!) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>13007</pre> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
For <code>primes::is_prime</code> see [[Prime decomposition#Tcl]] |
For <code>primes::is_prime</code> see [[Prime decomposition#Tcl]] |
||
< |
<syntaxhighlight lang="tcl">proc int2bits {n} { |
||
binary scan [binary format I1 $n] B* binstring |
binary scan [binary format I1 $n] B* binstring |
||
return [split [string trimleft $binstring 0] ""] |
return [split [string trimleft $binstring 0] ""] |
||
Line 2,781: | Line 3,360: | ||
} else { |
} else { |
||
puts "no factor found for M$exp" |
puts "no factor found for M$exp" |
||
}</ |
}</syntaxhighlight> |
||
=={{header|TI-83 BASIC}}== |
=={{header|TI-83 BASIC}}== |
||
Line 2,787: | Line 3,366: | ||
{{works with|TI-83 BASIC|TI-84Plus 2.55MP}} |
{{works with|TI-83 BASIC|TI-84Plus 2.55MP}} |
||
The program uses the new remainder function from OS 2.53MP, if not available it can be replaced by: |
The program uses the new remainder function from OS 2.53MP, if not available it can be replaced by: |
||
< |
<syntaxhighlight lang="ti83b">remainder(A,B) equivalent to iPart(B*fPart(A/B))</syntaxhighlight>Due to several problems, no Goto has been used. As a matter of fact the version is clearer. |
||
< |
<syntaxhighlight lang="ti83b">Prompt Q |
||
1→K:0→T |
1→K:0→T |
||
While K≤2^20 and T=0 |
While K≤2^20 and T=0 |
||
Line 2,820: | Line 3,399: | ||
End |
End |
||
If T=0:0→F |
If T=0:0→F |
||
Disp Q,F</ |
Disp Q,F</syntaxhighlight> |
||
{{in}} |
{{in}} |
||
<pre> |
<pre> |
||
Line 2,833: | Line 3,412: | ||
=={{header|uBasic/4tH}}== |
=={{header|uBasic/4tH}}== |
||
<lang>Print "A factor of M929 is "; FUNC(_FNmersenne_factor(929)) |
<syntaxhighlight lang="text">Print "A factor of M929 is "; FUNC(_FNmersenne_factor(929)) |
||
Print "A factor of M937 is "; FUNC(_FNmersenne_factor(937)) |
Print "A factor of M937 is "; FUNC(_FNmersenne_factor(937)) |
||
Line 2,886: | Line 3,465: | ||
Next |
Next |
||
Return (d@)</ |
Return (d@)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>A factor of M929 is 13007 |
<pre>A factor of M929 is 13007 |
||
Line 2,895: | Line 3,474: | ||
=={{header|VBScript}}== |
=={{header|VBScript}}== |
||
{{trans|REXX}} |
{{trans|REXX}} |
||
< |
<syntaxhighlight lang="vb">' Factors of a Mersenne number |
||
for i=1 to 59 |
for i=1 to 59 |
||
z=i |
z=i |
||
Line 2,957: | Line 3,536: | ||
loop |
loop |
||
testM=0 |
testM=0 |
||
end function</ |
end function</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,982: | Line 3,561: | ||
{{trans|BBC BASIC}} |
{{trans|BBC BASIC}} |
||
{{works with|Visual Basic|VB6 Standard}} |
{{works with|Visual Basic|VB6 Standard}} |
||
< |
<syntaxhighlight lang="vb">Sub mersenne() |
||
Dim q As Long, k As Long, p As Long, d As Long |
Dim q As Long, k As Long, p As Long, d As Long |
||
Dim factor As Long, i As Long, y As Long, z As Long |
Dim factor As Long, i As Long, y As Long, z As Long |
||
Line 3,017: | Line 3,596: | ||
okfactor: |
okfactor: |
||
Debug.Print "M" & q, "factor=" & factor |
Debug.Print "M" & q, "factor=" & factor |
||
End Sub</ |
End Sub</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
M47 factor=2351 |
M47 factor=2351 |
||
</pre> |
|||
=={{header|V (Vlang)}}== |
|||
{{trans|go}} |
|||
<syntaxhighlight lang="go">import math |
|||
const qlimit = int(2e8) |
|||
fn main() { |
|||
mtest(31) |
|||
mtest(67) |
|||
mtest(929) |
|||
} |
|||
fn mtest(m int) { |
|||
// the function finds odd prime factors by |
|||
// searching no farther than sqrt(N), where N = 2^m-1. |
|||
// the first odd prime is 3, 3^2 = 9, so M3 = 7 is still too small. |
|||
// M4 = 15 is first number for which test is meaningful. |
|||
if m < 4 { |
|||
println("$m < 4. M$m not tested.") |
|||
return |
|||
} |
|||
flimit := math.sqrt(math.pow(2, f64(m)) - 1) |
|||
mut qlast := 0 |
|||
if flimit < qlimit { |
|||
qlast = int(flimit) |
|||
} else { |
|||
qlast = qlimit |
|||
} |
|||
mut composite := []bool{len: qlast+1} |
|||
sq := int(math.sqrt(f64(qlast))) |
|||
loop: |
|||
for q := int(3); ; { |
|||
if q <= sq { |
|||
for i := q * q; i <= qlast; i += q { |
|||
composite[i] = true |
|||
} |
|||
} |
|||
q8 := q % 8 |
|||
if (q8 == 1 || q8 == 7) && mod_pow(2, m, q) == 1 { |
|||
println("M$m has factor $q") |
|||
return |
|||
} |
|||
for { |
|||
q += 2 |
|||
if q > qlast { |
|||
break loop |
|||
} |
|||
if !composite[q] { |
|||
break |
|||
} |
|||
} |
|||
} |
|||
println("No factors of M$m found.") |
|||
} |
|||
// base b to power p, mod m |
|||
fn mod_pow(b int, p int, m int) int { |
|||
mut pow := i64(1) |
|||
b64 := i64(b) |
|||
m64 := i64(m) |
|||
mut bit := u32(30) |
|||
for 1<<bit&p == 0 { |
|||
bit-- |
|||
} |
|||
for { |
|||
pow *= pow |
|||
if 1<<bit&p != 0 { |
|||
pow *= b64 |
|||
} |
|||
pow %= m64 |
|||
if bit == 0 { |
|||
break |
|||
} |
|||
bit-- |
|||
} |
|||
return int(pow) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
No factors of M31 found. |
|||
M67 has factor 193707721 |
|||
M929 has factor 13007 |
|||
</pre> |
|||
=={{header|Wren}}== |
|||
{{trans|JavaScript}} |
|||
{{libheader|Wren-math}} |
|||
{{libheader|Wren-fmt}} |
|||
<syntaxhighlight lang="wren">import "./math" for Int |
|||
import "./fmt" for Conv, Fmt |
|||
var trialFactor = Fn.new { |base, exp, mod| |
|||
var square = 1 |
|||
var bits = Conv.itoa(exp, 2).toList |
|||
var ln = bits.count |
|||
for (i in 0...ln) { |
|||
square = square * square * (bits[i] == "1" ? base : 1) % mod |
|||
} |
|||
return square == 1 |
|||
} |
|||
var mersenneFactor = Fn.new { |p| |
|||
var limit = (2.pow(p) - 1).sqrt.floor |
|||
var k = 1 |
|||
while ((2*k*p - 1) < limit) { |
|||
var q = 2*k*p + 1 |
|||
if (Int.isPrime(q) && (q%8 == 1 || q%8 == 7) && trialFactor.call(2, p, q)) { |
|||
return q // q is a factor of 2^p - 1 |
|||
} |
|||
k = k + 1 |
|||
} |
|||
return null |
|||
} |
|||
var m = [3, 5, 11, 17, 23, 29, 31, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 97, 929] |
|||
for (p in m) { |
|||
var f = mersenneFactor.call(p) |
|||
Fmt.write("2^$3d - 1 is ", p) |
|||
if (f) { |
|||
Fmt.print("composite (factor $d)", f) |
|||
} else { |
|||
System.print("prime") |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
2^ 3 - 1 is prime |
|||
2^ 5 - 1 is prime |
|||
2^ 11 - 1 is composite (factor 23) |
|||
2^ 17 - 1 is prime |
|||
2^ 23 - 1 is composite (factor 47) |
|||
2^ 29 - 1 is composite (factor 233) |
|||
2^ 31 - 1 is prime |
|||
2^ 37 - 1 is composite (factor 223) |
|||
2^ 41 - 1 is composite (factor 13367) |
|||
2^ 43 - 1 is composite (factor 431) |
|||
2^ 47 - 1 is composite (factor 2351) |
|||
2^ 53 - 1 is composite (factor 6361) |
|||
2^ 59 - 1 is composite (factor 179951) |
|||
2^ 67 - 1 is composite (factor 193707721) |
|||
2^ 71 - 1 is composite (factor 228479) |
|||
2^ 73 - 1 is composite (factor 439) |
|||
2^ 79 - 1 is composite (factor 2687) |
|||
2^ 83 - 1 is composite (factor 167) |
|||
2^ 97 - 1 is composite (factor 11447) |
|||
2^929 - 1 is composite (factor 13007) |
|||
</pre> |
</pre> |
||
=={{header|zkl}}== |
=={{header|zkl}}== |
||
{{trans|EchoLisp}} |
{{trans|EchoLisp}} |
||
< |
<syntaxhighlight lang="zkl">var [const] BN=Import("zklBigNum"); // libGMP |
||
// M = 2^P - 1 , P prime |
// M = 2^P - 1 , P prime |
||
Line 3,041: | Line 3,768: | ||
} |
} |
||
False |
False |
||
}</ |
}</syntaxhighlight> |
||
< |
<syntaxhighlight lang="zkl">m_divisor(929).println(); // 13007 |
||
m_divisor(4423).println(); // False |
m_divisor(4423).println(); // False |
||
(BN(2).pow(4423) - 1).probablyPrime().println(); // True</ |
(BN(2).pow(4423) - 1).probablyPrime().println(); // True</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Revision as of 16:09, 14 April 2024
You are encouraged to solve this task according to the task description, using any language you may know.
A Mersenne number is a number in the form of 2P-1.
If P is prime, the Mersenne number may be a Mersenne prime (if P is not prime, the Mersenne number is also not prime).
In the search for Mersenne prime numbers it is advantageous to eliminate exponents by finding a small factor before starting a, potentially lengthy, Lucas-Lehmer test.
There are very efficient algorithms for determining if a number divides 2P-1 (or equivalently, if 2P mod (the number) = 1). Some languages already have built-in implementations of this exponent-and-mod operation (called modPow or similar).
The following is how to implement this modPow yourself:
For example, let's compute 223 mod 47. Convert the exponent 23 to binary, you get 10111. Starting with square = 1, repeatedly square it. Remove the top bit of the exponent, and if it's 1 multiply square by the base of the exponentiation (2), then compute square modulo 47. Use the result of the modulo from the last step as the initial value of square in the next step:
remove optional square top bit multiply by 2 mod 47 ──────────── ─────── ───────────── ────── 1*1 = 1 1 0111 1*2 = 2 2 2*2 = 4 0 111 no 4 4*4 = 16 1 11 16*2 = 32 32 32*32 = 1024 1 1 1024*2 = 2048 27 27*27 = 729 1 729*2 = 1458 1
Since 223 mod 47 = 1, 47 is a factor of 2P-1. (To see this, subtract 1 from both sides: 223-1 = 0 mod 47.) Since we've shown that 47 is a factor, 223-1 is not prime. Further properties of Mersenne numbers allow us to refine the process even more. Any factor q of 2P-1 must be of the form 2kP+1, k being a positive integer or zero. Furthermore, q must be 1 or 7 mod 8. Finally any potential factor q must be prime. As in other trial division algorithms, the algorithm stops when 2kP+1 > sqrt(N).
These primality tests only work on Mersenne numbers where P is prime. For example, M4=15 yields no factors using these techniques, but factors into 3 and 5, neither of which fit 2kP+1.
- Task
Using the above method find a factor of 2929-1 (aka M929)
- Related tasks
- count in factors
- prime decomposition
- factors of an integer
- Sieve of Eratosthenes
- primality by trial division
- trial factoring of a Mersenne number
- partition an integer X into N primes
- sequence of primes by Trial Division
- See also
- Computers in 1948: 2127 - 1
(Note: This video is no longer available because the YouTube account associated with this video has been terminated.)
11l
F is_prime(a)
I a == 2 {R 1B}
I a < 2 | a % 2 == 0 {R 0B}
L(i) (3 .. Int(sqrt(a))).step(2)
I a % i == 0
R 0B
R 1B
F m_factor(p)
V max_k = 16384 I/ p
L(k) 0 .< max_k
V q = 2 * p * k + 1
I !is_prime(q)
L.continue
E I q % 8 != 1 & q % 8 != 7
L.continue
E I pow(2, p, q) == 1
R q
R 0
V exponent = Int(input(‘Enter exponent of Mersenne number: ’))
I !is_prime(exponent)
print(‘Exponent is not prime: #.’.format(exponent))
E
V factor = m_factor(exponent)
I factor == 0
print(‘No factor found for M#.’.format(exponent))
E
print(‘M#. has a factor: #.’.format(exponent, factor))
- Output:
Enter exponent of Mersenne number: 929 M929 has a factor: 13007
8086 Assembly
P: equ 929 ; P for 2^P-1
cpu 8086
bits 16
org 100h
section .text
mov ax,P ; Is P prime?
call prime
mov dx,notprm
jc msg ; If not, say so and stop.
xor bp,bp ; Let BP hold k
test_k: inc bp ; k += 1
mov ax,P ; Calculate 2kP + 1
mul bp ; AX = kP
shl ax,1 ; AX = 2kP
inc ax ; AX = 2kP + 1
mov dx,ovfl ; If AX overflows (16 bits), say so and stop
jc msg
mov bx,ax ; What is 2^P mod (2kP+1)?
mov cx,P
call modpow
dec ax ; If it is 1, we're done
jnz test_k ; If not, increment K and try again
mov dx,factor ; If so, we found a factor
call msg
prbx: mov ax,10 ; The factor is still in BX
xchg bx,ax ; Put factor in AX and divisor (10) in BX
mov di,number ; Generate ASCII representation of number
digit: xor dx,dx
div bx ; Divide current number by 10,
add dl,'0' ; add '0' to remainder,
dec di ; move pointer back,
mov [di],dl ; store digit,
test ax,ax ; and if there are more digits,
jnz digit ; find the next digit.
mov dx,di ; Finally, print the number.
jmp msg
;;; Calculate 2^CX mod BX
;;; Output: AX
;;; Destroyed: CX, DX
modpow: shl cx,1 ; Shift CX left until top bit in high bit
jnc modpow ; Keep shifting while carry zero
rcr cx,1 ; Put the top bit back into CX
mov ax,1 ; Start with square = 1
.step: mul ax ; Square (result is 32-bit, goes in DX:AX)
shl cx,1 ; Grab a bit from CX
jnc .nodbl ; If zero, don't multiply by two
shl ax,1 ; Shift DX:AX left (mul by two)
rcl dx,1
.nodbl: div bx ; Divide DX:AX by BX (putting modulus in DX)
mov ax,dx ; Continue with modulus
jcxz .done ; When CX reaches 0, we're done
jmp .step ; Otherwise, do the next step
.done: ret
;;; Is AX prime?
;;; Output: carry clear if prime, set if not prime.
;;; Destroyed: AX, BX, CX, DX, SI, DI, BP
prime: mov cx,[prcnt] ; See if AX is already in the list of primes
mov di,primes
repne scasw ; If so, return
je modpow.done ; Reuse the RET just above here (carry clear)
mov bp,ax ; Move AX out of the way
mov bx,[di-2] ; Start generating new primes
.sieve: inc bx ; BX = last prime + 2
inc bx
cmp bp,bx ; If BX higher than number to test,
jb modpow.done ; stop, number is not prime. (carry set)
mov cx,[prcnt] ; CX = amount of current primes
mov si,primes ; SI = start of primes
.try: mov ax,bx ; BX divisible by current prime?
xor dx,dx
div word [si]
test dx,dx ; If so, BX is not prime.
jz .sieve
inc si
inc si
loop .try ; Otherwise, try next prime.
mov ax,bx ; If we get here, BX _is_ prime
stosw ; So add it to the list
inc word [prcnt] ; We have one more prime
cmp ax,bp ; Is it the prime we are looking for?
jne .sieve ; If not, try next prime
ret
;;; Print message in DX
msg: mov ah,9
int 21h
ret
section .data
db "*****" ; Placeholder for number
number: db "$"
notprm: db "P is not prime.$"
ovfl: db "Range of factor exceeded (max 16 bits)."
factor: db "Found factor: $"
prcnt: dw 2 ; Amount of primes currently in list
primes: dw 2, 3 ; List of primes to be extended
- Output:
Found factor: 13007
360 Assembly
Use of bitwise operations (TM (Test under Mask), SLA (Shift Left Arithmetic),SRA (Shift Right Arithmetic)).
* Factors of a Mersenne number 11/09/2015
MERSENNE CSECT
USING MERSENNE,R15
MVC Q,=F'929' q=929 (M929=2**929-1)
LA R6,1 k=1
LOOPK C R6,=F'1048576' do k=1 to 2**20
BNL ELOOPK
LR R5,R6 k
M R4,Q *q
SLA R5,1 *2 by shift left 1
LA R5,1(R5) +1
ST R5,P p=k*q*2+1
L R2,P p
N R2,=F'7' p&7
C R2,=F'1' if ((p&7)=1) p='*001'
BE OK
C R2,=F'7' or if ((p&7)=7) p='*111'
BNE NOTOK
OK MVI PRIME,X'00' then prime=false is prime?
LA R2,2 loop count=2
LA R1,2 j=2 and after j=3
J2J3 L R4,P p
SRDA R4,32 r4>>r5
DR R4,R1 p/j
LTR R4,R4 if p//j=0
BZ NOTPRIME then goto notprime
LA R1,1(R1) j=j+1
BCT R2,J2J3
LA R7,5 d=5
WHILED LR R5,R7 d
MR R4,R7 *d
C R5,P do while(d*d<=p)
BH EWHILED
LA R2,2 loop count=2
LA R1,2 j=2 and after j=4
J2J4 L R4,P p
SRDA R4,32 r4>>r5
DR R4,R7 /d
LTR R4,R4 if p//d=0
BZ NOTPRIME then goto notprime
AR R7,R1 d=d+j
LA R1,2(R1) j=j+2
BCT R2,J2J4
B WHILED
EWHILED MVI PRIME,X'01' prime=true so is prime
NOTPRIME L R8,Q i=q
MVC Y,=F'1' y=1
MVC Z,=F'2' z=2
WHILEI LTR R8,R8 do while(i^=0)
BZ EWHILEI
ST R8,PG i
TM PG+3,B'00000001' if first bit of i not 1
BZ NOTFIRST
L R5,Y y
M R4,Z *z
LA R4,0
D R4,P /p
ST R4,Y y=(y*z)//p
NOTFIRST L R5,Z z
M R4,Z *z
LA R4,0
D R4,P /p
ST R4,Z z=(z*z)//p
SRA R8,1 i=i/2 by shift right 1
B WHILEI
EWHILEI CLI PRIME,X'01' if prime
BNE NOTOK
CLC Y,=F'1' and if y=1
BNE NOTOK
MVC FACTOR,P then factor=p
B OKFACTOR
NOTOK LA R6,1(R6) k=k+1
B LOOPK
ELOOPK MVC FACTOR,=F'0' factor=0
OKFACTOR L R1,Q
XDECO R1,PG edit q
L R1,FACTOR
XDECO R1,PG+12 edit factor
XPRNT PG,24 print
XR R15,R15
BR R14
PRIME DS X flag for prime
Q DS F
P DS F
Y DS F
Z DS F
FACTOR DS F a factor of q
PG DS CL24 buffer
YREGS
END MERSENNE
- Output:
929 13007
Ada
mersenne.adb:
with Ada.Text_IO;
-- reuse Is_Prime from [[Primality by Trial Division]]
with Is_Prime;
procedure Mersenne is
function Is_Set (Number : Natural; Bit : Positive) return Boolean is
begin
return Number / 2 ** (Bit - 1) mod 2 = 1;
end Is_Set;
function Get_Max_Bit (Number : Natural) return Natural is
Test : Natural := 0;
begin
while 2 ** Test <= Number loop
Test := Test + 1;
end loop;
return Test;
end Get_Max_Bit;
function Modular_Power (Base, Exponent, Modulus : Positive) return Natural is
Maximum_Bit : constant Natural := Get_Max_Bit (Exponent);
Square : Natural := 1;
begin
for Bit in reverse 1 .. Maximum_Bit loop
Square := Square ** 2;
if Is_Set (Exponent, Bit) then
Square := Square * Base;
end if;
Square := Square mod Modulus;
end loop;
return Square;
end Modular_Power;
Not_A_Prime_Exponent : exception;
function Get_Factor (Exponent : Positive) return Natural is
Factor : Positive;
begin
if not Is_Prime (Exponent) then
raise Not_A_Prime_Exponent;
end if;
for K in 1 .. 16384 / Exponent loop
Factor := 2 * K * Exponent + 1;
if Factor mod 8 = 1 or else Factor mod 8 = 7 then
if Is_Prime (Factor) and then Modular_Power (2, Exponent, Factor) = 1 then
return Factor;
end if;
end if;
end loop;
return 0;
end Get_Factor;
To_Test : constant Positive := 929;
Factor : Natural;
begin
Ada.Text_IO.Put ("2 **" & Integer'Image (To_Test) & " - 1 ");
begin
Factor := Get_Factor (To_Test);
if Factor = 0 then
Ada.Text_IO.Put_Line ("is prime.");
else
Ada.Text_IO.Put_Line ("has factor" & Integer'Image (Factor));
end if;
exception
when Not_A_Prime_Exponent =>
Ada.Text_IO.Put_Line ("is not a Mersenne number");
end;
end Mersenne;
- Output:
2 ** 929 - 1 has factor 13007
ALGOL 68
MODE ISPRIMEINT = INT;
PR READ "prelude/is_prime.a68" PR;
MODE POWMODSTRUCT = INT;
PR READ "prelude/pow_mod.a68" PR;
PROC m factor = (INT p)INT:BEGIN
INT m factor;
INT max k, msb, n, q;
FOR i FROM bits width - 2 BY -1 TO 0 WHILE ( BIN p SHR i AND 2r1 ) = 2r0 DO
msb := i
OD;
max k := ENTIER sqrt(max int) OVER p; # limit for k to prevent overflow of max int #
FOR k FROM 1 TO max k DO
q := 2*p*k + 1;
IF NOT is prime(q) THEN
SKIP
ELIF q MOD 8 /= 1 AND q MOD 8 /= 7 THEN
SKIP
ELSE
n := pow mod(2,p,q);
IF n = 1 THEN
m factor := q;
return
FI
FI
OD;
m factor := 0;
return:
m factor
END;
BEGIN
INT exponent, factor;
print("Enter exponent of Mersenne number:");
read(exponent);
IF NOT is prime(exponent) THEN
print(("Exponent is not prime: ", exponent, new line))
ELSE
factor := m factor(exponent);
IF factor = 0 THEN
print(("No factor found for M", exponent, new line))
ELSE
print(("M", exponent, " has a factor: ", factor, new line))
FI
FI
END
Example:
Enter exponent of Mersenne number:929 M +929 has a factor: +13007
Arturo
mersenneFactors: function [q][
if not? prime? q -> print "number not prime!"
r: new q
while -> r > 0
-> shl 'r 1
d: new 1 + 2 * q
while [true][
i: new 1
p: new r
while [p <> 0][
i: new (i * i) % d
if p < 0 -> 'i * 2
if i > d -> 'i - d
shl 'p 1
]
if? i <> 1 -> 'd + 2 * q
else -> break
]
print ["2 ^" q "- 1 = 0 ( mod" d ")"]
]
mersenneFactors 929
- Output:
2 ^ 929 - 1 = 0 ( mod 13007 )
AutoHotkey
ahk discussion
MsgBox % MFact(27) ;-1: 27 is not prime
MsgBox % MFact(2) ; 0
MsgBox % MFact(3) ; 0
MsgBox % MFact(5) ; 0
MsgBox % MFact(7) ; 0
MsgBox % MFact(11) ; 23
MsgBox % MFact(13) ; 0
MsgBox % MFact(17) ; 0
MsgBox % MFact(19) ; 0
MsgBox % MFact(23) ; 47
MsgBox % MFact(29) ; 233
MsgBox % MFact(31) ; 0
MsgBox % MFact(37) ; 223
MsgBox % MFact(41) ; 13367
MsgBox % MFact(43) ; 431
MsgBox % MFact(47) ; 2351
MsgBox % MFact(53) ; 6361
MsgBox % MFact(929) ; 13007
MFact(p) { ; blank if 2**p-1 can be prime, otherwise a prime divisor < 2**32
If !IsPrime32(p)
Return -1 ; Error (p must be prime)
Loop % 2.0**(p<64 ? p/2-1 : 31)/p ; test prime divisors < 2**32, up to sqrt(2**p-1)
If (((q:=2*p*A_Index+1)&7 = 1 || q&7 = 7) && IsPrime32(q) && PowMod(2,p,q)=1)
Return q
Return 0
}
IsPrime32(n) { ; n < 2**32
If n in 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
Return 1
If (!(n&1)||!mod(n,3)||!mod(n,5)||!mod(n,7)||!mod(n,11)||!mod(n,13)||!mod(n,17)||!mod(n,19))
Return 0
n1 := d := n-1, s := 0
While !(d&1)
d>>=1, s++
Loop 3 {
x := PowMod( A_Index=1 ? 2 : A_Index=2 ? 7 : 61, d, n)
If (x=1 || x=n1)
Continue
Loop % s-1
If (1 = x:=PowMod(x,2,n))
Return 0
Else If (x = n1)
Break
IfLess x,%n1%, Return 0
}
Return 1
}
PowMod(x,n,m) { ; x**n mod m
y := 1, i := n, z := x
While i>0
y := i&1 ? mod(y*z,m) : y, z := mod(z*z,m), i >>= 1
Return y
}
BBC BASIC
PRINT "A factor of M929 is "; FNmersenne_factor(929)
PRINT "A factor of M937 is "; FNmersenne_factor(937)
END
DEF FNmersenne_factor(P%)
LOCAL K%, Q%
IF NOT FNisprime(P%) THEN = -1
FOR K% = 1 TO 1000000
Q% = 2*K%*P% + 1
IF (Q% AND 7) = 1 OR (Q% AND 7) = 7 THEN
IF FNisprime(Q%) IF FNmodpow(2, P%, Q%) = 1 THEN = Q%
ENDIF
NEXT K%
= 0
DEF FNisprime(N%)
LOCAL D%
IF N% MOD 2=0 THEN = (N% = 2)
IF N% MOD 3=0 THEN = (N% = 3)
D% = 5
WHILE D% * D% <= N%
IF N% MOD D% = 0 THEN = FALSE
D% += 2
IF N% MOD D% = 0 THEN = FALSE
D% += 4
ENDWHILE
= TRUE
DEF FNmodpow(X%, N%, M%)
LOCAL I%, Y%, Z%
I% = N% : Y% = 1 : Z% = X%
WHILE I%
IF I% AND 1 THEN Y% = (Y% * Z%) MOD M%
Z% = (Z% * Z%) MOD M%
I% = I% >>> 1
ENDWHILE
= Y%
- Output:
A factor of M929 is 13007 A factor of M937 is 28111
Bracmat
( ( modPow
= square P divisor highbit log 2pow
. !arg:(?P.?divisor)
& 1:?square
& 2\L!P:#%?log+?
& 2^!log:?2pow
& whl
' ( mod
$ ( ( div$(!P.!2pow):1&2
| 1
)
* !square^2
. !divisor
)
: ?square
& mod$(!P.!2pow):?P
& 1/2*!2pow:~/:?2pow
)
& !square
)
& ( isPrime
= incs nextincs primeCandidate nextPrimeCandidate quotient
. 1 1 2 2 (4 2 4 2 4 6 2 6:?incs)
: ?nextincs
& 1:?primeCandidate
& ( nextPrimeCandidate
= ( !nextincs:&!incs:?nextincs
|
)
& !nextincs:%?inc ?nextincs
& !inc+!primeCandidate:?primeCandidate
)
& whl
' ( (!nextPrimeCandidate:?divisor)^2:~>!arg
& !arg*!divisor^-1:?quotient:/
)
& !quotient:/
)
& ( Factors-of-a-Mersenne-Number
= P k candidate bignum
. !arg:?P
& 2^!P+-1:?bignum
& 0:?k
& whl
' ( 2*(1+!k:?k)*!P+1:?candidate
& !candidate^2:~>!bignum
& ( ~(mod$(!candidate.8):(1|7))
| ~(isPrime$!candidate)
| modPow$(!P.!candidate):?mp:~1
)
)
& !mp:1
& (!candidate.(2^!P+-1)*!candidate^-1)
)
& ( Factors-of-a-Mersenne-Number$929:(?divisorA.?divisorB)
& out
$ ( str
$ ("found some divisors of 2^" !P "-1 : " !divisorA " and " !divisorB)
)
| out$"no divisors found"
)
);
- Output:
found some divisors of 2^!P-1 : 13007 and 348890248924938259750454781163390930305120269538278042934009621348894657205785 201247454118966026150852149399410259938217062100192168747352450719561908445272675574320888385228421992652298715687625495 638077382028762529439880103124705348782610789919949159935587158612289264184273
C
int isPrime(int n){
if (n%2==0) return n==2;
if (n%3==0) return n==3;
int d=5;
while(d*d<=n){
if(n%d==0) return 0;
d+=2;
if(n%d==0) return 0;
d+=4;}
return 1;}
main() {int i,d,p,r,q=929;
if (!isPrime(q)) return 1;
r=q;
while(r>0) r<<=1;
d=2*q+1;
do { for(p=r, i= 1; p; p<<= 1){
i=((long long)i * i) % d;
if (p < 0) i *= 2;
if (i > d) i -= d;}
if (i != 1) d += 2*q;
else break;
} while(1);
printf("2^%d - 1 = 0 (mod %d)\n", q, d);}
C#
using System;
namespace prog
{
class MainClass
{
public static void Main (string[] args)
{
int q = 929;
if ( !isPrime(q) ) return;
int r = q;
while( r > 0 )
r <<= 1;
int d = 2 * q + 1;
do
{
int i = 1;
for( int p=r; p!=0; p<<=1 )
{
i = (i*i) % d;
if (p < 0) i *= 2;
if (i > d) i -= d;
}
if (i != 1) d += 2 * q; else break;
}
while(true);
Console.WriteLine("2^"+q+"-1 = 0 (mod "+d+")");
}
static bool isPrime(int n)
{
if ( n % 2 == 0 ) return n == 2;
if ( n % 3 == 0 ) return n == 3;
int d = 5;
while( d*d <= n )
{
if ( n % d == 0 ) return false;
d += 2;
if ( n % d == 0 ) return false;
d += 4;
}
return true;
}
}
}
C++
#include <iostream>
#include <cstdint>
typedef uint64_t integer;
integer bit_count(integer n) {
integer count = 0;
for (; n > 0; count++)
n >>= 1;
return count;
}
integer mod_pow(integer p, integer n) {
integer square = 1;
for (integer bits = bit_count(p); bits > 0; square %= n) {
square *= square;
if (p & (1 << --bits))
square <<= 1;
}
return square;
}
bool is_prime(integer n) {
if (n < 2)
return false;
if (n % 2 == 0)
return n == 2;
for (integer p = 3; p * p <= n; p += 2)
if (n % p == 0)
return false;
return true;
}
integer find_mersenne_factor(integer p) {
for (integer k = 0, q = 1;;) {
q = 2 * ++k * p + 1;
if ((q % 8 == 1 || q % 8 == 7) && mod_pow(p, q) == 1 && is_prime(q))
return q;
}
return 0;
}
int main() {
std::cout << find_mersenne_factor(929) << std::endl;
return 0;
}
- Output:
13007
Clojure
(ns mersennenumber
(:gen-class))
(defn m* [p q m]
" Computes (p*q) mod m "
(mod (*' p q) m))
(defn power
"modular exponentiation (i.e. b^e mod m"
[b e m]
(loop [b b, e e, x 1]
(if (zero? e)
x
(if (even? e) (recur (m* b b m) (quot e 2) x)
(recur (m* b b m) (quot e 2) (m* b x m))))))
(defn divides? [k n]
" checks if k divides n "
(= (rem n k) 0))
(defn is-prime? [n]
" checks if n is prime "
(cond
(< n 2) false ; 0, 1 not prime (i.e. primes are greater than one)
(= n 2) true ; 2 is prime
(= 0 (mod n 2)) false ; all other evens are not prime
:else ; check for divisors up to sqrt(n)
(empty? (filter #(divides? % n) (take-while #(<= (* % %) n) (range 2 n))))))
;; Max k to check
(def MAX-K 16384)
(defn trial-factor [p k]
" check if k satisfies 2*k*P + 1 divides 2^p - 1 "
(let [q (+ (* 2 p k) 1)
mq (mod q 8)]
(cond
(not (is-prime? q)) nil
(and (not= 1 mq)
(not= 7 mq)) nil
(= 1 (power 2 p q)) q
:else nil)))
(defn m-factor [p]
" searches for k-factor "
(some #(trial-factor p %) (range 16384)))
(defn -main [p]
(if-not (is-prime? p)
(format "M%d = 2^%d - 1 exponent is not prime" p p)
(if-let [factor (m-factor p)]
(format "M%d = 2^%d - 1 is composite with factor %d" p p factor)
(format "M%d = 2^%d - 1 is prime" p p))))
;; Tests different p values
(doseq [p [2,3,4,5,7,11,13,17,19,23,29,31,37,41,43,47,53,929]
:let [s (-main p)]]
(println s))
- Output:
M2 = 2^2 - 1 is prime M3 = 2^3 - 1 is composite with factor 7 M4 = 2^4 - 1 exponent is not prime M5 = 2^5 - 1 is composite with factor 31 M7 = 2^7 - 1 is composite with factor 127 M11 = 2^11 - 1 is composite with factor 23 M13 = 2^13 - 1 is composite with factor 8191 M17 = 2^17 - 1 is composite with factor 131071 M19 = 2^19 - 1 is composite with factor 524287 M23 = 2^23 - 1 is composite with factor 47 M29 = 2^29 - 1 is composite with factor 233 M31 = 2^31 - 1 is prime M37 = 2^37 - 1 is composite with factor 223 M41 = 2^41 - 1 is composite with factor 13367 M43 = 2^43 - 1 is composite with factor 431 M47 = 2^47 - 1 is composite with factor 2351 M53 = 2^53 - 1 is composite with factor 6361 M929 = 2^929 - 1 is composite with factor 13007
CoffeeScript
mersenneFactor = (p) ->
limit = Math.sqrt(Math.pow(2,p) - 1)
k = 1
while (2*k*p - 1) < limit
q = 2*k*p + 1
if isPrime(q) and (q % 8 == 1 or q % 8 == 7) and trialFactor(2,p,q)
return q
k++
return null
isPrime = (value) ->
for i in [2...value]
return false if value % i == 0
return true if value % i != 0
trialFactor = (base, exp, mod) ->
square = 1
bits = exp.toString(2).split('')
for bit in bits
square = Math.pow(square, 2) * (if +bit is 1 then base else 1) % mod
return square == 1
checkMersenne = (p) ->
factor = mersenneFactor(+p)
console.log "M#{p} = 2^#{p}-1 is #{if factor is null then "prime" else "composite with #{factor}"}"
checkMersenne(prime) for prime in ["2","3","4","5","7","11","13","17","19","23","29","31","37","41","43","47","53","929"]
M2 = 2^2-1 is prime M3 = 2^3-1 is prime M4 = 2^4-1 is prime M5 = 2^5-1 is prime M7 = 2^7-1 is prime M11 = 2^11-1 is composite with 23 M13 = 2^13-1 is prime M17 = 2^17-1 is prime M19 = 2^19-1 is prime M23 = 2^23-1 is composite with 47 M29 = 2^29-1 is composite with 233 M31 = 2^31-1 is prime M37 = 2^37-1 is composite with 223 M41 = 2^41-1 is composite with 13367 M43 = 2^43-1 is composite with 431 M47 = 2^47-1 is composite with 2351 M53 = 2^53-1 is composite with 6361 M929 = 2^929-1 is composite with 13007
Common Lisp
(defun mersenne-fac (p &aux (m (1- (expt 2 p))))
(loop for k from 1
for n = (1+ (* 2 k p))
until (zerop (mod m n))
finally (return n)))
(print (mersenne-fac 929))
- Output:
13007
Version 2
We can use a primality test from the Primality by Trial Division task.
(defun primep (n)
"Is N prime?"
(and (> n 1)
(or (= n 2) (oddp n))
(loop for i from 3 to (isqrt n) by 2
never (zerop (rem n i)))))
Specific to this task, we define modulo-power and mersenne-prime-p.
(defun modulo-power (base power modulus)
(loop with square = 1
for bit across (format nil "~b" power)
do (setf square (* square square))
when (char= bit #\1) do (setf square (* square base))
do (setf square (mod square modulus))
finally (return square)))
(defun mersenne-prime-p (power)
(do* ((N (1- (expt 2 power)))
(sqN (isqrt N))
(k 1 (1+ k))
(q (1+ (* 2 power k)) (1+ (* 2 power k)))
(m (mod q 8) (mod q 8)))
((> q sqN) (values t))
(when (and (or (= 1 m) (= 7 m))
(primep q)
(= 1 (modulo-power 2 power q)))
(return (values nil q)))))
We can run the same tests from the Ruby entry.
> (loop for p in '(2 3 4 5 7 11 13 17 19 23 29 31 37 41 43 47 53 929) do (multiple-value-bind (primep factor) (mersenne-prime-p p) (format t "~&M~w = 2**~:*~w-1 is ~:[composite with factor ~w~;prime~]." p primep factor))) M2 = 2**2-1 is prime. M3 = 2**3-1 is prime. M4 = 2**4-1 is prime. M5 = 2**5-1 is prime. M7 = 2**7-1 is prime. M11 = 2**11-1 is composite with factor 23. M13 = 2**13-1 is prime. M17 = 2**17-1 is prime. M19 = 2**19-1 is prime. M23 = 2**23-1 is composite with factor 47. M29 = 2**29-1 is composite with factor 233. M31 = 2**31-1 is prime. M37 = 2**37-1 is composite with factor 223. M41 = 2**41-1 is composite with factor 13367. M43 = 2**43-1 is composite with factor 431. M47 = 2**47-1 is composite with factor 2351. M53 = 2**53-1 is composite with factor 6361. M929 = 2**929-1 is composite with factor 13007.
Crystal
require "big"
def prime?(n) # P3 Prime Generator primality test
return n | 1 == 3 if n < 5 # n: 0,1,4|false, 2,3|true
return false if n.gcd(6) != 1 # for n a P3 prime candidate (pc)
pc1, pc2 = -1, 1 # use P3's prime candidates sequence
until (pc1 += 6) > Math.sqrt(n).to_i # pcs are only 1/3 of all integers
return false if n % pc1 == 0 || n % (pc2 += 6) == 0 # if n is composite
end
true
end
# Compute b**e mod m
def powmod(b, e, m)
r, b = 1.to_big_i, b.to_big_i
while e > 0
r = (r * b) % m if e.odd?
b = (b * b) % m
e >>= 1
end
r
end
def mersenne_factor(p)
mers_num = 2.to_big_i ** p - 1
kp2 = p2 = 2.to_big_i * p
while (kp2 - 1) ** 2 < mers_num
q = kp2 + 1 # return q if it's a factor
return q if [1, 7].includes?(q % 8) && prime?(q) && (powmod(2, p, q) == 1)
kp2 += p2
end
true # could also set to `0` value to check for
end
def check_mersenne(p)
print "M#{p} = 2**#{p}-1 is "
f = mersenne_factor(p)
(puts "prime"; return) if f.is_a?(Bool) # or f == 0
puts "composite with factor #{f}"
end
(2..53).each { |p| check_mersenne(p) if prime?(p) }
check_mersenne 929
- Output:
M2 = 2**2-1 is prime M3 = 2**3-1 is prime M5 = 2**5-1 is prime M7 = 2**7-1 is prime M11 = 2**11-1 is composite with factor 23 M13 = 2**13-1 is prime M17 = 2**17-1 is prime M19 = 2**19-1 is prime M23 = 2**23-1 is composite with factor 47 M29 = 2**29-1 is composite with factor 233 M31 = 2**31-1 is prime M37 = 2**37-1 is composite with factor 223 M41 = 2**41-1 is composite with factor 13367 M43 = 2**43-1 is composite with factor 431 M47 = 2**47-1 is composite with factor 2351 M53 = 2**53-1 is composite with factor 6361 M929 = 2**929-1 is composite with factor 13007
D
import std.stdio, std.math, std.traits;
ulong mersenneFactor(in ulong p) pure nothrow @nogc {
static bool isPrime(T)(in T n) pure nothrow @nogc {
if (n < 2 || n % 2 == 0)
return n == 2;
for (Unqual!T i = 3; i ^^ 2 <= n; i += 2)
if (n % i == 0)
return false;
return true;
}
static ulong modPow(in ulong cb, in ulong ce,in ulong m)
pure nothrow @nogc {
ulong b = cb;
ulong result = 1;
for (ulong e = ce; e > 0; e >>= 1) {
if ((e & 1) == 1)
result = (result * b) % m;
b = (b ^^ 2) % m;
}
return result;
}
immutable ulong limit = p <= 64 ? cast(ulong)(real(2.0) ^^ p - 1).sqrt : uint.max; // prevents silent overflows
for (ulong k = 1; (2 * p * k + 1) < limit; k++) {
immutable ulong q = 2 * p * k + 1;
if ((q % 8 == 1 || q % 8 == 7) && isPrime(q) &&
modPow(2, p, q) == 1)
return q;
}
return 1; // returns a sensible smallest factor
}
void main() {
writefln("Factor of M929: %d", 929.mersenneFactor);
}
- Output:
Factor of M929: 13007
Delphi
See Pascal.
EasyLang
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
func bit_count n .
while n > 0
n = bitshift n -1
cnt += 1
.
return cnt
.
func mod_pow p n .
square = 1
bits = bit_count p
while bits > 0
square *= square
bits -= 1
if bitand p bitshift 1 bits > 0
square = bitshift square 1
.
square = square mod n
.
return square
.
func mersenne_factor p .
while 1 = 1
k += 1
q = 2 * k * p + 1
if (q mod 8 = 1 or q mod 8 = 7) and mod_pow p q = 1 and isprim q = 1
return q
.
.
.
print mersenne_factor 929
- Output:
13007
EchoLisp
;; M = 2^P - 1 , P prime
;; look for a prime divisor q such as : q < √ M, q = 1 or 7 modulo 8, q = 1 + 2kP
;; q is divisor if (powmod 2 P q) = 1
;; m-divisor returns q or #f
(define ( m-divisor P )
;; must limit the search as √ M may be HUGE
(define maxprime (min 1_000_000_000 (sqrt (expt 2 P))))
(for ((q (in-range 1 maxprime (* 2 P))))
#:when (member (modulo q 8) '(1 7))
#:when (prime? q)
#:break (= 1 (powmod 2 P q)) => q
#f ))
(m-divisor 929)
→ 13007
(m-divisor 4423)
→ #f
(lib 'bigint)
(prime? (1- (expt 2 4423))) ;; 2^4423 -1 is a Mersenne prime
→ #t
Elixir
defmodule Mersenne do
def mersenne_factor(p) do
limit = :math.sqrt(:math.pow(2, p) - 1)
mersenne_loop(p, limit, 1)
end
defp mersenne_loop(p, limit, k) when (2*k*p - 1) > limit, do: nil
defp mersenne_loop(p, limit, k) do
q = 2*k*p + 1
if prime?(q) and rem(q,8) in [1,7] and trial_factor(2,p,q),
do: q, else: mersenne_loop(p, limit, k+1)
end
defp trial_factor(base, exp, mod) do
Integer.digits(exp, 2)
|> Enum.reduce(1, fn bit,square ->
(square * square * (if bit==1, do: base, else: 1)) |> rem(mod)
end) == 1
end
def check_mersenne(p) do
IO.write "M#{p} = 2**#{p}-1 is "
f = mersenne_factor(p)
IO.puts if f, do: "composite with factor #{f}", else: "prime"
end
def prime?(n), do: prime?(n, :math.sqrt(n), 2)
defp prime?(_, limit, i) when limit < i, do: true
defp prime?(n, limit, i) do
if rem(n,i) == 0, do: false, else: prime?(n, limit, i+1)
end
end
[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,929]
|> Enum.each(fn p -> Mersenne.check_mersenne(p) end)
- Output:
M2 = 2**2-1 is prime M3 = 2**3-1 is prime M5 = 2**5-1 is prime M7 = 2**7-1 is prime M11 = 2**11-1 is composite with factor 23 M13 = 2**13-1 is prime M17 = 2**17-1 is prime M19 = 2**19-1 is prime M23 = 2**23-1 is composite with factor 47 M29 = 2**29-1 is composite with factor 233 M31 = 2**31-1 is prime M37 = 2**37-1 is composite with factor 223 M41 = 2**41-1 is composite with factor 13367 M43 = 2**43-1 is composite with factor 431 M47 = 2**47-1 is composite with factor 2351 M53 = 2**53-1 is composite with factor 6361 M929 = 2**929-1 is composite with factor 13007
Erlang
The modpow function is not my original. This is a translation of python, more or less.
-module(mersene2).
-export([prime/1,modpow/3,mf/1]).
mf(P) -> merseneFactor(P,math:sqrt(math:pow(2,P)-1),2).
merseneFactor(P,Limit,Acc) when Acc >= Limit -> io:write("None found");
merseneFactor(P,Limit,Acc) ->
Q = 2 * P * Acc + 1,
Isprime = prime(Q),
Mod = modpow(2,P,Q),
if
Isprime == false ->
merseneFactor(P,Limit,Acc+1);
Q rem 8 =/= 1 andalso Q rem 8 =/= 7 ->
merseneFactor(P,Limit,Acc+1);
Mod == 1 ->
io:format("M~w is composite with Factor: ~w~n",[P,Q]);
true -> merseneFactor(P,Limit,Acc+1)
end.
modpow(B, E, M) -> modpow(B, E, M, 1).
modpow(_B, E, _M, R) when E =< 0 -> R;
modpow(B, E, M, R) ->
R1 = case E band 1 =:= 1 of
true -> (R * B) rem M;
false -> R
end,
modpow( (B*B) rem M, E bsr 1, M, R1).
prime(N) -> divisors(N, N-1).
divisors(N, 1) -> true;
divisors(N, C) ->
case N rem C =:= 0 of
true -> false;
false -> divisors(N, C-1)
end.
- Output:
30> [ mersene2:mf(X) || X <- [37,41,43,47,53,92,929]]. M37 is composite with Factor: 223 M41 is composite with Factor: 13367 M43 is composite with Factor: 431 M47 is composite with Factor: 2351 M53 is composite with Factor: 6361 M92 is composite with Factor: 1657 M929 is composite with Factor: 13007 [ok,ok,ok,ok,ok,ok,ok]
Factor
USING: combinators.short-circuit interpolate io kernel locals
math math.bits math.functions math.primes sequences ;
IN: rosetta-code.mersenne-factors
: mod-pow-step ( square bit m -- square' )
[ [ sq ] [ [ 2 * ] when ] bi* ] dip mod ;
:: mod-pow ( m q -- n )
1 :> s! m make-bits <reversed>
[ s swap q mod-pow-step s! ] each s ;
: halt-search? ( m q N -- ? )
dupd > [
{
[ nip 8 mod [ 1 ] [ 7 ] bi [ = ] 2bi@ or ]
[ mod-pow 1 = ] [ nip prime? ]
} 2&&
] dip or ;
:: find-mersenne-factor ( m -- factor/f )
1 :> k!
2 m * 1 + :> q! ! the tentative factor.
2 m ^ sqrt :> N ! upper bound on search.
[ m q N halt-search? ] [ k 1 + k! 2 k * m * 1 + q! ] until
q N > f q ? ;
: test-mersenne ( m -- )
dup find-mersenne-factor
[ [I M${1} is not prime: factor ${0} found.I] ]
[ [I No factor found for M${}.I] ] if* nl ;
929 test-mersenne
- Output:
M929 is not prime: factor 13007 found.
Forth
: prime? ( odd -- ? )
3
begin 2dup dup * >=
while 2dup mod 0=
if 2drop false exit
then 2 +
repeat 2drop true ;
: 2-exp-mod { e m -- 2^e mod m }
1
0 30 do
e 1 i lshift >= if
dup *
e 1 i lshift and if 2* then
m mod
then
-1 +loop ;
: factor-mersenne ( exponent -- factor )
16384 over / dup 2 < abort" Exponent too large!"
1 do
dup i * 2* 1+ ( q )
dup prime? if
dup 7 and dup 1 = swap 7 = or if
2dup 2-exp-mod 1 = if
nip unloop exit
then
then
then drop
loop drop 0 ;
929 factor-mersenne . \ 13007
4423 factor-mersenne . \ 0
Fortran
PROGRAM EXAMPLE
IMPLICIT NONE
INTEGER :: exponent, factor
WRITE(*,*) "Enter exponent of Mersenne number"
READ(*,*) exponent
factor = Mfactor(exponent)
IF (factor == 0) THEN
WRITE(*,*) "No Factor found"
ELSE
WRITE(*,"(A,I0,A,I0)") "M", exponent, " has a factor: ", factor
END IF
CONTAINS
FUNCTION isPrime(number)
! code omitted - see [[Primality by Trial Division]]
END FUNCTION
FUNCTION Mfactor(p)
INTEGER :: Mfactor
INTEGER, INTENT(IN) :: p
INTEGER :: i, k, maxk, msb, n, q
DO i = 30, 0 , -1
IF(BTEST(p, i)) THEN
msb = i
EXIT
END IF
END DO
maxk = 16384 / p ! limit for k to prevent overflow of 32 bit signed integer
DO k = 1, maxk
q = 2*p*k + 1
IF (.NOT. isPrime(q)) CYCLE
IF (MOD(q, 8) /= 1 .AND. MOD(q, 8) /= 7) CYCLE
n = 1
DO i = msb, 0, -1
IF (BTEST(p, i)) THEN
n = MOD(n*n*2, q)
ELSE
n = MOD(n*n, q)
ENDIF
END DO
IF (n == 1) THEN
Mfactor = q
RETURN
END IF
END DO
Mfactor = 0
END FUNCTION
END PROGRAM EXAMPLE
- Output:
M929 has a factor: 13007
FreeBASIC
' FB 1.05.0 Win64
Function isPrime(n As Integer) As Boolean
If n Mod 2 = 0 Then Return n = 2
If n Mod 3 = 0 Then Return n = 3
Dim d As Integer = 5
While d * d <= n
If n Mod d = 0 Then Return False
d += 2
If n Mod d = 0 Then Return False
d += 4
Wend
Return True
End Function
' test 929 plus all prime numbers below 100 which are known not to be Mersenne primes
Dim q(1 To 16) As Integer = {11, 23, 29, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 97, 929}
For k As Integer = 1 To 16
If isPrime(q(k)) Then
Dim As Integer d, i, p, r = q(k)
While r > 0 : r Shl= 1 : Wend
d = 2 * q(k) + 1
Do
i = 1
p = r
While p <> 0
i = (i * i) Mod d
If p < 0 Then i *= 2
If i > d Then i -= d
p Shl= 1
Wend
If i <> 1 Then
d += 2 * q(k)
Else
Exit Do
End If
Loop
Print "2^"; Str(q(k)); Tab(6); " - 1 = 0 (mod"; d; ")"
Else
Print Str(q(k)); " is not prime"
End If
Next
Print
Print "Press any key to quit"
Sleep
- Output:
2^11 - 1 = 0 (mod 23) 2^23 - 1 = 0 (mod 47) 2^29 - 1 = 0 (mod 233) 2^37 - 1 = 0 (mod 223) 2^41 - 1 = 0 (mod 13367) 2^43 - 1 = 0 (mod 431) 2^47 - 1 = 0 (mod 2351) 2^53 - 1 = 0 (mod 6361) 2^59 - 1 = 0 (mod 179951) 2^67 - 1 = 0 (mod 193707721) 2^71 - 1 = 0 (mod 228479) 2^73 - 1 = 0 (mod 439) 2^79 - 1 = 0 (mod 2687) 2^83 - 1 = 0 (mod 167) 2^97 - 1 = 0 (mod 11447) 2^929 - 1 = 0 (mod 13007)
Frink
Frink has built-in routines for iterating through prime numbers and modular exponentiation. The following program will find all of the factors of the number given enough runtime.
for p = primes[]
if modPow[2, 929, p] - 1 == 0
println[p]
- Output:
13007
GAP
MersenneSmallFactor := function(n)
local k, m, d;
if IsPrime(n) then
d := 2*n;
m := 1;
for k in [1 .. 1000000] do
m := m + d;
if PowerModInt(2, n, m) = 1 then
return m;
fi;
od;
fi;
return fail;
end;
# If n is not prime, fail immediately
MersenneSmallFactor(15);
# fail
MersenneSmallFactor(929);
# 13007
MersenneSmallFactor(1009);
# 3454817
# We stop at k = 1000000 in 2*k*n + 1, so it may fail if 2^n - 1 has only larger factors
MersenneSmallFactor(101);
# fail
FactorsInt(2^101-1);
# [ 7432339208719, 341117531003194129 ]
Go
package main
import (
"fmt"
"math"
)
// limit search to small primes. really this is higher than
// you'd want it, but it's fun to factor M67.
const qlimit = 2e8
func main() {
mtest(31)
mtest(67)
mtest(929)
}
func mtest(m int32) {
// the function finds odd prime factors by
// searching no farther than sqrt(N), where N = 2^m-1.
// the first odd prime is 3, 3^2 = 9, so M3 = 7 is still too small.
// M4 = 15 is first number for which test is meaningful.
if m < 4 {
fmt.Printf("%d < 4. M%d not tested.\n", m, m)
return
}
flimit := math.Sqrt(math.Pow(2, float64(m)) - 1)
var qlast int32
if flimit < qlimit {
qlast = int32(flimit)
} else {
qlast = qlimit
}
composite := make([]bool, qlast+1)
sq := int32(math.Sqrt(float64(qlast)))
loop:
for q := int32(3); ; {
if q <= sq {
for i := q * q; i <= qlast; i += q {
composite[i] = true
}
}
if q8 := q % 8; (q8 == 1 || q8 == 7) && modPow(2, m, q) == 1 {
fmt.Printf("M%d has factor %d\n", m, q)
return
}
for {
q += 2
if q > qlast {
break loop
}
if !composite[q] {
break
}
}
}
fmt.Printf("No factors of M%d found.\n", m)
}
// base b to power p, mod m
func modPow(b, p, m int32) int32 {
pow := int64(1)
b64 := int64(b)
m64 := int64(m)
bit := uint(30)
for 1<<bit&p == 0 {
bit--
}
for {
pow *= pow
if 1<<bit&p != 0 {
pow *= b64
}
pow %= m64
if bit == 0 {
break
}
bit--
}
return int32(pow)
}
- Output:
No factors of M31 found. M67 has factor 193707721 M929 has factor 13007
Haskell
Using David Amos module Primes [1] for prime number testing:
import Data.List
import HFM.Primes (isPrime)
import Control.Monad
import Control.Arrow
int2bin = reverse.unfoldr(\x -> if x==0 then Nothing
else Just ((uncurry.flip$(,))$divMod x 2))
trialfac m = take 1. dropWhile ((/=1).(\q -> foldl (((`mod` q).).pm) 1 bs)) $ qs
where qs = filter (liftM2 (&&) (liftM2 (||) (==1) (==7) .(`mod`8)) isPrime ).
map (succ.(2*m*)). enumFromTo 1 $ m `div` 2
bs = int2bin m
pm n b = 2^b*n*n
*Main> trialfac 929
[13007]
Icon and Unicon
The following works in both languages:
procedure main(A)
p := integer(A[1]) | 929
write("M",p," has a factor ",mfactor(p))
end
procedure mfactor(p)
if isPrime(p) then {
limit := sqrt(2^p)-1
k := 1
while 2*p*k-1 < limit do {
q := 2*p*k+1
if isPrime(q) & (q%8 = (1|7)) & btest(p,q) then return q
k +:= 1
}
}
end
procedure btest(p, q)
return (2^p % q) = 1
end
procedure isPrime(n)
if n%(i := 2|3) = 0 then return n = i;
d := 5
while d*d <= n do {
if n%d = 0 then fail
d +:= 2
if n%d = 0 then fail
d +:= 4
}
return
end
Sample runs:
->fmn M929 has a factor 13007 ->fmn 41 M41 has a factor 13367 ->
J
trialfac=: 3 : 0
qs=. (#~8&(1=|+.7=|))(#~1&p:)1+(*(1x+i.@<:@<.)&.-:)y
qs#~1=qs&|@(2&^@[**:@])/ 1,~ |.#: y
)
- Examples:
trialfac 929
13007
trialfac 44497
Empty output --> No factors found.
Java
import java.math.BigInteger;
class MersenneFactorCheck
{
private final static BigInteger TWO = BigInteger.valueOf(2);
public static boolean isPrime(long n)
{
if (n == 2)
return true;
if ((n < 2) || ((n & 1) == 0))
return false;
long maxFactor = (long)Math.sqrt((double)n);
for (long possibleFactor = 3; possibleFactor <= maxFactor; possibleFactor += 2)
if ((n % possibleFactor) == 0)
return false;
return true;
}
public static BigInteger findFactorMersenneNumber(int primeP)
{
if (primeP <= 0)
throw new IllegalArgumentException();
BigInteger bigP = BigInteger.valueOf(primeP);
BigInteger m = BigInteger.ONE.shiftLeft(primeP).subtract(BigInteger.ONE);
// There are more complicated ways of getting closer to sqrt(), but not that important here, so go with simple
BigInteger maxFactor = BigInteger.ONE.shiftLeft((primeP + 1) >>> 1);
BigInteger twoP = BigInteger.valueOf(primeP << 1);
BigInteger possibleFactor = BigInteger.ONE;
int possibleFactorBits12 = 0;
int twoPBits12 = primeP & 3;
while ((possibleFactor = possibleFactor.add(twoP)).compareTo(maxFactor) <= 0)
{
possibleFactorBits12 = (possibleFactorBits12 + twoPBits12) & 3;
// "Furthermore, q must be 1 or 7 mod 8". We know it's odd due to the +1 done above, so bit 0 is set. Therefore, we only care about bits 1 and 2 equaling 00 or 11
if ((possibleFactorBits12 == 0) || (possibleFactorBits12 == 3))
if (TWO.modPow(bigP, possibleFactor).equals(BigInteger.ONE))
return possibleFactor;
}
return null;
}
public static void checkMersenneNumber(int p)
{
if (!isPrime(p))
{
System.out.println("M" + p + " is not prime");
return;
}
BigInteger factor = findFactorMersenneNumber(p);
if (factor == null)
System.out.println("M" + p + " is prime");
else
System.out.println("M" + p + " is not prime, has factor " + factor);
return;
}
public static void main(String[] args)
{
for (int p = 1; p <= 50; p++)
checkMersenneNumber(p);
checkMersenneNumber(929);
return;
}
}
- Output:
M1 is not prime M2 is prime M3 is prime M4 is not prime M5 is prime M6 is not prime M7 is prime M8 is not prime M9 is not prime M10 is not prime M11 is not prime, has factor 23 M12 is not prime M13 is prime M14 is not prime ... M47 is not prime, has factor 2351 M48 is not prime M49 is not prime M50 is not prime M929 is not prime, has factor 13007
JavaScript
function mersenne_factor(p){
var limit, k, q
limit = Math.sqrt(Math.pow(2,p) - 1)
k = 1
while ((2*k*p - 1) < limit){
q = 2*k*p + 1
if (isPrime(q) && (q % 8 == 1 || q % 8 == 7) && trial_factor(2,p,q)){
return q // q is a factor of 2**p-1
}
k++
}
return null
}
function isPrime(value){
for (var i=2; i < value; i++){
if (value % i == 0){
return false
}
if (value % i != 0){
return true;
}
}
}
function trial_factor(base, exp, mod){
var square, bits
square = 1
bits = exp.toString(2).split('')
for (var i=0,ln=bits.length; i<ln; i++){
square = Math.pow(square, 2) * (bits[i] == 1 ? base : 1) % mod
}
return (square == 1)
}
function check_mersenne(p){
var f, result
console.log("M"+p+" = 2^"+p+"-1 is ")
f = mersenne_factor(p)
console.log(f == null ? "prime" : "composite with factor "+f)
}
> check_mersenne(3) "M3 = 2**3-1 is prime" > check_mersenne(23) "M23 = 2**23-1 is composite with factor 47" > check_mersenne(929) "M929 = 2**929-1 is composite with factor 13007"
Julia
# v0.6
using Primes
function mersennefactor(p::Int)::Int
q = 2p + 1
while true
if log2(q) > p / 2
return -1
elseif q % 8 in (1, 7) && Primes.isprime(q) && powermod(2, p, q) == 1
return q
end
q += 2p
end
end
for i in filter(Primes.isprime, push!(collect(1:20), 929))
mf = mersennefactor(i)
if mf != -1 println("M$i = ", mf, " × ", (big(2) ^ i - 1) ÷ mf)
else println("M$i is prime") end
end
- Output:
M2 is prime M3 is prime M5 is prime M7 is prime M11 = 23 × 89 M13 is prime M17 is prime M19 is prime M929 = 13007 × 34889024892493825975045478116339093030512026953827804293400962134 88946572057852012474541189660261508521493994102599382170621001921687473524507195 61908445272675574320888385228421992652298715687625495638077382028762529439880103 124705348782610789919949159935587158612289264184273
Kotlin
// version 1.0.6
fun isPrime(n: Int): Boolean {
if (n < 2) return false
if (n % 2 == 0) return n == 2
if (n % 3 == 0) return n == 3
var d = 5
while (d * d <= n) {
if (n % d == 0) return false
d += 2
if (n % d == 0) return false
d += 4
}
return true
}
fun main(args: Array<String>) {
// test 929 plus all prime numbers below 100 which are known not to be Mersenne primes
val q = intArrayOf(11, 23, 29, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 97, 929)
for (k in 0 until q.size) {
if (isPrime(q[k])) {
var i: Long
var d: Int
var p: Int
var r: Int = q[k]
while (r > 0) r = r shl 1
d = 2 * q[k] + 1
while (true) {
i = 1L
p = r
while (p != 0) {
i = (i * i) % d
if (p < 0) i *= 2
if (i > d) i -= d
p = p shl 1
}
if (i != 1L)
d += 2 * q[k]
else
break
}
println("2^${"%3d".format(q[k])} - 1 = 0 (mod $d)")
} else {
println("${q[k]} is not prime")
}
}
}
- Output:
2^ 11 - 1 = 0 (mod 23) 2^ 23 - 1 = 0 (mod 47) 2^ 29 - 1 = 0 (mod 233) 2^ 37 - 1 = 0 (mod 223) 2^ 41 - 1 = 0 (mod 13367) 2^ 43 - 1 = 0 (mod 431) 2^ 47 - 1 = 0 (mod 2351) 2^ 53 - 1 = 0 (mod 6361) 2^ 59 - 1 = 0 (mod 179951) 2^ 67 - 1 = 0 (mod 193707721) 2^ 71 - 1 = 0 (mod 228479) 2^ 73 - 1 = 0 (mod 439) 2^ 79 - 1 = 0 (mod 2687) 2^ 83 - 1 = 0 (mod 167) 2^ 97 - 1 = 0 (mod 11447) 2^929 - 1 = 0 (mod 13007)
Lingo
on modPow (b, e, m)
bits = getBits(e)
sq = 1
repeat while TRUE
tb = bits[1]
bits.deleteAt(1)
sq = sq*sq
if tb then sq=sq*b
sq = sq mod m
if bits.count=0 then return sq
end repeat
end
on getBits (n)
bits = []
f = 1
repeat while TRUE
bits.addAt(1, bitAnd(f, n)>0)
f = f * 2
if f>n then exit repeat
end repeat
return bits
end
repeat with i = 2 to the maxInteger
if modPow(2, 929, i)=1 then
put "M929 has a factor: " & i
exit repeat
end if
end repeat
- Output:
-- "M929 has a factor: 13007"
Mathematica/Wolfram Language
Believe it or not, this type of test runs faster in Mathematica than the squaring version described above.
For[i = 2, i < Prime[1000000], i = NextPrime[i],
If[Mod[2^44497, i] == 1,
Print["divisible by "<>i]]]; Print["prime test passed; call Lucas and Lehmer"]
Maxima
mersenne_fac(p) := block([m: 2^p - 1, k: 1],
while mod(m, 2 * k * p + 1) # 0 do k: k + 1,
2 * k * p + 1
)$
mersenne_fac(929);
/* 13007 */
Nim
import math
proc isPrime(a: int): bool =
if a == 2: return true
if a < 2 or a mod 2 == 0: return false
for i in countup(3, int sqrt(float a), 2):
if a mod i == 0:
return false
return true
const q = 929
if not isPrime q: quit 1
var r = q
while r > 0: r = r shl 1
var d = 2 * q + 1
while true:
var i = 1
var p = r
while p != 0:
i = (i * i) mod d
if p < 0: i *= 2
if i > d: i -= d
p = p shl 1
if i != 1: d += 2 * q
else: break
echo "2^",q," - 1 = 0 (mod ",d,")"
- Output:
2^929 - 1 = 0 (mod 13007)
Octave
(GNU Octave has a isprime
built-in test)
% test a bit; lsb is 1 (like built-in bit* ops)
function b = bittst(n, p)
b = bitand(n, 2^(p-1)) > 0;
endfunction
function f = Mfactor(p)
% msb is the index of the first non-zero bit
[b, msb] = max(bitand(p, 2 .^ [32:-1:1]) > 0);
maxk = floor(sqrt(intmax()) / p);
for k = 1 : maxk
q = 2*p*k + 1;
if ( ! isprime(q) )
continue;
endif
if ( (mod(q, 8) != 1) && ( mod(q, 8) != 7) )
continue;
endif
n = 1;
for i = msb:-1:1
if ( bittst(p, i) )
n = mod(n*n*2, q);
else
n = mod(n*n, q);
endif
endfor
if ( n==1 )
f = q;
return
endif
endfor
f = 0;
endfunction
printf("%d\n", Mfactor(929));
PARI/GP
This version takes about 15 microseconds to find a factor of 2929 − 1.
factorMersenne(p)={
forstep(q=2*p+1,sqrt(2)<<(p\2),2*p,
[1,0,0,0,0,0,1][q%8] && Mod(2, q)^p==1 && return(q)
);
1<<p-1
};
factorMersenne(929)
This implementation seems to be broken:
TM(p) = local(status=1, i=1, len=0, S=0);{
printp("Test TM \t...");
S=2*p+1;
len = length(binary(p));
B = Vecsmall(binary(p));
q = B[i]*B[i];
while( i<=len & status ==1,
if( B[i] != 0,
q = q*2;
);
r = q%S;
q = r*r;
if( i == len & r == 1,
status = 0;
printp("Not Prime!");
);
i++;
);
return(status);
}
Pascal
program FactorsMersenneNumber(input, output);
function isPrime(n: longint): boolean;
var
d: longint;
begin
isPrime := true;
if (n mod 2) = 0 then
begin
isPrime := (n = 2);
exit;
end;
if (n mod 3) = 0 then
begin
isPrime := (n = 3);
exit;
end;
d := 5;
while d*d <= n do
begin
if (n mod d) = 0 then
begin
isPrime := false;
exit;
end;
d := d + 2;
end;
end;
function btest(n, pos: longint): boolean;
begin
btest := (n shr pos) mod 2 = 1;
end;
function MFactor(p: longint): longint;
var
i, k, maxk, msb, n, q: longint;
begin
for i := 30 downto 0 do
if btest(p, i) then
begin
msb := i;
break;
end;
maxk := 16384 div p; // limit for k to prevent overflow of 32 bit signed integer
for k := 1 to maxk do
begin
q := 2*p*k + 1;
if not isprime(q) then
continue;
if ((q mod 8) <> 1) and ((q mod 8) <> 7) then
continue;
n := 1;
for i := msb downto 0 do
if btest(p, i) then
n := (n*n*2) mod q
else
n := (n*n) mod q;
if n = 1 then
begin
mfactor := q;
exit;
end;
end;
mfactor := 0;
end;
var
exponent, factor: longint;
begin
write('Enter the exponent of the Mersenne number (suggestion: 929): ');
readln(exponent);
if not isPrime(exponent) then
begin
writeln('M', exponent, ' (2**', exponent, ' - 1) is not prime.');
exit;
end;
factor := MFactor(exponent);
if factor = 0 then
writeln('M', exponent, ' (2**', exponent, ' - 1) has no factor.')
else
writeln('M', exponent, ' (2**', exponent, ' - 1) has the factor: ', factor);
end.
- Output:
:> ./FactorsMersenneNumber Enter the exponent of the Mersenne number (suggestion: 929): 929 M929 (2**929 - 1) has the factor: 13007
Perl
use strict;
use utf8;
sub factors {
my $n = shift;
my $p = 2;
my @out;
while ($n >= $p * $p) {
while ($n % $p == 0) {
push @out, $p;
$n /= $p;
}
$p = next_prime($p);
}
push @out, $n if $n > 1 || !@out;
@out;
}
sub next_prime {
my $p = shift;
do { $p = $p == 2 ? 3 : $p + 2 } until is_prime($p);
$p;
}
my %pcache;
sub is_prime {
my $x = shift;
$pcache{$x} //= (factors($x) == 1)
}
sub mtest {
my @bits = split "", sprintf("%b", shift);
my $p = shift;
my $sq = 1;
while (@bits) {
$sq = $sq * $sq;
$sq *= 2 if shift @bits;
$sq %= $p;
}
$sq == 1;
}
for my $m (2 .. 60, 929) {
next unless is_prime($m);
use bigint;
my ($f, $k, $x) = (0, 0, 2**$m - 1);
my $q;
while (++$k) {
$q = 2 * $k * $m + 1;
next if (($q & 7) != 1 && ($q & 7) != 7);
next unless is_prime($q);
last if $q * $q > $x;
last if $f = mtest($m, $q);
}
print $f? "M$m = $x = $q × @{[$x / $q]}\n"
: "M$m = $x is prime\n";
}
- Output:
M2 = 3 is prime M2 = 3 is prime M3 = 7 is prime M5 = 31 is prime M7 = 127 is prime M11 = 2047 = 23 × 89 M13 = 8191 is prime ... M53 = 9007199254740991 = 6361 × 1416003655831 M59 = 576460752303423487 = 179951 × 3203431780337 M929 = 4538..<yadda yadda>..8911 = 13007 × 348890..<blah blah>..84273
Following the task introduction, this uses GMP's modular exponentiation and simple probable prime test for the small numbers, then looks for small factors before doing a Lucas-Lehmer test. For ranges above about p=2000, looking for small factors this way saves time (the amount of testing should be adjusted based on the input size and platform -- this example just uses a fixed amount). Note as well that the Lucas-Lehmer test shown here is ignoring the large speedup we can get by optimizing the modulo operation, but that's a different task.
use Math::GMP;
# Use GMP's simple probable prime test.
sub is_prime { Math::GMP->new(shift)->probab_prime(20); }
# Lucas-Lehmer test, deterministic for 2^p-1 given p
sub is_mersenne_prime {
my($p, $mp, $s) = ($_[0], Math::GMP->new(2)**$_[0]-1, Math::GMP->new(4));
return 1 if $p == 2;
$s = ($s * $s - 2) % $mp for 3 .. $p;
$s == 0;
}
for my $p (2 .. 100, 929) {
next unless is_prime($p);
my $mp = Math::GMP->new(2) ** $p - 1;
my $lim = $mp->bsqrt();
$lim = 1000000 if $lim > 1000000; # We're using it as a pre-test
my $factor;
for (my $q = Math::GMP->new(2*$p+1); $q <= $lim && !$factor; $q += 2*$p) {
next unless ($q & 7) == 1 || ($q & 7) == 7;
next unless is_prime($q);
$factor = $q if Math::GMP->new(2)->powm_gmp($p,$q) == 1; # $mp % $q == 0
}
if ($factor) {
print "M$p = $factor * ",$mp/$factor,"\n";
} else {
print "M$p is ", is_mersenne_prime($p) ? "prime" : "composite", "\n";
}
}
- Output:
M2 is prime M3 is prime M5 is prime M7 is prime M11 = 23 * 89 M13 is prime M17 is prime M19 is prime M23 = 47 * 178481 M29 = 233 * 2304167 M31 is prime M37 = 223 * 616318177 M41 = 13367 * 164511353 M43 = 431 * 20408568497 M47 = 2351 * 59862819377 M53 = 6361 * 1416003655831 M59 = 179951 * 3203431780337 M61 is prime M67 is composite M71 = 228479 * 10334355636337793 M73 = 439 * 21514198099633918969 M79 = 2687 * 224958284260258499201 M83 = 167 * 57912614113275649087721 M89 is prime M97 = 11447 * 13842607235828485645766393 M929 = 13007 * 348890248924[.....]64184273
Phix
Translation/Amalgamation of BBC BASIC, D, and Go
with javascript_semantics function modpow(atom x, atom n, atom m) atom i = n, y = 1, z = x while i do if and_bits(i,1) then y = mod(y*z,m) end if z = mod(z*z,m) i = floor(i/2) end while return y end function function mersenne_factor(integer p) if not is_prime(p) then return -1 end if atom limit = sqrt(power(2,p))-1 integer k = 1 while 1 do atom q = 2*p*k + 1 if q>=limit then exit end if if find(mod(q,8),{1,7}) and is_prime(q) and modpow(2,p,q)=1 then return q end if k += 1 end while return 0 end function sequence tests = {11, 23, 29, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 97, 929, 937} for i=1 to length(tests) do integer ti = tests[i] printf(1,"A factor of M%d is %d\n",{ti,mersenne_factor(ti)}) end for
- Output:
A factor of M11 is 23 A factor of M23 is 47 A factor of M29 is 233 A factor of M37 is 223 A factor of M41 is 13367 A factor of M43 is 431 A factor of M47 is 2351 A factor of M53 is 6361 A factor of M59 is 179951 A factor of M67 is 193707721 A factor of M71 is 228479 A factor of M73 is 439 A factor of M79 is 2687 A factor of M83 is 167 A factor of M97 is 11447 A factor of M929 is 13007 A factor of M937 is 28111
PHP
Requires bcmath
echo 'M929 has a factor: ', mersenneFactor(929), '</br>';
function mersenneFactor($p) {
$limit = sqrt(pow(2, $p) - 1);
for ($k = 1; 2 * $p * $k - 1 < $limit; $k++) {
$q = 2 * $p * $k + 1;
if (isPrime($q) && ($q % 8 == 1 || $q % 8 == 7) && bcpowmod("2", "$p", "$q") == "1") {
return $q;
}
}
return 0;
}
function isPrime($n) {
if ($n < 2 || $n % 2 == 0) return $n == 2;
for ($i = 3; $i * $i <= $n; $i += 2) {
if ($n % $i == 0) {
return false;
}
}
return true;
}
- Output:
M929 has a factor: 13007
PicoLisp
(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
(de prime? (N)
(or
(= N 2)
(and
(> N 1)
(bit? 1 N)
(let S (sqrt N)
(for (D 3 T (+ D 2))
(T (> D S) T)
(T (=0 (% N D)) NIL) ) ) ) ) )
(de mFactor (P)
(let (Lim (sqrt (dec (** 2 P))) K 0 Q)
(loop
(setq Q (inc (* 2 (inc 'K) P)))
(T (>= Q Lim) NIL)
(T
(and
(member (% Q 8) (1 7))
(prime? Q)
(= 1 (**Mod 2 P Q)) )
Q ) ) ) )
- Output:
: (for P (2 3 4 5 7 11 13 17 19 23 29 31 37 41 43 47 53 929) (prinl "M" P " = 2**" P "-1 is " (cond ((not (prime? P)) "not prime") ((mFactor P) (pack "composite with factor " @)) (T "prime") ) ) ) M2 = 2**2-1 is prime M3 = 2**3-1 is prime M4 = 2**4-1 is not prime M5 = 2**5-1 is prime M7 = 2**7-1 is prime M11 = 2**11-1 is composite with factor 23 M13 = 2**13-1 is prime M17 = 2**17-1 is prime M19 = 2**19-1 is prime M23 = 2**23-1 is composite with factor 47 M29 = 2**29-1 is composite with factor 233 M31 = 2**31-1 is prime M37 = 2**37-1 is composite with factor 223 M41 = 2**41-1 is composite with factor 13367 M43 = 2**43-1 is composite with factor 431 M47 = 2**47-1 is composite with factor 2351 M53 = 2**53-1 is composite with factor 6361 M929 = 2**929-1 is composite with factor 13007
Prolog
mersenne_factor(P, F) :-
prime(P),
once((
between(1, 100_000, K), % Fail if we can't find a small factor
Q is 2*K*P + 1,
test_factor(Q, P, F))).
test_factor(Q, P, prime) :- Q*Q > (1 << P - 1), !.
test_factor(Q, P, Q) :-
R is Q /\ 7, member(R, [1, 7]),
prime(Q),
powm(2, P, Q) =:= 1.
wheel235(L) :-
W = [4, 2, 4, 2, 4, 6, 2, 6 | W],
L = [1, 2, 2 | W].
prime(N) :-
N >= 2,
wheel235(W),
prime(N, 2, W).
prime(N, D, _) :- D*D > N, !.
prime(N, D, [A|As]) :-
N mod D =\= 0,
D2 is D + A, prime(N, D2, As).
- Output:
?- mersenne_factor(23, X). X = 47. ?- mersenne_factor(5,X). X = prime. ?- mersenne_factor(25,X). false. ?- mersenne_factor(929,X). X = 13007. ?- mersenne_factor(127,X). false.
Python
def is_prime(number):
return True # code omitted - see Primality by Trial Division
def m_factor(p):
max_k = 16384 / p # arbitrary limit; since Python automatically uses long's, it doesn't overflow
for k in xrange(max_k):
q = 2*p*k + 1
if not is_prime(q):
continue
elif q % 8 != 1 and q % 8 != 7:
continue
elif pow(2, p, q) == 1:
return q
return None
if __name__ == '__main__':
exponent = int(raw_input("Enter exponent of Mersenne number: "))
if not is_prime(exponent):
print "Exponent is not prime: %d" % exponent
else:
factor = m_factor(exponent)
if not factor:
print "No factor found for M%d" % exponent
else:
print "M%d has a factor: %d" % (exponent, factor)
- Example:
Enter exponent of Mersenne number: 929 M929 has a factor: 13007
Racket
#lang racket
(define (number->digits n)
(map (compose1 string->number string)
(string->list (number->string n 2))))
(define (modpow exp base)
(for/fold ([square 1])
([d (number->digits exp)])
(modulo (* (if (= d 1) 2 1) square square) base)))
; Search through all integers from 1 on to find the first divisor.
; Returns #f if 2^p-1 is prime.
(define (mersenne-factor p)
(for/first ([i (in-range 1 (floor (expt 2 (quotient p 2))) (* 2 p))]
#:when (and (member (modulo i 8) '(1 7))
(= 1 (modpow p i))))
i))
(mersenne-factor 929)
- Output:
13007
Raku
(formerly Perl 6)
sub mtest($bits, $p) {
my @bits = $bits.base(2).comb;
loop (my $sq = 1; @bits; $sq %= $p) {
$sq ×= $sq;
$sq += $sq if 1 == @bits.shift;
}
$sq == 1;
}
for flat 2 .. 60, 929 -> $m {
next unless is-prime($m);
my $f = 0;
my $x = 2**$m - 1;
my $q;
for 1..* -> $k {
$q = 2 × $k × $m + 1;
next unless $q % 8 == 1|7 or is-prime($q);
last if $q × $q > $x or $f = mtest($m, $q);
}
say $f ?? "M$m = $x\n\t= $q × { $x div $q }"
!! "M$m = $x is prime";
}
- Output:
M2 = 3 is prime M3 = 7 is prime M5 = 31 is prime M7 = 127 is prime M11 = 2047 = 23 × 89 M13 = 8191 is prime M17 = 131071 is prime M19 = 524287 is prime M23 = 8388607 = 47 × 178481 M29 = 536870911 = 233 × 2304167 M31 = 2147483647 is prime M37 = 137438953471 = 223 × 616318177 M41 = 2199023255551 = 13367 × 164511353 M43 = 8796093022207 = 431 × 20408568497 M47 = 140737488355327 = 2351 × 59862819377 M53 = 9007199254740991 = 6361 × 1416003655831 M59 = 576460752303423487 = 179951 × 3203431780337 M929 = 4538015467766671944574165338592225830478699345884382504442663144885072806275648112625635725391102144133907238129251016389326737199538896813326509341743147661691195191795226666084858428449394948944821764472508048114220424520501343042471615418544488778723282182172070046459244838911 = 13007 × 348890248924938259750454781163390930305120269538278042934009621348894657205785201247454118966026150852149399410259938217062100192168747352450719561908445272675574320888385228421992652298715687625495638077382028762529439880103124705348782610789919949159935587158612289264184273
REXX
REXX practically has no limit (well, up to around 8 million) on the number of decimal digits (precision).
This REXX version automatically adjusts the numeric digits to whatever is needed.
/*REXX program uses exponent─and─mod operator to test possible Mersenne numbers. */
numeric digits 20 /*this will be increased if necessary. */
parse arg N spec /*obtain optional arguments from the CL*/
if N=='' | N=="," then N= 88 /*Not specified? Then use the default.*/
if spec=='' | spec=="," then spec= 920 970 /* " " " " " " */
do j=1; z= j /*process a range, & then do some more.*/
if j==N then j= word(spec, 1) /*now, use the high range of numbers. */
if j>word(spec, 2) then leave /*done with " " " " " */
if \isPrime(z) then iterate /*if Z isn't a prime, keep plugging.*/
r= commas( testMer(z) ); L= length(r) /*add commas; get its new length. */
if r==0 then say right('M'z, 10) "──────── is a Mersenne prime."
else say right('M'z, 50) "is composite, a factor:"right(r, max(L, 13) )
end /*j*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do jc=length(_)-3 to 1 by -3; _=insert(',', _, jc); end; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
isPrime: procedure; parse arg x; if wordpos(x, '2 3 5 7') \== 0 then return 1
if x<11 then return 0; if x//2 == 0 | x//3 == 0 then return 0
do j=5 by 6; if x//j == 0 | x//(j+2) == 0 then return 0
if j*j>x then return 1 /*◄─┐ ___ */
end /*j*/ /* └─◄ Is j>√ x ? Then return 1*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
iSqrt: procedure; parse arg x; #= 1; r= 0; do while #<=x; #= # * 4
end /*while*/
do while #>1; #= # % 4; _= x-r-#; r= r % 2
if _>=0 then do; x= _; r= r + #
end
end /*while*/ /*iSqrt ≡ integer square root.*/
return r /*───── ─ ── ─ ─ */
/*──────────────────────────────────────────────────────────────────────────────────────*/
testMer: procedure; parse arg x; p= 2**x /* [↓] do we have enough digits?*/
$$=x2b( d2x(x) ) + 0
if pos('E',p)\==0 then do; parse var p "E" _; numeric digits _ + 2; p= 2**x
end
!.= 1; !.1= 0; !.7= 0 /*array used for a quicker test. */
R= iSqrt(p) /*obtain integer square root of P*/
do k=2 by 2; q= k*x + 1 /*(shortcut) compute value of Q. */
m= q // 8 /*obtain the remainder when ÷ 8. */
if !.m then iterate /*M must be either one or seven.*/
parse var q '' -1 _; if _==5 then iterate /*last digit a five ? */
if q// 3==0 then iterate /*divisible by three? */
if q// 7==0 then iterate /* " " seven? */
if q//11==0 then iterate /* " " eleven?*/
/* ____ */
if q>R then return 0 /*Is q>√2**x ? A Mersenne prime*/
sq= 1; $= $$ /*obtain binary version from $. */
do until $==''; sq= sq*sq
parse var $ _ 2 $ /*obtain 1st digit and the rest. */
if _ then sq= (sq+sq) // q
end /*until*/
if sq==1 then return q /*Not a prime? Return a factor.*/
end /*k*/
Program note: the iSqrt function computes the integer square root of a non-negative integer without using any floating point, just integers.
- output when using the default (two) ranges of numbers:
M2 ──────── is a Mersenne prime. M3 ──────── is a Mersenne prime. M5 ──────── is a Mersenne prime. M7 ──────── is a Mersenne prime. M11 is composite, a factor: 23 M13 ──────── is a Mersenne prime. M17 ──────── is a Mersenne prime. M19 ──────── is a Mersenne prime. M23 is composite, a factor: 47 M29 is composite, a factor: 233 M31 ──────── is a Mersenne prime. M37 is composite, a factor: 223 M41 is composite, a factor: 13,367 M43 is composite, a factor: 431 M47 is composite, a factor: 2,351 M53 is composite, a factor: 6,361 M59 is composite, a factor: 179,951 M61 ──────── is a Mersenne prime. M67 is composite, a factor: 193,707,721 M71 is composite, a factor: 228,479 M73 is composite, a factor: 439 M79 is composite, a factor: 2,687 M83 is composite, a factor: 167 M929 is composite, a factor: 13,007 M937 is composite, a factor: 28,111 M941 is composite, a factor: 7,529 M947 is composite, a factor: 295,130,657 M953 is composite, a factor: 343,081 M967 is composite, a factor: 23,209
Ring
# Project : Factors of a Mersenne number
see "A factor of M929 is " + mersennefactor(929) + nl
see "A factor of M937 is " + mersennefactor(937) + nl
func mersennefactor(p)
if not isprime(p)
return -1
ok
for k = 1 to 50
q = 2*k*p + 1
if (q && 7) = 1 or (q && 7) = 7
if isprime(q)
if modpow(2, p, q) = 1
return q
ok
ok
ok
next
return 0
func isprime(num)
if (num <= 1) return 0 ok
if (num % 2 = 0) and num != 2 return 0 ok
for i = 3 to floor(num / 2) -1 step 2
if (num % i = 0) return 0 ok
next
return 1
func modpow(x,n,m)
i = n
y = 1
z = x
while i > 0
if i & 1
y = (y * z) % m
ok
z = (z * z) % m
i = (i >> 1)
end
return y
Output:
A factor of M929 is 13007 A factor of M937 is 28111
RPL
PRIM?
is defined at Primality by trial division
RPL code | Comment |
---|---|
≪ SWAP R→B → quotient power ≪ 2 power B→R LN 2 LN / FLOOR ^ R→B 1 WHILE OVER B→R REPEAT SQ IF OVER power AND B→R THEN DUP + END quotient MOD SWAP SR SWAP END SWAP DROP ≫ ≫ 'MODPOW' STO ≪ 2 OVER ^ 1 - √ 0 → power max k ≪ 1 WHILE 'k' INCR 2 * 1 + DUP max ≤ REPEAT IF { 1 7 } OVER 8 MOD POS THEN IF DUP PRIM? THEN IF power OVER MODPOW 1 == THEN SWAP max 'k' STO END END END DROP END DROP ≫ 'MFACT' STO |
MODPOW ( power quotient → remainder ) create top-bit mask square = 1 while mask is not zero square *= square if unmasked bit = 1 then square += square square = square mod quotient shift mask right clean stack return square MFACT ( N → factor ) factor = 1 while 2k+1 ≤ sqrt(M(N)) if 2k+1 mod 8 is equal to 1 or 7 if 2k+1 prime is 2K+1 a factor of M(N) ? if yes, exit loop return factor |
929 MFACT
- Output:
1: 13007
Factor found in 69 minutes on a 4-bit HP-48SX calculator.
Ruby
require 'prime'
def mersenne_factor(p)
limit = Math.sqrt(2**p - 1)
k = 1
while (2*k*p - 1) < limit
q = 2*k*p + 1
if q.prime? and (q % 8 == 1 or q % 8 == 7) and trial_factor(2,p,q)
# q is a factor of 2**p-1
return q
end
k += 1
end
nil
end
def trial_factor(base, exp, mod)
square = 1
("%b" % exp).each_char {|bit| square = square**2 * (bit == "1" ? base : 1) % mod}
(square == 1)
end
def check_mersenne(p)
print "M#{p} = 2**#{p}-1 is "
f = mersenne_factor(p)
if f.nil?
puts "prime"
else
puts "composite with factor #{f}"
end
end
Prime.each(53) { |p| check_mersenne p }
check_mersenne 929
- Output:
M2 = 2**2-1 is prime M3 = 2**3-1 is prime M5 = 2**5-1 is prime M7 = 2**7-1 is prime M11 = 2**11-1 is composite with factor 23 M13 = 2**13-1 is prime M17 = 2**17-1 is prime M19 = 2**19-1 is prime M23 = 2**23-1 is composite with factor 47 M29 = 2**29-1 is composite with factor 233 M31 = 2**31-1 is prime M37 = 2**37-1 is composite with factor 223 M41 = 2**41-1 is composite with factor 13367 M43 = 2**43-1 is composite with factor 431 M47 = 2**47-1 is composite with factor 2351 M53 = 2**53-1 is composite with factor 6361 M929 = 2**929-1 is composite with factor 13007
Rust
fn bit_count(mut n: usize) -> usize {
let mut count = 0;
while n > 0 {
n >>= 1;
count += 1;
}
count
}
fn mod_pow(p: usize, n: usize) -> usize {
let mut square = 1;
let mut bits = bit_count(p);
while bits > 0 {
square = square * square;
bits -= 1;
if (p & (1 << bits)) != 0 {
square <<= 1;
}
square %= n;
}
return square;
}
fn is_prime(n: usize) -> bool {
if n < 2 {
return false;
}
if n % 2 == 0 {
return n == 2;
}
if n % 3 == 0 {
return n == 3;
}
let mut p = 5;
while p * p <= n {
if n % p == 0 {
return false;
}
p += 2;
if n % p == 0 {
return false;
}
p += 4;
}
true
}
fn find_mersenne_factor(p: usize) -> usize {
let mut k = 0;
loop {
k += 1;
let q = 2 * k * p + 1;
if q % 8 == 1 || q % 8 == 7 {
if mod_pow(p, q) == 1 && is_prime(p) {
return q;
}
}
}
}
fn main() {
println!("{}", find_mersenne_factor(929));
}
- Output:
13007
Scala
Full-blown version
/** Find factors of a Mersenne number
*
* The implementation finds factors for M929 and further.
*
* @example M59 = 2^059 - 1 = 576460752303423487 ( 2 msec)
* @example = 179951 × 3203431780337.
*/
object FactorsOfAMersenneNumber extends App {
val two: BigInt = 2
// An infinite stream of primes, lazy evaluation and memo-ized
val oddPrimes = sieve(LazyList.from(3, 2))
def sieve(nums: LazyList[Int]): LazyList[Int] =
LazyList.cons(nums.head, sieve((nums.tail) filter (_ % nums.head != 0)))
def primes: LazyList[Int] = sieve(2 #:: oddPrimes)
def factorMersenne(p: Int): Option[Long] = {
val limit = (mersenne(p) - 1 min Int.MaxValue).toLong
def factorTest(p: Long, q: Long): Boolean = {
(List(1, 7) contains (q % 8)) && two.modPow(p, q) == 1 && BigInt(q).isProbablePrime(7)
}
// Build a stream of factors from (2*p+1) step-by (2*p)
def s(a: Long): LazyList[Long] = a #:: s(a + (2 * p)) // Build stream of possible factors
// Limit and Filter Stream and then take the head element
val e = s(2 * p + 1).takeWhile(_ < limit).filter(factorTest(p, _))
e.headOption
}
def mersenne(p: Int): BigInt = (two pow p) - 1
// Test
(primes takeWhile (_ <= 97)) ++ List(929, 937) foreach { p => { // Needs some intermediate results for nice formatting
val nMersenne = mersenne(p);
val lit = s"${nMersenne}"
val preAmble = f"${s"M${p}"}%4s = 2^$p%03d - 1 = ${lit}%s"
val datum = System.nanoTime
val result = factorMersenne(p)
val mSec = ((System.nanoTime - datum) / 1.0e+6).round
def decStr = {
if (lit.length > 30) f"(M has ${lit.length}%3d dec)" else ""
}
def sPrime: String = {
if (result.isEmpty) " is a Mersenne prime number." else " " * 28
}
println(f"$preAmble${sPrime} ${f"($mSec%,1d"}%13s msec)")
if (result.isDefined)
println(f"${decStr}%-17s = ${result.get} × ${nMersenne / result.get}")
}
}
}
- Output:
M2 = 2^002 - 1 = 3 is a Mersenne prime number. (63 msec) M3 = 2^003 - 1 = 7 is a Mersenne prime number. (0 msec) M5 = 2^005 - 1 = 31 is a Mersenne prime number. (1 msec) M7 = 2^007 - 1 = 127 is a Mersenne prime number. (2 msec) M11 = 2^011 - 1 = 2047 (2.097 msec) = 23 × 89 M13 = 2^013 - 1 = 8191 is a Mersenne prime number. (33 msec) M17 = 2^017 - 1 = 131071 is a Mersenne prime number. (254 msec) M19 = 2^019 - 1 = 524287 is a Mersenne prime number. (524 msec) M23 = 2^023 - 1 = 8388607 (0 msec) = 47 × 178481 M29 = 2^029 - 1 = 536870911 (0 msec) = 233 × 2304167 M31 = 2^031 - 1 = 2147483647 is a Mersenne prime number. (31.484 msec) M37 = 2^037 - 1 = 137438953471 (0 msec) = 223 × 616318177 M41 = 2^041 - 1 = 2199023255551 (0 msec) = 13367 × 164511353 M43 = 2^043 - 1 = 8796093022207 (0 msec) = 431 × 20408568497 M47 = 2^047 - 1 = 140737488355327 (0 msec) = 2351 × 59862819377 M53 = 2^053 - 1 = 9007199254740991 (0 msec) = 6361 × 1416003655831 M59 = 2^059 - 1 = 576460752303423487 (1 msec) = 179951 × 3203431780337 M61 = 2^061 - 1 = 2305843009213693951 is a Mersenne prime number. (16.756 msec) M67 = 2^067 - 1 = 147573952589676412927 (1.435 msec) = 193707721 × 761838257287 M71 = 2^071 - 1 = 2361183241434822606847 (2 msec) = 228479 × 10334355636337793 M73 = 2^073 - 1 = 9444732965739290427391 (0 msec) = 439 × 21514198099633918969 M79 = 2^079 - 1 = 604462909807314587353087 (0 msec) = 2687 × 224958284260258499201 M83 = 2^083 - 1 = 9671406556917033397649407 (0 msec) = 167 × 57912614113275649087721 M89 = 2^089 - 1 = 618970019642690137449562111 is a Mersenne prime number. (11.097 msec) M97 = 2^097 - 1 = 158456325028528675187087900671 (0 msec) = 11447 × 13842607235828485645766393 M929 = 2^929 - 1 = 4538015467766671944574165338592225830478699345884382504442663144885072806275648112625635725391102144133907238129251016389326737199538896813326509341743147661691195191795226666084858428449394948944821764472508048114220424520501343042471615418544488778723282182172070046459244838911 (0 msec) (M has 280 dec) = 13007 × 348890248924938259750454781163390930305120269538278042934009621348894657205785201247454118966026150852149399410259938217062100192168747352450719561908445272675574320888385228421992652298715687625495638077382028762529439880103124705348782610789919949159935587158612289264184273 M937 = 2^937 - 1 = 1161731959748268017810986326679609812602547032546401921137321765090578638406565916832162745700122148898280252961088260195667644723081957584211586391486245801392945969099578026517723757683045106929874371704962060317240428677248343818872733547147389127353160238636049931893566678761471 (0 msec) (M has 283 dec) = 28111 × 41326596696960905617409068573854000661753300577937530544531385048222355604801178073784737138491058621119143856891902109340387916583613446131819799775399160520541637405271175928203328152077304504637841830776637626453716647477796727931156257235508844486256634009321971181870679761
Scheme
This works with PLT Scheme, other implementations only need to change the inclusion.
#lang scheme
;;; this needs to be changed for other R6RS implementations
(require rnrs/arithmetic/bitwise-6)
;;; modpow, as per the task description.
(define (modpow exponent base)
(let loop ([square 1] [index (- (bitwise-length exponent) 1)])
(if (< index 0)
square
(loop (modulo (* (if (bitwise-bit-set? exponent index) 2 1)
square square) base)
(- index 1)))))
;;; search through all integers from 1 on to find the first divisor
;;; returns #f if 2^p-1 is prime
(define (mersenne-factor p)
(for/first ((i (in-range 1 (floor (expt 2 (quotient p 2))) (* 2 p)))
#:when (and (or (= 1 (modulo i 8)) (= 7 (modulo i 8)))
(= 1 (modpow p i))))
i))
- Output:
> (mersenne-factor 929) 13007 > (mersenne-factor 23) 47 > (mersenne-factor 3) #f
Seed7
$ include "seed7_05.s7i";
const func boolean: isPrime (in integer: number) is func
result
var boolean: prime is FALSE;
local
var integer: upTo is 0;
var integer: testNum is 3;
begin
if number = 2 then
prime := TRUE;
elsif odd(number) and number > 2 then
upTo := sqrt(number);
while number rem testNum <> 0 and testNum <= upTo do
testNum +:= 2;
end while;
prime := testNum > upTo;
end if;
end func;
const func integer: modPow (in var integer: base,
in var integer: exponent, in integer: modulus) is func
result
var integer: power is 1;
begin
if exponent < 0 or modulus < 0 then
raise RANGE_ERROR;
else
while exponent > 0 do
if odd(exponent) then
power := (power * base) mod modulus;
end if;
exponent >>:= 1;
base := base ** 2 mod modulus;
end while;
end if;
end func;
const func integer: mersenneFactor (in integer: exponent) is func
result
var integer: factor is 0;
local
var integer: maxk is 0;
var integer: k is 1;
var boolean: searching is TRUE;
begin
maxk := 16384 div exponent; # Limit for k to prevent overflow of 32 bit signed integer
while k <= maxk and searching do
factor := 2 * exponent * k + 1;
if (factor mod 8 = 1 or factor mod 8 = 7) and
isPrime(factor) and modPow(2, exponent, factor) = 1 then
searching := FALSE;
end if;
incr(k);
end while;
if searching then
factor := 0;
end if;
end func;
const proc: main is func
begin
writeln("Factor of M929: " <& mersenneFactor(929));
end func;
Original source: isPrime, modPow (modified to use integer instead of bigInteger).
- Output:
Factor of M929: 13007
Sidef
func mtest(b, p) {
var bits = b.base(2).digits
for (var sq = 1; bits; sq %= p) {
sq *= sq
sq += sq if bits.shift==1
}
sq == 1
}
for m (2..60 -> grep{ .is_prime }, 929) {
var f = 0
var x = (2**m - 1)
var q
{ |k|
q = (2*k*m + 1)
q%8 ~~ [1,7] || q.is_prime || next
q*q > x || (f = mtest(m, q)) && break
} << 1..Inf
say (f ? "M#{m} is composite with factor #{q}"
: "M#{m} is prime")
}
- Output:
M2 is prime M3 is prime M5 is prime M7 is prime M11 is composite with factor 23 M13 is prime M17 is prime M19 is prime M23 is composite with factor 47 M29 is composite with factor 233 M31 is prime M37 is composite with factor 223 M41 is composite with factor 13367 M43 is composite with factor 431 M47 is composite with factor 2351 M53 is composite with factor 6361 M59 is composite with factor 179951 M929 is composite with factor 13007
Swift
import Foundation
extension BinaryInteger {
var isPrime: Bool {
if self == 0 || self == 1 {
return false
} else if self == 2 {
return true
}
let max = Self(ceil((Double(self).squareRoot())))
for i in stride(from: 2, through: max, by: 1) where self % i == 0 {
return false
}
return true
}
func modPow(exp: Self, mod: Self) -> Self {
guard exp != 0 else {
return 1
}
var res = Self(1)
var base = self % mod
var exp = exp
while true {
if exp & 1 == 1 {
res *= base
res %= mod
}
if exp == 1 {
return res
}
exp >>= 1
base *= base
base %= mod
}
}
}
func mFactor(exp: Int) -> Int? {
for k in 0..<16384 {
let q = 2*exp*k + 1
if !q.isPrime {
continue
} else if q % 8 != 1 && q % 8 != 7 {
continue
} else if 2.modPow(exp: exp, mod: q) == 1 {
return q
}
}
return nil
}
print(mFactor(exp: 929)!)
- Output:
13007
Tcl
For primes::is_prime
see Prime decomposition#Tcl
proc int2bits {n} {
binary scan [binary format I1 $n] B* binstring
return [split [string trimleft $binstring 0] ""]
# another method
if {$n == 0} {return 0}
set bits [list]
while {$n > 0} {
lappend bits [expr {$n % 2}]
set n [expr {$n / 2}]
}
return [lreverse $bits]
}
proc trial_factor {base exp mod} {
set square 1
foreach bit [int2bits $exp] {
set square [expr {($square ** 2) * ($bit == 1 ? $base : 1) % $mod}]
}
return [expr {$square == 1}]
}
proc m_factor p {
set limit [expr {sqrt(2**$p - 1)}]
for {set k 1} {2 * $k * $p - 1 < $limit} {incr k} {
set q [expr {2 * $k * $p + 1}]
if { ! [primes::is_prime $q]} {
continue
} elseif { ! ($q % 8 == 1 || $q % 8 == 7)} {
# optimization
continue
} elseif {[trial_factor 2 $p $q]} {
# $q is a factor of 2**$p-1
return $q
}
}
return -1
}
set exp 929
if {[set fact [m_factor 929]] > 0} {
puts "M$exp has a factor: $fact"
} else {
puts "no factor found for M$exp"
}
TI-83 BASIC
The program uses the new remainder function from OS 2.53MP, if not available it can be replaced by:
remainder(A,B) equivalent to iPart(B*fPart(A/B))
Due to several problems, no Goto has been used. As a matter of fact the version is clearer.
Prompt Q
1→K:0→T
While K≤2^20 and T=0
2KQ+1→P
remainder(P,8)→W
If W=1 or W=7
Then
0→E:0→M
If remainder(P,2)=0:1→M
If remainder(P,3)=0:1→M
5→D
While M=0 and DD≤P
If remainder(P,D)=0:1→M
D+2→D
If remainder(P,D)=0:1→M
D+4→D
End
If M=0:1→E
Q→I:1→Y:2→Z
While I≠0
If remainder(I,2)=1:remainder(YZ,P)→Y
remainder(ZZ,P)→Z
iPart(I/2)→I
End
If E=1 and Y=1
Then
P→F:1→T
End
End
K+1→K
End
If T=0:0→F
Disp Q,F
- Input:
Q=?929
- Output:
929 13007 Done
uBasic/4tH
Print "A factor of M929 is "; FUNC(_FNmersenne_factor(929))
Print "A factor of M937 is "; FUNC(_FNmersenne_factor(937))
End
_FNmersenne_factor Param(1)
Local(2)
If (FUNC(_FNisprime(a@)) = 0) Then Return (-1)
For b@ = 1 TO 99999
c@ = (2*a@*b@) + 1
If (FUNC(_FNisprime(c@))) Then
If (AND (c@, 7) = 1) + (AND (c@, 7) = 7) Then
Until FUNC(_ModPow2 (a@, c@)) = 1
EndIf
EndIf
Next
Return (c@ * (b@<100000))
_FNisprime Param(1)
Local (1)
If ((a@ % 2) = 0) Then Return (a@ = 2)
If ((a@ % 3) = 0) Then Return (a@ = 3)
b@ = 5
Do Until ((b@ * b@) > a@) + ((a@ % b@) = 0)
b@ = b@ + 2
Until (a@ % b@) = 0
b@ = b@ + 4
Loop
Return ((b@ * b@) > a@)
_ModPow2 Param(2)
Local(2)
d@ = 1
For c@ = 30 To 0 Step -1
If ((a@+1) > SHL(1,c@)) Then
d@ = d@ * d@
If AND (a@, SHL(1,c@)) Then
d@ = d@ * 2
EndIf
d@ = d@ % b@
EndIf
Next
Return (d@)
- Output:
A factor of M929 is 13007 A factor of M937 is 28111 0 OK, 0:123
VBScript
' Factors of a Mersenne number
for i=1 to 59
z=i
if z=59 then z=929 ':) 61 turns into 929.
if isPrime(z) then
r=testM(z)
zz=left("M" & z & space(4),4)
if r=0 then
Wscript.echo zz & " prime."
else
Wscript.echo zz & " not prime, a factor: " & r
end if
end if
next
function modPow(base,n,div)
dim i,y,z
i = n : y = 1 : z = base
do while i
if i and 1 then y = (y * z) mod div
z = (z * z) mod div
i = i \ 2
loop
modPow= y
end function
function isPrime(x)
dim i
if x=2 or x=3 or _
x=5 or x=7 _
then isPrime=1: exit function
if x<11 then isPrime=0: exit function
if x mod 2=0 then isPrime=0: exit function
if x mod 3=0 then isPrime=0: exit function
i=5
do
if (x mod i) =0 or _
(x mod (i+2)) =0 _
then isPrime=0: exit function
if i*i>x then isPrime=1: exit function
i=i+6
loop
end function
function testM(x)
dim sqroot,k,q
sqroot=Sqr(2^x)
k=1
do
q=2*k*x+1
if q>sqroot then exit do
if (q and 7)=1 or (q and 7)=7 then
if isPrime(q) then
if modPow(2,x,q)=1 then
testM=q
exit function
end if
end if
end if
k=k+1
loop
testM=0
end function
- Output:
M2 prime. M3 prime. M5 prime. M7 prime. M11 not prime, a factor: 23 M13 prime. M17 prime. M19 prime. M23 not prime, a factor: 47 M29 not prime, a factor: 233 M31 prime. M37 not prime, a factor: 223 M41 not prime, a factor: 13367 M43 not prime, a factor: 431 M47 not prime, a factor: 2351 M53 not prime, a factor: 6361 M929 not prime, a factor: 13007
Visual Basic
Sub mersenne()
Dim q As Long, k As Long, p As Long, d As Long
Dim factor As Long, i As Long, y As Long, z As Long
Dim prime As Boolean
q = 929 'input value
For k = 1 To 1048576 '2**20
p = 2 * k * q + 1
If (p And 7) = 1 Or (p And 7) = 7 Then 'p=*001 or p=*111
'p is prime?
prime = False
If p Mod 2 = 0 Then GoTo notprime
If p Mod 3 = 0 Then GoTo notprime
d = 5
Do While d * d <= p
If p Mod d = 0 Then GoTo notprime
d = d + 2
If p Mod d = 0 Then GoTo notprime
d = d + 4
Loop
prime = True
notprime: 'modpow
i = q: y = 1: z = 2
Do While i 'i <> 0
On Error GoTo okfactor
If i And 1 Then y = (y * z) Mod p 'test first bit
z = (z * z) Mod p
On Error GoTo 0
i = i \ 2
Loop
If prime And y = 1 Then factor = p: GoTo okfactor
End If
Next k
factor = 0
okfactor:
Debug.Print "M" & q, "factor=" & factor
End Sub
- Output:
M47 factor=2351
V (Vlang)
import math
const qlimit = int(2e8)
fn main() {
mtest(31)
mtest(67)
mtest(929)
}
fn mtest(m int) {
// the function finds odd prime factors by
// searching no farther than sqrt(N), where N = 2^m-1.
// the first odd prime is 3, 3^2 = 9, so M3 = 7 is still too small.
// M4 = 15 is first number for which test is meaningful.
if m < 4 {
println("$m < 4. M$m not tested.")
return
}
flimit := math.sqrt(math.pow(2, f64(m)) - 1)
mut qlast := 0
if flimit < qlimit {
qlast = int(flimit)
} else {
qlast = qlimit
}
mut composite := []bool{len: qlast+1}
sq := int(math.sqrt(f64(qlast)))
loop:
for q := int(3); ; {
if q <= sq {
for i := q * q; i <= qlast; i += q {
composite[i] = true
}
}
q8 := q % 8
if (q8 == 1 || q8 == 7) && mod_pow(2, m, q) == 1 {
println("M$m has factor $q")
return
}
for {
q += 2
if q > qlast {
break loop
}
if !composite[q] {
break
}
}
}
println("No factors of M$m found.")
}
// base b to power p, mod m
fn mod_pow(b int, p int, m int) int {
mut pow := i64(1)
b64 := i64(b)
m64 := i64(m)
mut bit := u32(30)
for 1<<bit&p == 0 {
bit--
}
for {
pow *= pow
if 1<<bit&p != 0 {
pow *= b64
}
pow %= m64
if bit == 0 {
break
}
bit--
}
return int(pow)
}
- Output:
No factors of M31 found. M67 has factor 193707721 M929 has factor 13007
Wren
import "./math" for Int
import "./fmt" for Conv, Fmt
var trialFactor = Fn.new { |base, exp, mod|
var square = 1
var bits = Conv.itoa(exp, 2).toList
var ln = bits.count
for (i in 0...ln) {
square = square * square * (bits[i] == "1" ? base : 1) % mod
}
return square == 1
}
var mersenneFactor = Fn.new { |p|
var limit = (2.pow(p) - 1).sqrt.floor
var k = 1
while ((2*k*p - 1) < limit) {
var q = 2*k*p + 1
if (Int.isPrime(q) && (q%8 == 1 || q%8 == 7) && trialFactor.call(2, p, q)) {
return q // q is a factor of 2^p - 1
}
k = k + 1
}
return null
}
var m = [3, 5, 11, 17, 23, 29, 31, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 97, 929]
for (p in m) {
var f = mersenneFactor.call(p)
Fmt.write("2^$3d - 1 is ", p)
if (f) {
Fmt.print("composite (factor $d)", f)
} else {
System.print("prime")
}
}
- Output:
2^ 3 - 1 is prime 2^ 5 - 1 is prime 2^ 11 - 1 is composite (factor 23) 2^ 17 - 1 is prime 2^ 23 - 1 is composite (factor 47) 2^ 29 - 1 is composite (factor 233) 2^ 31 - 1 is prime 2^ 37 - 1 is composite (factor 223) 2^ 41 - 1 is composite (factor 13367) 2^ 43 - 1 is composite (factor 431) 2^ 47 - 1 is composite (factor 2351) 2^ 53 - 1 is composite (factor 6361) 2^ 59 - 1 is composite (factor 179951) 2^ 67 - 1 is composite (factor 193707721) 2^ 71 - 1 is composite (factor 228479) 2^ 73 - 1 is composite (factor 439) 2^ 79 - 1 is composite (factor 2687) 2^ 83 - 1 is composite (factor 167) 2^ 97 - 1 is composite (factor 11447) 2^929 - 1 is composite (factor 13007)
zkl
var [const] BN=Import("zklBigNum"); // libGMP
// M = 2^P - 1 , P prime
// Look for a prime divisor q such as:
// q < M.sqrt(), q = 1 or 7 modulo 8, q = 1 + 2kP
// q is divisor if 2.powmod(P,q) == 1
// m-divisor returns q or False
fcn m_divisor(P){
// must limit the search as M.sqrt() may be HUGE and I'm slow
maxPrime:='wrap{ BN(2).pow(P).sqrt().min(0d5_000_000) };
t,b2:=BN(0),BN(2); // so I can do some in place BigInt math
foreach q in (maxPrime(P*2)){ // 0..maxPrime -1, faster than just odd #s
if((q%8==1 or q%8==7) and t.set(q).probablyPrime() and
b2.powm(P,q)==1) return(q);
}
False
}
m_divisor(929).println(); // 13007
m_divisor(4423).println(); // False
(BN(2).pow(4423) - 1).probablyPrime().println(); // True
- Output:
13007 False True
- Programming Tasks
- Prime Numbers
- Arithmetic
- Arithmetic operations
- GUISS/Omit
- 11l
- 8086 Assembly
- 360 Assembly
- Ada
- ALGOL 68
- Arturo
- AutoHotkey
- BBC BASIC
- Bracmat
- C
- C sharp
- C++
- Clojure
- CoffeeScript
- Common Lisp
- Crystal
- D
- Delphi
- EasyLang
- EchoLisp
- Elixir
- Erlang
- Factor
- Forth
- Fortran
- FreeBASIC
- Frink
- GAP
- Go
- Haskell
- Icon
- Unicon
- J
- Java
- JavaScript
- Julia
- Kotlin
- Lingo
- Mathematica
- Wolfram Language
- Maxima
- Nim
- Octave
- PARI/GP
- Pascal
- Perl
- Phix
- PHP
- PicoLisp
- Prolog
- Python
- Racket
- Raku
- REXX
- Ring
- RPL
- Ruby
- Rust
- Scala
- Scheme
- Seed7
- Sidef
- Swift
- Tcl
- TI-83 BASIC
- UBasic/4tH
- VBScript
- Visual Basic
- V (Vlang)
- Wren
- Wren-math
- Wren-fmt
- Zkl