diff options
Diffstat (limited to 'macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua')
-rw-r--r-- | macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua b/macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua new file mode 100644 index 0000000000..be19298e64 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua @@ -0,0 +1,220 @@ +-- Methods related to the Zassenhaus factorization algorithm. + + +-- Square-free factorization in the rational field. +function PolynomialRing:rationalsquarefreefactorization(keeplc) + local monic = self / self:lc() + local terms = {} + terms[0] = PolynomialRing.gcd(monic, monic:derivative()) + local b = monic // terms[0] + local c = monic:derivative() // terms[0] + local d = c - b:derivative() + local i = 1 + while b.degree ~= Integer.zero() or b.coefficients[0] ~= Integer.one() do + terms[i] = PolynomialRing.gcd(b, d) + b, c = b // terms[i], d // terms[i] + i = i + 1 + d = c - b:derivative() + end + if keeplc and terms[1] then + terms[1] = terms[1] * self:lc() + end + return terms +end + +-- Factors the largest possible constant out of a polynomial whos underlying ring is a Euclidean domain but not a field +function PolynomialRing:factorconstant() + local gcd = Integer.zero() + for i = 0, self.degree:asnumber() do + gcd = self.ring.gcd(gcd, self.coefficients[i]) + end + if gcd == Integer.zero() then + return Integer.one(), self + end + return gcd, self / gcd +end + +-- Converts a polynomial in the rational polynomial ring to the integer polynomial ring +function PolynomialRing:rationaltointeger() + local lcm = Integer.one() + for i = 0, self.degree:asnumber() do + if self.coefficients[i]:getring() == Rational:getring() then + lcm = lcm * self.coefficients[i].denominator / Integer.gcd(lcm, self.coefficients[i].denominator) + end + end + return Integer.one() / lcm, self * lcm +end + +-- Uses Zassenhaus's Algorithm to factor sqaure-free polynomials over the intergers +function PolynomialRing:zassenhausfactor() + + -- Creates a monic polynomial V with related roots + local V = {} + local n = self.degree:asnumber() + local l = self:lc() + for i = 0, n - 1 do + V[i] = l ^ Integer(n - 1 - i) * self.coefficients[i] + end + V[n] = Integer.one() + V = PolynomialRing(V, "y", self.degree) + + -- Performs Berlekamp Factorization in a sutable prime base + local p = V:findprime() + local S = V:inring(PolynomialRing.R("y", p)):berlekampfactor() + + -- If a polynomial is irreducible with coefficients in mod p, it is also irreducible over the integers + if #S == 1 then + return {self} + end + + -- Performs Hensel lifting on the factors mod p + local k = V:findmaxlifts(p) + local W = V:henselift(S, k) + local M = {} + + -- Returns the solutions back to the original from the monic transformation + for i, factor in ipairs(W) do + local w = {} + for j = 0, factor.degree:asnumber() do + w[j] = factor.coefficients[j]:inring(Integer.getring()) * l ^ Integer(j) + end + _, M[i] = PolynomialRing(w, self.symbol, factor.degree):factorconstant() + end + + return M + +end + +-- Finds the smallest prime such that this polynomial with coefficients in mod p is square-free +function PolynomialRing:findprime() + + local smallprimes = {Integer(2), Integer(3), Integer(5), Integer(7), Integer(11), Integer(13), Integer(17), Integer(19), Integer(23), + Integer(29), Integer(31), Integer(37), Integer(41), Integer(43), Integer(47), Integer(53), Integer(59)} + + for _, p in pairs(smallprimes) do + local P = PolynomialRing({IntegerModN(Integer.one(), p)}, self.symbol) + local s = self:inring(P:getring()) + if PolynomialRing.gcd(s, s:derivative()) == P then + return p + end + end + + error("Execution error: No suitable prime found for factoring.") +end + +-- Finds the maximum number of times Hensel Lifting will be applied to raise solutions to the appropriate power +function PolynomialRing:findmaxlifts(p) + local n = self.degree:asnumber() + local h = self.coefficients[0] + for i=0 , n do + if self.coefficients[i] > h then + h = self.coefficients[i] + end + end + + local B = 2^n * math.sqrt(n) * h:asnumber() + return Integer(math.ceil(math.log(2*B, p:asnumber()))) +end + +-- Uses Hensel lifting on the factors of a polynomial S mod p to find them in the integers +function PolynomialRing:henselift(S, k) + local p = S[1].ring.modulus + if k == Integer.one() then + return self:truefactors(S, p, k) + end + G = self:genextendsigma(S) + local V = S + for j = 2, k:asnumber() do + local Vp = V[1]:inring(PolynomialRing.R("y")) + for i = 2, #V do + Vp = Vp * V[i]:inring(PolynomialRing.R("y")) + end + local E = self - Vp:inring(PolynomialRing.R("y")) + if E == Integer.zero() then + return V + end + E = E:inring(PolynomialRing.R("y", p ^ Integer(j))):inring(PolynomialRing.R("y")) + F = E / p ^ (Integer(j) - Integer.one()) + R = self:genextendR(V, G, F) + local Vnew = {} + for i, v in ipairs(V) do + local vnew = v:inring(PolynomialRing.R("y", p ^ Integer(j))) + local rnew = R[i]:inring(PolynomialRing.R("y", p ^ Integer(j))) + Vnew[i] = vnew + (p) ^ (Integer(j) - Integer.one()) * rnew + end + V = Vnew + end + return self:truefactors(V, p, k) +end + +-- Gets a list of sigma polynomials for use in hensel lifting +function PolynomialRing:genextendsigma(S) + local v = S[1] * S[2] + local _, A, B = PolynomialRing.extendedgcd(S[2], S[1]) + local SIGMA = {A, B} + for i, _ in ipairs(S) do + if i >= 3 then + v = v * S[i] + local sum = SIGMA[1] * (v // S[1]) + for j = 2, i-1 do + sum = sum + SIGMA[j] * (v // S[j]) + end + _, A, B = PolynomialRing.extendedgcd(sum, v // S[i]) + for j = 1, i-1 do + SIGMA[j] = SIGMA[j] * A + end + SIGMA[i] = B + end + end + + return SIGMA +end + +-- Gets a list of r polynomials for use in hensel lifting +function PolynomialRing:genextendR(V, G, F) + R = {} + for i, v in ipairs(V) do + local pring = G[1]:getring() + R[i] = F:inring(pring) * G[i] % v:inring(pring) + end + return R +end + +-- Updates factors of the polynomial to the correct ones in the integer ring +function PolynomialRing:truefactors(l, p, k) + local U = self + local L = l + local factors = {} + local m = 1 + while m <= #L / 2 do + local C = Subarrays(L, m) + while #C > 0 do + local t = C[1] + local prod = t[1] + for i = 2, #t do + prod = prod * t[i] + end + local T = prod:inring(PolynomialRing.R("y", p ^ k)):inring(PolynomialRing.R("y")) + -- Convert to symmetric representation - this is the only place it actually matters + for i = 0, T.degree:asnumber() do + if T.coefficients[i] > p ^ k / Integer(2) then + T.coefficients[i] = T.coefficients[i] - p^k + end + end + local Q, R = U:divremainder(T) + if R == Integer.zero() then + factors[#factors+1] = T + U = Q + L = RemoveAll(L, t) + C = RemoveAny(C, t) + else + C = Remove(C, t) + end + end + m = m + 1 + end + if U ~= Integer.one() then + factors[#factors+1] = U + end + return factors +end
\ No newline at end of file |