Factors of a Mersenne number: Difference between revisions
| Line 548: | Line 548: | ||
</lang> |
</lang> |
||
Output: |
Output: |
||
<pre style="height: |
<pre style="height:42ex;overflow:scroll"> |
||
M2 ──────── is a Mersenne prime. |
M2 ──────── is a Mersenne prime. |
||
M3 ──────── is a Mersenne prime. |
M3 ──────── is a Mersenne prime. |
||
Revision as of 09:40, 13 November 2010
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 where P is 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. 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).
This method only works for Mersenne numbers where P is prime (M27 yields no factors).
Task: Using the above method find a factor of 2929-1 (aka M929)
ALGOL 68
<lang algol68>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</lang> Example:
Enter exponent of Mersenne number:929 M +929 has a factor: +13007
AutoHotkey
ahk discussion <lang autohotkey>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
}</lang>
Common Lisp
We can use a primality test from the Primality by Trial Division task.
<lang lisp>(defun primep (a)
(cond ((= a 2) T)
((or (<= a 1) (= (mod a 2) 0)) nil)
((loop for i from 3 to (sqrt a) by 2 do
(if (= (mod a i) 0)
(return nil))) nil)
(T T)))
(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)))))</lang>
Specific to this task, we define modulo-power and mersenne-prime-p.
<lang lisp>(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)))))</lang>
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.
Forth
<lang 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</lang>
Fortran
<lang 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</lang> Output
M929 has a factor: 13007
Haskell
Using David Amos module Primes [1] for prime number testing
<lang haskell>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</lang>
<lang haskell>*Main> trialfac 929 [13007]</lang>
J
<lang j>trialfac=: 3 : 0
qs=. (#~8&(1=|+.7=|))(#~1&p:)1+(*(1x+i.@<:@<.)&.-:)y qs#~1=qs&|@(2&^@[**:@])/ 1,~ |.#: y
)</lang> Examples: <lang j>trialfac 929 13007</lang> <lang j>trialfac 44497</lang>Empty output --> No factors found.
Mathematica
Believe it or not, this type of test runs faster in Mathematica than the squaring version described above.
<lang mathematica> 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"]</lang>
Octave
(GNU Octave has a isprime built-in test)
<lang octave>% 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));</lang>
PARI/GP
<lang parigp>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); }</lang>
PicoLisp
<lang 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)
(for (D 3 T (+ D 2))
(T (> D (sqrt N)) 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 ) ) ) )</lang>
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 "
(if (mFactor P)
(pack "composite with factor " @)
"prime" ) ) )
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
Python
<lang 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)</lang>
Example:
Enter exponent of Mersenne number: 929 M929 has a factor: 13007
REXX
REXX has no limit on numeric digits (precision). <lang rexx> /*REXX program uses exponent-&-mod operator to test possible Mersenne #s*/
numeric digits 500 /*we're dealing with some biggies*/
do j=1 by 2 to 61 /*when J=63, it turns into 929. */
z=j
if z== 1 then z= 2 /*oops, let's use 2 instead of 1.*/
if z==61 then z=929 /*switcheroo, 61 turns into 929.*/
if \isPrime(z) then iterate /*if not prime, keep pluging. */
r=testM(z) /*not, give it the 3rd degree. */
if r==0 then say right('M'z,5) "──────── is a Mersenne prime."
else say right('M'z,45) "is composite, a factor:" r
end
exit
testM:procedure; parse arg merc /*test a possible Mersenne Prime.*/
sqroot=iSqrt(2**merc) /*iSqrt = integer square root. */
do k=1 q=2*k*merc+1 if q>sqroot then leave /*if q>√(2^merc), then we're done*/ _=q//8 /*perform modulus arithmetic. */ if _\==1 & _\==7 then iterate /*must be either one or seven. */ if \isPrime(q) then iterate /*if not prime, keep on trukin'. */ if powmod(2,merc,q)==1 then return q /*Not prime? Return a factor.*/ end
return 0 /*it's a Mersenne prime, by gum. */
powmod: procedure; parse arg base,n,div
bits=strip(x2b(d2x(n)),'L',0) /*dec───>hex, then hex───>binary*/
sq=1
do forever until bits=
sq=sq**2
highbit=left(bits,1)
bits=substr(bits,2)
if highbit then do
sq=sq*base
sq=sq//div
end
end
return sq
/*─────────────────────────────────────ISQRT subroutine─────────────────*/ iSqrt: procedure; parse arg x; r=0; q=1; do while q<=x; q=q*4; end; do while q>1;q=q%4;_=x-r-q;r=r%2;if _>=0 then do;x=_;r=r+q;end;end;return r
/*─────────────────────────────────────ISPRIME subroutine───────────────*/ isPrime: procedure; arg x; if wordpos(x,'2 3 5 7')\==0 then return 1; if x<2 | 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 </lang> Output:
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: 13367
M43 is composite, a factor: 431
M47 is composite, a factor: 2351
M53 is composite, a factor: 6361
M59 is composite, a factor: 179951
M929 is composite, a factor: 13007
Ruby
<lang ruby>require 'mathn'
def mersenne_factor(p)
limit = Math.sqrt(2**p - 1)
k = 1
while (2*k*p - 1) < limit
q = 2*k*p + 1
if prime?(q) 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 prime?(value)
return false if value < 2 png = Prime.new for prime in png q,r = value.divmod prime return true if q < prime return false if r == 0 end
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
png = Prime.new for p in png
check_mersenne p break if p == 53
end p = 929 check_mersenne p</lang>
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
Scheme
This works with PLT Scheme, other implementations only need to change the inclusion.
<lang scheme>
- 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))
</lang>
> (mersenne-factor 929) 13007 > (mersenne-factor 23) 47 > (mersenne-factor 3) #f
Tcl
For primes::is_prime see Prime decomposition#Tcl
<lang 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"
}</lang>