Strassen’s Factoring Algorithm

January 27, 2018

In 1976, Volker Strassen, working with John Pollard, developed a factoring algorithm that is still the fastest proven deterministic factoring algorithm, with time complexity O(n1/4 log n); unfortunately, it is generally slower than other similar algorithms, even Pollard’s rho algorithm. The basic idea is that if you have the product fi of a consecutive set of integers modulo the number to be factored n, and that set of integers contains one of the factors of n, then gcd(fi, n) will be greater than 1. Here’s the algorithm:

  1. Compute c = n1/4 and define an array f containing c integers, each initially 1.
  2. In each array element fi, compute the product modulo n of the integers from i c + 1 to c(i+1).
  3. Compute the greatest common divisor of each of the fi with n, stopping when you find a (possibly-composite) factor.

The second loop, in Step 3, takes time O(n1/4 log n), so the trick is to compute the products of Step 2 in less time than O(n²). That can be done using product trees, as Berstein suggests, but that’s messy and a little bit complicated, and even if you manage to implement it, Strassen’s algorithm still isn’t competitive, so we won’t bother.

Your task is to implement Strassen’s factoring algorithm as described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Advertisements

Pages: 1 2

5 Responses to “Strassen’s Factoring Algorithm”

  1. chaw said

    Here is a naive implementation in very simple, standard Scheme,
    favoring clarity over efficiency.

    (import (scheme base)
            (scheme write)
            (only (srfi 1) list-tabulate))
    
    (define (factor/strassen n)
      (let g ((rp (range-products n)))
        (if (null? rp)
            n
            (let ((m (gcd n (car rp))))
              (if (> m 1)
                  m
                  (g (cdr rp)))))))
    
    (define (range-products n)
      (let ((c (exact (floor (expt n 1/4)))))
        (list-tabulate c (lambda (i)
                           (range-product-mod (+ (* i c) 1)
                                              (* c (+ i 1))
                                              n)))))
    
    (define (range-product-mod a b n)
      (let g ((a a)
              (p 1))
        (if (> a b)
            p
            (g (+ a 1)
               (modulo (* p a) n)))))
    
    (display (map factor/strassen '(23 24 101 17391743 302930532221)))
    (newline)
    

  2. Paul said

    In Python. I increased the size of the array by 1, as I did not get a factor for n= 97*101 otherwise.

    from functools import reduce
    from math import  gcd
    from ma.mymath import iroot  # integer root
    
    def prod(i, c, n):
        return reduce(lambda acc, x: acc * x % n,
                      range(i * c + 1, i * c + c + 1), 1)
    
    def strassen_factor(n):
        c = iroot(4, n)
        a = (gcd(prod(i, c, n), n) for i in range(c+2))
        return next((g for g in a if g > 1), -1)
    
  3. Daniel said

    Here’s a solution in C using GMP.

    @Paul, my implementation has the same problem you identified. While your modification worked for 97101, it did not work for another problematic case, 5961. With the modification, the input value is returned as the factor.

    /* strassen.c */
    
    // Build
    //   $ gcc -o strassen strassen.c -lgmp
    
    #include <gmp.h>
    
    void calc_strassen_factor(mpz_t factor_out, const mpz_t n) {
      mpz_t c;
      mpz_t i;
      mpz_t j;
      mpz_t next;
      mpz_t prod_mod;
      mpz_t gcd;
    
      mpz_init(c);
      mpz_init(i);
      mpz_init(j);
      mpz_init(next);
      mpz_init(prod_mod);
      mpz_init(gcd);
    
      mpz_root(c, n, 4);
      mpz_set_ui(i, 0);
      mpz_set_ui(next, 1);
      mpz_set_ui(factor_out, 1);
    
      while (mpz_cmp(c, i) >= 0) {
        // Calculates step 2 inefficiently (no product trees)
        mpz_set_ui(prod_mod, 1);
        mpz_set_ui(j, 0);
        while (mpz_cmp(j, c) < 0) {
          mpz_mul(prod_mod, prod_mod, next);
          mpz_mod(prod_mod, prod_mod, n);
          mpz_add_ui(next, next, 1);
          mpz_add_ui(j, j, 1);
        }
        mpz_gcd(gcd, prod_mod, n);
        if (mpz_cmp_ui(gcd, 1)) {
          mpz_set(factor_out, gcd);
          break;
        }
        mpz_add_ui(i, i, 1);
      }
    
      mpz_clear(c);
      mpz_clear(i);
      mpz_clear(j);
      mpz_clear(next);
      mpz_clear(prod_mod);
      mpz_clear(gcd);
    }
    
    int main(void) {
      mpz_t number;
      mpz_t factor;
      mpz_init(number);
      mpz_init(factor);
    
      long numbers[] = {13290059, 11111111111111111, 2071723};
      size_t n = sizeof(numbers) / sizeof(long);
    
      for (size_t i = 0; i < n; ++i) {
        mpz_set_ui(number, numbers[i]);
        calc_strassen_factor(factor, number);
        gmp_printf("strassen_factor(%Zd) = %Zd\n", number, factor);
      }
    
      mpz_clear(number);
      mpz_clear(factor);
    
      return 0;
    }
    

    Output:

    strassen_factor(13290059) = 3119
    strassen_factor(11111111111111111) = 2071723
    strassen_factor(2071723) = 1
    
  4. Daniel said

    My asterisks were used as markup for italics. Let me try again.

    @Paul, my implementation has the same problem you identified. While your modification worked for

    97*101

    , it did not work for another problematic case,

    59*61

    . With the modification, the input value is returned as the factor.

  5. Globules said

    A Haskell version.

    import Data.List (foldl')
    import Math.NumberTheory.Powers (integerFourthRoot)
    
    -- Compute (a0 * a1 * ... * an) mod m.  Take the remainder after each
    -- multiplication to keep the intermediate values small.
    prodRem :: Integral a => a -> [a] -> a
    prodRem m = foldl' (\a x -> (a * x) `rem` m) 1
    
    -- Return a factor of n, if one exists.  The result is only valid for n ≥ 16.
    factor :: Integral a => a -> Maybe a
    factor n = let c  = integerFourthRoot n
                   fs = [prodRem n [c*i+1 .. c*(i+1)] | i <- [0 .. c-1]]
               in case dropWhile (== 1) $ map (gcd n) fs of
                    []    -> Nothing
                    (i:_) -> Just i
    
    type I = Integer
    
    main :: IO ()
    main = do
      print $ factor (13290059 :: I)
      print $ factor (11111111111111111 :: I)
      print $ factor (2071723 :: I)
    
    $ ./strassfact 
    Just 3119
    Just 2071723
    Nothing
    

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

%d bloggers like this: