Factors of a Mersenne number: Difference between revisions
No edit summary |
Corrected - this trial factoring works for all Mersenne numbers, not just those over M_13. |
||
Line 12:
27*27 = 729 1 729*2 = 1458 1
Since 2<sup>23</sup> mod 47 = 1, 47 is a factor of 2<sup>P</sup>-1. (To see this, subtract 1 from both sides: 2<sup>23</sup>-1 = 0 mod 47.) Since we've shown that 47 is a factor, 2<sup>23</sup>-1 is not prime. Further properties of Mersenne numbers allow us to refine the process even more. Any factor q of 2<sup>P</sup>-1 must be of the form 2kp+1. Furthermore, q must be 1 or 7 mod 8. Finally any potential factor q must be [[Primality by Trial Division|prime]]. As in other trial division algorithms, the algorithm stops when 2kp+1 > sqrt(N).
This method only works for
'''Task:''' Using the above method find a factor of 2<sup>929</sup>-1 (aka M929)
|
Revision as of 13:34, 13 February 2009
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
PROC is prime = (INT number)BOOL:BEGIN SKIP # code omitted - see Primality by Trial Division # END; 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 := 1; FOR i FROM msb BY -1 TO 0 DO IF ( BIN p SHR i AND 2r1 ) = 2r1 THEN n := n*n*2 MOD q ELSE n := n*n MOD q FI OD; IF n = 1 THEN m factor := q; GO TO 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
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
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