polynomial_factorization_in_finite_field.sf

#!/usr/bin/ruby

# Factorize a polynomial in Z_n[x], using the Cantor-Zassenhaus algorithm.

# Reference:
#   Lecture 15, Week 8, 1hr - Polynomial factorization
#   https://youtube.com/watch?v=tMqKqsMb-Ro

func distinct_degree_factorization(f) {

    var p = f.modulus
    var x = PolyMod(1, p)
    var w = x
    var g = []

    for k in (1..f.deg) {
        w = (w**p % f)
        g << gcd(w - x, f)
        f = (f / g[-1])
    }

    if (f != 1) {
        g << f
    }

    return g
}

func distinct_degree_split(f, d) {

    f.deg == 0 && return []
    f.deg == 1 && return [f]
    f.deg == d && return [f]

    var p = f.modulus
    var w = PolyMod(irand(d), p)       # random monic of degree <= d
    var n_power = idiv(p**d - 1, 2)
    var g = gcd(w**n_power - 1, f)

    if (g == f) {
        return __FUNC__(f, d)
    }

    if (g.deg == 0) {   # workaround for infinite recursion
        return [f/g]
    }

    __FUNC__(g, d) + __FUNC__(f/g, d)
}

func cantor_zassenhaus_polynomial_factorization(f) {

    var a = f

    f /= f.leading_coefficient      # make f monic
    f /= gcd(f, f.derivative)       # make f squarefree
    f /= f.content                  # make f primitive

    var ddf = distinct_degree_factorization(f)
    var dds = ddf.map_kv {|k,v| distinct_degree_split(v, k+1)... }

    var extra_factor = a/dds.prod

    if (extra_factor.deg > 0) {
        if (dds.none {|g| gcd(g, extra_factor) == g }) {
            dds << extra_factor
        }
    }

    var valuations = []

    var p = a
    for g in dds {
        var e = 0
        loop {
            var(q, r) = divmod(p, g)
            r == 0 || break
            ++e
            p = q
        }
        valuations << e
    }

    [[leading_coefficient(a), 1]] + (dds ~Z valuations)
}

for f in (
    PolyMod("8*x^4", 101),
    PolyMod("8*x^4 + 18*x^3 - 63*x^2 - 162*x - 81", 5),
    PolyMod("7*x^3 + 2*x^2 + 8*x + 1", 17)*PolyMod("x^2+x+1", 17)
) {

    say "\nFactoring: #{f} in Z_#{f.modulus}[x]:"

    cantor_zassenhaus_polynomial_factorization(f).each_2d{|factor, exp|
        if (exp > 1) {
            say "(#{factor})^#{exp}"
        }
        else {
            say factor
        }
    }
}

__END__
Factoring: 8*x^4 in Z_101[x]:
8
(x)^4

Factoring: 3*x^4 + 3*x^3 + 2*x^2 + 3*x + 4 in Z_5[x]:
3
x + 4
x^2 + 1
3*x + 1

Factoring: 7*x^5 + 9*x^4 + 11*x^2 + 9*x + 1 in Z_17[x]:
7
11*x + 3
4*x^2 + 4*x + 4
11*x^2 + 5*x + 9