In [1]:
R.<y,x> = PolynomialRing(QQ)

# There is a reason for the order of variables
# I'll later want to do f % y^i and for that sage
# needs y to have "higher precedence" than x.
#
# ... just ignore this for now.

R

Out[1]:
Multivariate Polynomial Ring in y, x over Rational Field
In [2]:
def extended_euclid(f , g, R = R):
old_d, d = f       , g
old_u, u = R.one() , R.zero()
old_v, v = R.zero(), R.one()

while not d.is_zero():
q         = old_d // d
old_d , d = d , (old_d - q * d)
old_u , u = u , (old_u - q * u)
old_v , v = v , (old_v - q * v)
l = old_d.lc()
return (old_d / l, old_u / l, old_v / l)

In [3]:
def HenselLift(f, g, h, a, b, y_mon, R = R, monic = False):
'''
Expecting polynomials f, g, h, a, b, y_mon
from the ring R such that
*  f - gh  = 1 mod y_mon
*  ag + bh = 1 mod y_mon
and y_mon is y^i for some i.

It will return polynomials g', h', a', b' such that
*  g = g'   mod y_mon^2
*  h = h'   mod y_mon^2
*  f = g'h' mod y_mon^2
*  a'g' + b' h' = 1 mod y_mon^2
'''

if (f - g*h) % y_mon != R.zero(): raise ValueError("f - gh is not 0 mod I")
if (a*g + b*h) % y_mon != R.one(): raise ValueError("ag + bh is not 1 mod I")

gamma = f - g * h
g1 = g + (b * gamma)
h1 = h + (a * gamma)

if monic:
q , r = (g1 - g) // g , (g1 - g) % g
g1 = (g + r) % y_mon^2
h1 = (h1 * ( 1 + q)) % y_mon^2

delta = a*g1 + b*h1 - R.one()
a1 = a*(1 - delta)
b1 = b*(1 - delta)
return (g1, h1, a1, b1)

def iteratedHenselLift(f, g, h, a, b, y_mon, k, R = R, monic = False):
if k == 1:
return HenselLift(f, g, h, a, b, y_mon, R, monic)
else:
g1, h1, a1, b1 = iteratedHenselLift(f, g, h, a, b, y_mon, k-1, R, monic)
y_mon_1 = y^(2**(k-1))
return HenselLift(f,g1, h1, a1, b1, y_mon_1, R, monic)


## Example 1¶

\begin{align*} f(x,y) & = x^2 - 2xy - 3y^2 + 3x - 5y + 2\\ & = (x + y + 2)(x - 3y + 1) \end{align*}
In [4]:
f = x^2 - 2*x*y - 3*y^2 + 3*x - 5*y + 2
f % y

Out[4]:
x^2 + 3*x + 2
In [5]:
factor(f % y)

Out[5]:
(x + 1) * (x + 2)
In [6]:
g = x + 1
h = x + 2
_ , a, b = extended_euclid(g, h)
a*g + b*h == 1

Out[6]:
True
In [7]:
gk, hk, ak, bk = iteratedHenselLift(f, g, h, a, b, y, k = 5, monic = True)

In [8]:
print(gk)
print(hk)

-3*y + x + 1
y + x + 2


Great, we found our factors!

## Example 2¶

Let's take $f(x,y) = x^3 + x - y$ which is actually irreducible.

However, $f \bmod{y} = x^3 + x = x (x^2 + 1)$.

In [9]:
f = x^3 + x - y
f0 = (f % y)
factor(f0)

Out[9]:
x * (x^2 + 1)
In [10]:
g = x
h = x^2 + 1
_, a,b = extended_euclid(g, h)

In [11]:
gk, hk, ak, bk = iteratedHenselLift(f, g, h, a, b, y, k = 5, monic = True)

In [12]:
print(gk)

11124755664*y^31 - 1822766520*y^29 + 300830572*y^27 - 50067108*y^25 + 8414640*y^23 - 1430715*y^21 + 246675*y^19 - 43263*y^17 + 7752*y^15 - 1428*y^13 + 273*y^11 - 55*y^9 + 12*y^7 - 3*y^5 + y^3 - y + x


We appear to be getting some garbage. Perhaps this suggests that we started off with an $f(x,y)$ that was irreducible.

(Spoiler: This intuition is wrong.)

## Example 3¶

$$f(x,y) = x^4 - y^2 - 2y - 1$$
In [13]:
f = x^4 - y^2 - 2*y - 1
f

Out[13]:
x^4 - y^2 - 2*y - 1
In [14]:
f0 = f % y
factor(f0)

Out[14]:
(x - 1) * (x + 1) * (x^2 + 1)
In [15]:
g = x - 1
h = (x + 1)*(x^2 + 1)
_, a,b = extended_euclid(g, h)

In [16]:
gk, hk, ak, bk = iteratedHenselLift(f, g, h, a, b, y, k = 5, monic = True)

In [17]:
print(f"gk: {gk}")
print("")
print(f"hk: {hk}")

gk: -238436656380769/144115188075855872*y^31 + 125280277081421/72057594037927936*y^30 - 32968493968795/18014398509481984*y^29 + 17383387729001/9007199254740992*y^28 - 2295919134019/1125899906842624*y^27 + 1215486600363/562949953421312*y^26 - 322476036831/140737488355328*y^25 + 171529806825/70368744177664*y^24 - 11435320455/4398046511104*y^23 + 6116566755/2199023255552*y^22 - 1641030105/549755813888*y^21 + 883631595/274877906944*y^20 - 119409675/34359738368*y^19 + 64822395/17179869184*y^18 - 17678835/4294967296*y^17 + 9694845/2147483648*y^16 - 334305/67108864*y^15 + 185725/33554432*y^14 - 52003/8388608*y^13 + 29393/4194304*y^12 - 4199/524288*y^11 + 2431/262144*y^10 - 715/65536*y^9 + 429/32768*y^8 - 33/2048*y^7 + 21/1024*y^6 - 7/256*y^5 + 5/128*y^4 - 1/16*y^3 + 1/8*y^2 - 1/2*y + x - 1

hk: 238436656380769/144115188075855872*y^31*x^2 - 125280277081421/72057594037927936*y^30*x^2 - 12123897782073/144115188075855872*y^31 + 32968493968795/18014398509481984*y^29*x^2 + 6593698793759/72057594037927936*y^30 - 17383387729001/9007199254740992*y^28*x^2 - 1798281489207/18014398509481984*y^29 + 2295919134019/1125899906842624*y^27*x^2 + 983965343151/9007199254740992*y^28 - 1215486600363/562949953421312*y^26*x^2 - 135054066707/1125899906842624*y^27 + 322476036831/140737488355328*y^25*x^2 + 74417546961/562949953421312*y^26 - 171529806825/70368744177664*y^24*x^2 - 20583576819/140737488355328*y^25 + 11435320455/4398046511104*y^23*x^2 + 11435320455/70368744177664*y^24 - 6116566755/2199023255552*y^22*x^2 - 797813055/4398046511104*y^23 + 1641030105/549755813888*y^21*x^2 + 447553665/2199023255552*y^22 - 883631595/274877906944*y^20*x^2 - 126233085/549755813888*y^21 + 119409675/34359738368*y^19*x^2 + 71645805/274877906944*y^20 - 64822395/17179869184*y^18*x^2 - 10235115/34359738368*y^19 + 17678835/4294967296*y^17*x^2 + 5892945/17179869184*y^18 - 9694845/2147483648*y^16*x^2 - 1710855/4294967296*y^17 + 334305/67108864*y^15*x^2 + 1002915/2147483648*y^16 - 185725/33554432*y^14*x^2 - 37145/67108864*y^15 + 52003/8388608*y^13*x^2 + 22287/33554432*y^14 - 29393/4194304*y^12*x^2 - 6783/8388608*y^13 + 4199/524288*y^11*x^2 + 4199/4194304*y^12 - 2431/262144*y^10*x^2 - 663/524288*y^11 + 715/65536*y^9*x^2 + 429/262144*y^10 - 429/32768*y^8*x^2 - 143/65536*y^9 + 33/2048*y^7*x^2 + 99/32768*y^8 - 21/1024*y^6*x^2 - 9/2048*y^7 + 7/256*y^5*x^2 + 7/1024*y^6 - 5/128*y^4*x^2 - 3/256*y^5 + 1/16*y^3*x^2 + 3/128*y^4 - 1/8*y^2*x^2 - 1/16*y^3 + 1/2*y*x^2 + x^3 + 3/8*y^2 + y*x + x^2 + 3/2*y + x + 1


Garbage again. Suggests that $f(x,y)$ is irreducible, right?

In [18]:
factor(f)

Out[18]:
(x^2 - y - 1) * (x^2 + y + 1)

Hunh?!?!

### Making sense of this¶

The polynomial $f(x,y)$ actually has two irreducible factors, $g^\ast = (x^2 - y - 1)$ and $h^\ast = (x^2 + y + 1)$. However, when we went modulo $y$, the polynomial $g^\ast = x^2 - 1 \bmod{y}$ factorised further into $(x - 1)$ and $(x + 1)$. We picked up $g = x - 1$ to start the entire lifting process, and $g$ was only part of an actual factor of $f(x,y)$. Perhaps it isn't surprising that the lifted version looked messy.

But how do we really distinguish the mess in Example 2 and Example 3?

Suppose on the other hand, we did the entire lifting process starting with $g_1 = x +1$ instead, what would we end up with?

