Factors of a Mersenne number

Factors of a Mersenne number
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 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. For example, let's see if 47 divides 2²³-1. 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 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

Once you get to the end of the binary number, if the modulo result is 1 then 47 is a factor of 2^P-1. Thus, 2²³ = 1 mod 47. 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.

This method only works for larger exponents (it reports M₁₃ is its own factor) and 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

Fortran

Works with: Fortran version 90 and later

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

<python>def is_prime(number):

   pass # 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 factor == 0:
           print "No factor found for M%d" % exponent
       else:
           print "M%d has a factor: %d" % (exponent, factor)</python>

Example:

Enter exponent of Mersenne number: 929
M929 has a factor: 13007