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