Factors of a Mersenne number
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:
You are encouraged to solve this task according to the task description, using any language you may know.
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
Template:In progress - last edit 9.III.2009
<lang mathematica>M = {929}; (* one or many numbers to check*) prime = Table[Prime[n], {n, 1000000}]; (* gives first 1000000 prime numbers *) For[i = 1, i < 1 + Length[M], i++,
Print["searching for factors on M", Mi, " in range <", First[prime], ", ", Last[prime], ">"]; For[l = 1, l < 1 + Length[prime], l++, S = 1; binary = IntegerDigits[Mi, 2]; For[n = 1, n < 1 + Length[binary], n++, S = Mod[(S*S) + (S*S*binaryn), primel]; ]; If[S == 1, Print[" ( ", primel, " is factor )"]; Break[];]; (* first factor BREAKS the search*) ]; Print["search enDEaD."];
];</lang> Example: <lang mathematica>searching for factors on M929 in range <2, 15485863>
( 13007 is factor )
search enDEaD.</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>
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
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
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>