/* Package: bifactor Description: Naive bivariate factorization, neither efficient nor foolproof Authors: K.B. (1998-05-04, initial implementation) author2 (date2, what...) author3 (date3, what...) Maintainer: Karim Belabas \\ empty = orphaned Requires: \\ e.g. elldata, other script, optional if empty Files: bifactor.gp \\ optional if only 1 file Function: bifactor Help: bifactor(T) returns the list of irreducible factors of the bivariate polynomial T in K[x,y]. Doc: bifactor(T) returns the list of irreducible factors of the bivariate polynomial T in K[x,y], using Kronecker substitution to reduce to the univariate case. Infer the base field from the coefficients of T, as factor(). Assumption: let S be the squarefree part of T and n = deg_y S; there must exist a small integer i such that the univariate polynomial S((y+i)^(n+1), y) is squarefree. FIXME: uses global variables and local(), should use Hensel mod (x-i)^k rather than Kronecker substitution. Example: bifactor(x^2-y^2) %1 = [x-y, x+y] Example: bifactor(x^5+y*x^4+y^2*x^3+(y^2-1)*x^2+(y^3-y)*x+(y^4-y^2)) %2 = [x^2 + y*x + y^2, x^3 + (y^2 - 1)] Example: bifactor((x^2+y^2) * Mod(1,13)) %3 = [Mod(1, 13)*x + Mod(5, 13)*y, Mod(1, 13)*x + Mod(8, 13)*y] Example: bifactor(x^3+I*y*x^2+(y^2+I)*x+(I*y^3-y)) %4 = [x^2 + (y^2 + I), x + I*y] Function: function2 Doc: ... Example: textin textout */ \\ recover tentative factor from Kronecker; use global var BLOC split(p) = { my(v); Polrev( vector(1 + ceil(poldegree(p)/poldegree(BLOC)), k, v = divrem(p, BLOC); p = v[1]; v[2])); } \\ recursive recombination of local factors; use global var Lmod, Lfac \\ update global var GLOBALf cmbf(p, i) = { my(q, v); if (i > #Lmod, return); if (cmbf(p, i+1) && p != 1, return (1)); if (!Lmod[i], return); p *= Lmod[i]; q = split(p); if (poldegree(q, x), v = divrem(GLOBALf, q); if (!v[2], GLOBALf = v[1]; listput(Lfac, q); return (1)); ); if (cmbf(p, i+1), Lmod[i] = 0; return (1)); } /* Factor f(x, y) */ bifactor(f) = { local(Lfac, Lmod, BLOC); my (n,F,p,v); F = gcd(f, f'); if (poldegree(F), f /= F; \\ should do squarefree fact... f /= content(f); ); n = poldegree(f, 'y); for (i = 0, 10000, F = subst(f, 'x, BLOC = ('y + i)^(n+1)); if (issquarefree(F), break); ); if (!BLOC, error("bifactor fails")); Lmod = factor(F)[,1]; Lfac = listcreate(#Lmod); GLOBALf = f; cmbf(1, 1); if (poldegree(GLOBALf), listput(Lfac,GLOBALf)); Vec(Lfac); }