Strassen’s Factoring Algorithm

January 27, 2018

We’re going to cheat. The only reason to compute the f_i is to take advantage of the fast product trees, and we’re not going to implement that, so instead of pre-computing the f_i 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.

Posted by programmingpraxis

Filed in Exercises

5 Comments »

5 Responses to “Strassen’s Factoring Algorithm”

chaw said

January 27, 2018 at 6:00 PM

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)

Paul said

January 29, 2018 at 2:54 PM

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)

Daniel said

February 12, 2018 at 12:00 AM

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

Daniel said
February 12, 2018 at 12:02 AM
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.

Globules said

February 19, 2018 at 12:00 AM

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

Programming Praxis