Monte Carlo Factorization

June 19, 2009

Though the explanation is lengthy, the code is agreeably short:

(define (factor n . c)
  (define (f x c) (modulo (+ (* x x) c) n))
  (let ((c (if (pair? c) (car c) 1)))
    (let loop ((x 2) (y 2) (d 1))
      (cond ((= d 1)
              (let ((x (f x c)) (y (f (f y c) c)))
                (loop x y (gcd (- x y) n))))
            ((= d n) (factor n (+ c 1)))
            (else d)))))

Factor finds a single factor of n. The pseudo-random sequence is generated by f. The main loop starts with x = y = 2; it loops if the greatest common divisor d is 1, restarts with the next c if d is n, and otherwise reports d as a factor of n.

Pollard’s method works only if n is composite, and the factor function finds only a single factor. Factors, shown below, finds all the factors of n; factors is recursive, stopping when n is prime:

(define (factors n)
  (sort < (let fact ((n n) (fs '()))
    (cond ((= n 1) fs)
          ((even? n) (fact (/ n 2) (cons 2 fs)))
          ((prime? n) (cons n fs))
          (else (let ((f (factor n)))
                  (append fs (fact f '()) (fact (/ n f) '()))))))))

Factors uses the prime? function, and its companion check?, from the exercise on the Rabin-Miller primality checker. The factors of 298 – 1 = 316912650057057350374175801343 are 3, 43, 127, 4363953127297 and 4432676798593. You can run the program at http://programmingpraxis.codepad.org/4mQExWYY.

About these ads

Pages: 1 2

5 Responses to “Monte Carlo Factorization”

  1. [...] Praxis – Monte Carlo factorization By Remco Niemeijer In today’s Programming Praxis problem we have to implement John Pollard’s factorization algorithm. Our [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2009/06/19/programming-praxis-monte-carlo-factorization/ for a version with comments):

    import Control.Arrow  
    import Data.Bits  
    import Data.List  
    import System.Random 
    
    isPrime :: Integer -> StdGen -> Bool  
    isPrime n g =  
        let (s, d) = (length *** head) . span even $ iterate (`div` 2) (n-1)  
            xs = map (expm n d) . take 50 $ randomRs (2, n - 2) g  
        in all (\x -> elem x [1, n - 1] ||  
                      any (== n-1) (take s $ iterate (expm n 2) x)) xs  
      
    expm :: Integer -> Integer -> Integer -> Integer  
    expm m e b = foldl' (\r (b', _) -> mod (r * b') m) 1 .  
                 filter (flip testBit 0 . snd) .  
                 zip (iterate (flip mod m . (^ 2)) b) $  
                 takeWhile (> 0) $ iterate (`shiftR` 1) e
    
    factor :: Integer -> Integer -> Integer
    factor c n = factor' 2 2 1 where
        f x = mod (x * x + c) n
        factor' x y 1 = factor' x' y' (gcd (x' - y') n) where
                            (x', y') = (f x, f $ f y)
        factor' _ _ d = if d == n then factor (c + 1) n else d
    
    factors :: Integer -> StdGen -> [Integer]
    factors n g = sort $ fs n where
        fs x | even x      = 2 : fs (div x 2)
             | isPrime x g = [x]
             | otherwise   = f : fs (div x f) where f = factor 1 x
    
    main :: IO ()
    main = print . factors (2^98 - 1) =<< getStdGen
    
  3. Here’s my attempt in Python. A couple of issues in the code remain. The factors that it discovers aren’t guaranteed to be prime. I cribbed the Miller-Rabin test from one of the python code repositories. And, I don’t really understand exactly how this works. :-) Back to the reference books.

    #!/usr/bin/env python
    
    #
    # a basic implementation of the Pollard rho factorization
    # Written by Mark VandeWettering <markv@pixar.com>
    #
    
    import sys
    import locale
    import random
    
    class FactorError(Exception):
        def __init__(self, value):
            self.value = value
        def __str__(self):
            return repr(self.value)
    
    def miller_rabin_pass(a, n):
        d = n - 1
        s = 0
        while d % 2 == 0:
            d >>= 1
            s += 1
    
        a_to_power = pow(a, d, n)
        if a_to_power == 1:
            return True
        for i in xrange(s-1):
            if a_to_power == n - 1:
                return True
            a_to_power = (a_to_power * a_to_power) % n
        return a_to_power == n - 1
    
    def isprime(n):
        for repeat in xrange(20):
            a = 0
            while a == 0:
                a = random.randrange(n)
            if not miller_rabin_pass(a, n):
                return False
        return True
    
    def gcd(a, b):
        while b != 0:
            a, b = b, a%b
        return a
    
    def findfactor(n):
        for c in range(1, 50):
            x = y = random.randint(1, n-1)
            x = (x * x + c) % n
            y = (y * y + c) % n
            y = (y * y + c) % n
            while True:
                t = gcd(n, abs(x-y))
                if t == 1:
                    x = (x * x + c) % n
                    y = (y * y + c) % n
                    y = (y * y + c) % n
                elif t == n:
                    break
                else:
                    return t
        raise FactorError("couldn't find a factor.")
    
    def factor(n):
        r = []
        while True:
            if isprime(n):
                r.append(n)
                break
            try:
                f = findfactor(n)
                r.append(f)
                n = n / f
            except FactorError, msg:
                r.append(n)
                break
        r.sort()
        return r
    
    def doit(n):
            flist = factor(n)
            print locale.format("%d", n, True), "="
            for f in flist:
                    print "\t%s" % locale.format("%d", f, True)
    
    locale.setlocale(locale.LC_ALL, "")
    
    doit(2**98-1)
    
    
  4. Okay, I fixed a couple of things, and extended the program a tiny bit. It now is a numeric calculator of sorts. It’s not industrial strength or anything, but you can basically type any python numeric expression, and it will use eval() (at least with a predefined environment) to evaluate the number. I’ve also predefined a couple of built in functions. prime(n) will return an n digit prime. rsa(n) will return an rsa key which is the combination of two n/2 digit primes. factor(n) factors n. I’ve also added code to do some trial division as well, to get rid of small factors, and it collapses multiple occurrences of a factor (instead of printing 128 copies of 2 when factoring 2^128, it outputs “2**128″).

    #!/usr/bin/env python
    
    #
    # a basic implementation of the Pollard rho factorization
    # Written by Mark VandeWettering <markv@pixar.com>
    #
    
    import sys
    import locale
    import random
    try:
        import readline
    except ImportError, msg:
        print msg
        print "Line editing disabled."
    
    
    # an inefficient but straightforward way to find primes...
    
    def primes(n):
        primes = [2]
        for x in range(3, n, 2):
            prime = True
            for p in primes:
                if p * p > n:
                    break
                if x % p == 0:
                    # it's composite..
                    prime = False
                    break
            if prime:
                primes.append(x)
        return primes
    
    class FactorError(Exception):
        def __init__(self, value):
            self.value = value
        def __str__(self):
            return repr(self.value)
    
    def miller_rabin_pass(a, n):
        d = n - 1
        s = 0
        while d % 2 == 0:
            d >>= 1
            s += 1
    
        a_to_power = pow(a, d, n)
        if a_to_power == 1:
            return True
        for i in xrange(s-1):
            if a_to_power == n - 1:
                return True
            a_to_power = (a_to_power * a_to_power) % n
        return a_to_power == n - 1
    
    def isprime(n):
        for repeat in xrange(20):
            a = 0
            while a == 0:
                a = random.randrange(n)
            if not miller_rabin_pass(a, n):
                return False
        return True
    
    def gcd(a, b):
        while b != 0:
            a, b = b, a%b
        return a
    
    def findfactor(n):
        for c in range(1, 10):
            x = y = random.randint(1, n-1)
            x = (x * x + c) % n
            y = (y * y + c) % n
            y = (y * y + c) % n
            while True:
                t = gcd(n, abs(x-y))
                if t == 1:
                    x = (x * x + c) % n
                    y = (y * y + c) % n
                    y = (y * y + c) % n
                elif t == n:
                    break
                else:
                    return t
        raise FactorError("couldn't find a factor.")
    
    def factor(n):
        r = []
        for p in primes(10000):
            while n % p == 0:
                r.append(p)
                n = n / p
        if n == 1:
            return r
        while True:
            if isprime(n):
                r.append(n)
                break
            try:
                f = findfactor(n)
                r.append(f)
                n = n / f
            except FactorError, msg:
                r.append(n)
                break
        r.sort()
        return r
    
    # this function would be easier to write recursively, but
    # python isn't good at tail recursion, so in theory, it could
    # fail.  Too bad.
    
    def shorten(flist):
        slist = []
        idx = 0
        while flist[idx:] != []:
            hd = flist[idx]
            idx = idx + 1
            exp = 1
            while flist[idx:] != [] and flist[idx] == hd:
                exp = exp + 1
                idx = idx + 1
            if exp > 1:
                    slist.append(locale.format("%d", hd, True) + "**"+str(exp))
            else:
                    slist.append(locale.format("%d", hd, True))
        return slist
    
    def factorit(n):
            flist = factor(n)
            print locale.format("%d", n, True), "="
            for f in shorten(flist):
                    print "\t%s" % f
    
    locale.setlocale(locale.LC_ALL, "")
    
    import string
    
    def prime(n):
            "generate an n bit prime"
            while True:
                    x = int(''.join([random.choice(string.digits) for i in range(n)]))
                    if isprime(x):
                            return x
    
    def rsapair(n):
            return prime(n//2)*prime(n//2)
    
    env = { "prime" : prime,
            "rsa" : rsapair,
            "factor" : factorit}
    
    while True:
            try:
                    num = raw_input("enter a python numexpr >> ")
                    n = eval(num, env)
                    if n:
                            print n
            except NameError, msg:
                    print >> sys.stderr, msg
            except SyntaxError, msg:
                    print >> sys.stderr, msg
            except KeyboardInterrupt, msg:
                    print >> sys.stderr, "**Interrupted**"
                    continue
            except EOFError, msg:
                    print
                    break
    
  5. programmingpraxis said

    Instead of eval, you might want to build your own calculator. See the very first Programming Praxis exercise for an RPN calculator.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 621 other followers

%d bloggers like this: