polynomial_factorization_in_finite_field.sf

7#!/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_poly  = PolyMod(1 => 1, p)   # monomial x
    var w       = x_poly
    var factors = []

    for k in (1..f.degree) {
        w = w.powmod(p, f)              # w = x^(p^k) mod f
        var diff = (w - x_poly)
        var gk   = diff.monic_gcd(f)

        factors << [k, gk]
        f /= gk                   # exact division
        break if (f==1)
    }

    # Safety net for non-squarefree or degenerate input
    if (f != 1) {
        factors << [f.degree, f]
    }

    return factors
}

func distinct_degree_split(f, d) {

    var p = f.modulus
    var deg_f = f.degree

    return []   if (deg_f == 0)
    return [f]  if (deg_f == d)

    var is_char2 = (p == 2)
    var exp = (is_char2 ? nil : ((p**d - 1)/2))

    var g = nil
    var h = nil

    for (1..200) {

        # Random non-zero polynomial of degree < deg_f
        var t = PolyMod(deg_f.of { irand(1, p-1) }, p)
        next if (t == 0)

        if (is_char2) {

            # Absolute trace map
            h = t%f
            var ti = h

            for (1 .. d-1) {
                ti = ti.powmod(p, f)
                h  = ((h+ti) % f)
            }
        } else {
            # Quadratic-residue split
            h = (t.powmod(exp, f) - 1)
        }

        g = h.monic_gcd(f)
        break if ((g != 1) && (g != f))
    }

    return [f] if !(defined(g) && (g != 1) && (g != f))

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

func cantor_zassenhaus_polynomial_factorization(original) {

    # Stage 1 – preprocessing
    var lc = original.lift.leading_coefficient
    var f  = original.normalize_to_monic

    var df = f.derivative
    var sq = f.monic_gcd(df)

    f = f/sq
    f = f.normalize_to_monic if !(f==0 || f==1)

    # Stage 2 – distinct degree factorization
    var ddf_groups = (f==1 ? [] : distinct_degree_factorization(f))

    # Stage 3 – equal degree factorization
    var irreducibles = []

    ddf_groups.each_2d {|d, fk|
        next if (fk == 1)
        irreducibles << distinct_degree_split(fk, d)...
    }

    # Stage 4 – multiplicities in the original polynomial
    var remainder = original
    var factors = []

    for irred in (irreducibles) {
        var e = 0

        loop {
            var (q,r) = remainder.divmod(irred)
            break if (r != 0)
            e += 1
            remainder = q
        }

        next if (e == 0)
        factors << [irred, e]
    }

    [[lc, 1], factors.sort...]
}

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]:"
    var F = cantor_zassenhaus_polynomial_factorization(f)

    assert_eq(F.prod_2d{|a,b| a**b }, f)
    assert_eq(f.factor_exp, F)

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

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

Factoring: 3*x^4 + 3*x^3 + 2*x^2 + 3*x + 4 (mod 5) in Z_5[x]:
3
(x + 2 (mod 5))^2
x + 3 (mod 5)
x + 4 (mod 5)

Factoring: 7*x^5 + 9*x^4 + 11*x^2 + 9*x + 1 (mod 17) in Z_17[x]:
7
x + 8 (mod 17)
x^2 + x + 1 (mod 17)
x^2 + 2*x + 7 (mod 17)