Strassen’s Factoring Algorithm

January 27, 2018

We’re going to cheat. The only reason to compute the fi is to take advantage of the fast product trees, and we’re not going to implement that, so instead of pre-computing the fi we compute them one at a time and immediately check the greatest common divisor, returning as soon as we find a factor:

(define (strassen n)
  (call-with-current-continuation
    (lambda (return)
      (let ((c (iroot 4 n)))
        (do ((i 0 (+ i 1))) ((= i c) #f)
          (let* ((jmin (+ 1 (* i c))) (jmax (+ jmin c -1)))
            (do ((j jmin (+ j 1)) (f 1 (modulo (* f j) n)))
                ((< jmax j) (let ((g (gcd f n)))
                              (when (< 1 g) (return g)))))))))))

As a practical matter, it is probably better to remove small factors with trial division before calling strassen, reducing the likelihood of finding a composite factor. Here are some examples:

> (strassen 1329059)
3119
> (strassen 11111111111111111)
2071723
> (strassen 2071723)
#f

Strassen’s factoring algorithm has no practical utility, as it is slower than other algorithms (Pollard rho, SQUFOF) with the same time complexity. It’s value is that it serves as an proven upper bound on the complexity of deterministically factoring integers.

You can run the program at https://ideone.com/HPO1aI.

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 )

Google+ photo

You are commenting using your Google+ 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 )

w

Connecting to %s

%d bloggers like this: