polynomial_factorization_monte_carlo.sf

#!/usr/bin/ruby

# Factorize a polynomial using evaluation with random integer values and integer factorization.

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

# This method is not recommended, as it is quite slow.

func polynomial_factorization(f, tries = 500) {

    tries.times {
        var x = tries.irand
        var v = f(x).abs

        v.is_int || next

        v.divisors.each {|d|
            d > 1 || next
            d < v || next

            var e = d.ilog(x)
            e > 0 || break

            var (c,r) = divmod(d, x**e)
            var g = Poly([[0, r], [e, c]])

            if (f % g == 0) {
                return __FUNC__(g)+__FUNC__(f/g)
            }

            var (c,r) = divmod(x**(e+1), d)
            var g = Poly([[0, x-r], [e+1, c], [e, -1]])

            if (f % g == 0) {
                return __FUNC__(g)+__FUNC__(f/g)
            }
        }
    }

    return [f]
}

func display_factorization(f) {
    say (f.pretty, ' = ', polynomial_factorization(f).map{"(#{.pretty})"}.join(' * '))
}

display_factorization(Poly("9*x^4 - 1"))
display_factorization(Poly("2*x^4 - 6*x^3 + 7*x^2 - 5*x + 1"))
display_factorization(Poly("8*x^4 + 18*x^3 - 63*x^2 - 162*x - 81"))

__END__
9*x^4 - 1 = (3*x^2 + 1) * (3*x^2 - 1)
2*x^4 - 6*x^3 + 7*x^2 - 5*x + 1 = (x^2 - x + 1) * (2*x^2 - 4*x + 1)
8*x^4 + 18*x^3 - 63*x^2 - 162*x - 81 = (x + 3) * (2*x + 3) * (4*x^2 - 9*x - 9)