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

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_gcd(gcd, prod_mod, n);
if (mpz_cmp_ui(gcd, 1)) {
mpz_set(factor_out, gcd);
break;
}
}

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

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