Factors of a Mersenne number: Difference between revisions

Content added Content deleted

Inline

Revision as of 22:40, 21 August 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 - with prelude inserted manually

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

<lang algol>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

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

Usage:

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

: 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.

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>

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>