Factors of a Mersenne number: Difference between revisions

Content added Content deleted

Inline

Revision as of 14:22, 9 March 2009

A Mersenne number is a number in the form of 2^P-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 2^P-1 (or equivalently, if 2^P 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 2²³ 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 2²³ mod 47 = 1, 47 is a factor of 2^P-1. (To see this, subtract 1 from both sides: 2²³-1 = 0 mod 47.) Since we've shown that 47 is a factor, 2²³-1 is not prime. Further properties of Mersenne numbers allow us to refine the process even more. Any factor q of 2^P-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 (M₂₇ yields no factors).

Task: Using the above method find a factor of 2⁹²⁹-1 (aka M929)

ALGOL 68

Translation of: Fortran

Works with: ALGOL 68 version Standard - no extensions to language used

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

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

Works with: Fortran version 90 and later

<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

Mathematica

Template:In progress - last edit 9.III.2009

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", M[[i]], " in range <", First[prime], ", ", Last[prime], ">"];
  For[l = 1, l < 1 + Length[prime], l++,
    S = 1;
    binary = IntegerDigits[M[[i]], 2];
    For[n = 1, n < 1 + Length[binary], n++, 
      S = Mod[(S*S) + (S*S*binary[[n]]), prime[[l]]];
    ];
    If[S == 1, Print[" ( ", prime[[l]], " is factor )"];  Break[];]; (* first factor BREAKS the search*)
  ];
  Print["search enDEaD."];
];

Example:

searching for factors on M929 in range <2, 15485863>

 ( 13007 is factor )

search enDEaD.

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