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:
- Compute c = n1/4 and define an array f containing c integers, each initially 1.
- In each array element fi, compute the product modulo n of the integers from i c + 1 to c(i+1).
- 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.
Programming Praxis 
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)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)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:
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
, it did not work for another problematic case,
. With the modification, the input value is returned as the factor.
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)