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 (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 - with prelude inserted manually

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

<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

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

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

<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

Translation of: Fortran

(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>