In [19]:
g1 = (x + 1)
h1 = (x - 1) * (x^2 + 1)
_, a1, b1 = extended_euclid(g1, h1)
gk_prime, _, _, _ = iteratedHenselLift(f, g1, h1, a1, b1, y, k=5, monic = True)
gk_prime

Out[19]:
238436656380769/144115188075855872*y^31 - 125280277081421/72057594037927936*y^30 + 32968493968795/18014398509481984*y^29 - 17383387729001/9007199254740992*y^28 + 2295919134019/1125899906842624*y^27 - 1215486600363/562949953421312*y^26 + 322476036831/140737488355328*y^25 - 171529806825/70368744177664*y^24 + 11435320455/4398046511104*y^23 - 6116566755/2199023255552*y^22 + 1641030105/549755813888*y^21 - 883631595/274877906944*y^20 + 119409675/34359738368*y^19 - 64822395/17179869184*y^18 + 17678835/4294967296*y^17 - 9694845/2147483648*y^16 + 334305/67108864*y^15 - 185725/33554432*y^14 + 52003/8388608*y^13 - 29393/4194304*y^12 + 4199/524288*y^11 - 2431/262144*y^10 + 715/65536*y^9 - 429/32768*y^8 + 33/2048*y^7 - 21/1024*y^6 + 7/256*y^5 - 5/128*y^4 + 1/16*y^3 - 1/8*y^2 + 1/2*y + x + 1

That's messy too, but that's also expected for the same reasons mentioned above --- $x+1$ is only part of a real factor of $f(x,y)$ and so its lift looks messy.

However, $g_k'$ and $g_k$ together are part of a proper factor of $f$. What is their product going to be?

In [20]:
gk_prime * gk

Out[20]:
-56852039106040910281913031361/20769187434139310514121985316880384*y^62 + 29871410377750308792191592749/5192296858534827628530496329220096*y^61 - 47138737696913734391959792661/5192296858534827628530496329220096*y^60 + 1034392363255160315811029933/81129638414606681695789005144064*y^59 - 2727221176208987724407432445/162259276829213363391578010288128*y^58 + 1728004005766030877481684087/81129638414606681695789005144064*y^57 - 2131952994126921212477402445/81129638414606681695789005144064*y^56 + 5040078000300050147700715/158456325028528675187087900672*y^55 - 24059240265583258252231715/633825300114114700748351602688*y^54 + 7101564657998636156741025/158456325028528675187087900672*y^53 - 8316064668376595166616965/158456325028528675187087900672*y^52 + 151201175788665366665763/2475880078570760549798248448*y^51 - 350189595773086527650049/4951760157141521099596496896*y^50 + 404064918199715224211595/4951760157141521099596496896*y^49 - 465026356535455625543325/4951760157141521099596496896*y^48 + 130442175746270862705/1208925819614629174706176*y^47 - 1197968947370924359785/9671406556917033397649408*y^46 + 343802139976325885785/2417851639229258349412352*y^45 - 394972691042569738553/2417851639229258349412352*y^44 + 7101503107989975685/37778931862957161709568*y^43 - 16385419366240285361/75557863725914323419136*y^42 + 9486295422560165209/37778931862957161709568*y^41 - 11039378753773223931/37778931862957161709568*y^40 + 25261736278657263/73786976294838206464*y^39 - 119481185101757325/295147905179352825856*y^38 + 35771942691071587/73786976294838206464*y^37 - 43633908117680727/73786976294838206464*y^36 + 855566825836877/1152921504606846976*y^35 - 2255585268115403/2305843009213693952*y^34 + 6611198199648595/4611686018427387904*y^33 - 14544636039226909/4611686018427387904*y^32 + x^2 - y - 1
In [21]:
(gk_prime * gk) % y^(2**5)

Out[21]:
x^2 - y - 1

Aha! When we combine the two factors together and go modulo $y^{2^k}$, we do recover our original factor.

Therefore, this perhaps suggest the following:

If $f(x,y)$ was reducible, then this $g_k$ is supposed to be part of a proper factor $g^\ast(x,y)$ of $f$. Note that $\deg_x(g^\ast) < \deg_x(f)$ and $\deg_y(g^\ast) \leq \deg_y(f)$ (since $f$ is monic in $x$, each factor must be monic in $x$ too so the degree in $x$ must be smaller).

Claim(?): $f(x,y)$ is reducible if and only if we can find some $g_k'(x,y)$ such that $\tilde{g}(x,y) = g_k \cdot g_k' \bmod{y^{2^k}}$ satisfies $\deg_x(\tilde{g}) < \deg_x(f)$ and $\deg_y(\tilde{g}) \leq \deg_y(f)$.

Turns out, the above claim is indeed true, and we can actually test if such a $g_k'$ exists by just writing down a system of linear equations to solve.