{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R.<y,x> = PolynomialRing(QQ) \n",
    "\n",
    "# There is a reason for the order of variables\n",
    "# I'll later want to do f % y^i and for that sage\n",
    "# needs y to have \"higher precedence\" than x.\n",
    "#\n",
    "# ... just ignore this for now.\n",
    "\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extended_euclid(f , g, R = R):\n",
    "    old_d, d = f       , g\n",
    "    old_u, u = R.one() , R.zero()\n",
    "    old_v, v = R.zero(), R.one()\n",
    "    \n",
    "    while not d.is_zero():\n",
    "        q         = old_d // d\n",
    "        old_d , d = d , (old_d - q * d)\n",
    "        old_u , u = u , (old_u - q * u)\n",
    "        old_v , v = v , (old_v - q * v)\n",
    "    l = old_d.lc()\n",
    "    return (old_d / l, old_u / l, old_v / l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HenselLift(f, g, h, a, b, y_mon, R = R, monic = False):\n",
    "    '''\n",
    "    Expecting polynomials f, g, h, a, b, y_mon\n",
    "    from the ring R such that \n",
    "    *  f - gh  = 1 mod y_mon\n",
    "    *  ag + bh = 1 mod y_mon\n",
    "    and y_mon is y^i for some i. \n",
    "    \n",
    "    It will return polynomials g', h', a', b' such that\n",
    "    *  g = g'   mod y_mon^2\n",
    "    *  h = h'   mod y_mon^2\n",
    "    *  f = g'h' mod y_mon^2\n",
    "    *  a'g' + b' h' = 1 mod y_mon^2\n",
    "    '''\n",
    "    \n",
    "    if (f - g*h) % y_mon != R.zero(): raise ValueError(\"f - gh is not 0 mod I\")\n",
    "    if (a*g + b*h) % y_mon != R.one(): raise ValueError(\"ag + bh is not 1 mod I\")\n",
    "        \n",
    "    gamma = f - g * h\n",
    "    g1 = g + (b * gamma)\n",
    "    h1 = h + (a * gamma)\n",
    "    \n",
    "    if monic:\n",
    "        q , r = (g1 - g) // g , (g1 - g) % g\n",
    "        g1 = (g + r) % y_mon^2\n",
    "        h1 = (h1 * ( 1 + q)) % y_mon^2\n",
    "    \n",
    "    delta = a*g1 + b*h1 - R.one()\n",
    "    a1 = a*(1 - delta)\n",
    "    b1 = b*(1 - delta)\n",
    "    return (g1, h1, a1, b1)\n",
    "\n",
    "def iteratedHenselLift(f, g, h, a, b, y_mon, k, R = R, monic = False):\n",
    "    if k == 1:\n",
    "        return HenselLift(f, g, h, a, b, y_mon, R, monic)\n",
    "    else:\n",
    "        g1, h1, a1, b1 = iteratedHenselLift(f, g, h, a, b, y_mon, k-1, R, monic)\n",
    "        y_mon_1 = y^(2**(k-1))\n",
    "        return HenselLift(f,g1, h1, a1, b1, y_mon_1, R, monic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1\n",
    "\n",
    "\\begin{align*}\n",
    "f(x,y) & = x^2 - 2xy - 3y^2 + 3x - 5y + 2\\\\\n",
    "       & = (x + y + 2)(x - 3y + 1)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x^2 - 2*x*y - 3*y^2 + 3*x - 5*y + 2\n",
    "f % y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor(f % y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = x + 1\n",
    "h = x + 2\n",
    "_ , a, b = extended_euclid(g, h)\n",
    "a*g + b*h == 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gk, hk, ak, bk = iteratedHenselLift(f, g, h, a, b, y, k = 5, monic = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gk)\n",
    "print(hk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, we found our factors!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2\n",
    "\n",
    "Let's take $f(x,y) = x^3 + x - y$ which is actually irreducible. \n",
    "\n",
    "However, $f \\bmod{y} = x^3 + x = x (x^2 + 1)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x^3 + x - y\n",
    "f0 = (f % y)\n",
    "factor(f0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = x\n",
    "h = x^2 + 1\n",
    "_, a,b = extended_euclid(g, h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gk, hk, ak, bk = iteratedHenselLift(f, g, h, a, b, y, k = 5, monic = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We appear to be getting some garbage. Perhaps this suggests that we started off with an $f(x,y)$ that was irreducible. \n",
    "\n",
    "(Spoiler: This intuition is wrong.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3\n",
    "\n",
    "\n",
    "$$f(x,y)  = x^4 - y^2 - 2y - 1$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x^4 - y^2 - 2*y - 1\n",
    "factor(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0 = f % y\n",
    "factor(f0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = x - 1\n",
    "h = (x + 1)*(x^2 + 1)\n",
    "_, a,b = extended_euclid(g, h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gk, hk, ak, bk = iteratedHenselLift(f, g, h, a, b, y, k = 5, monic = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"gk: {gk}\")\n",
    "print(\"\")\n",
    "print(f\"hk: {hk}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Garbage again. Suggests that $f(x,y)$ is irreducible, right?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hunh?!?!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Making sense of this\n",
    "\n",
    "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. \n",
    "\n",
    "But how do we really distinguish the mess in Example 2 and Example 3?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose on the other hand, we did the entire lifting process starting with $g_1 = x +1$ instead, what would we end up with?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = (x + 1)\n",
    "h1 = (x - 1) * (x^2 + 1)\n",
    "_, a1, b1 = extended_euclid(g1, h1)\n",
    "gk_prime, _, _, _ = iteratedHenselLift(f, g1, h1, a1, b1, y, k=5, monic = True)\n",
    "gk_prime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "However, $g_k'$ and $g_k$ together are part of a proper factor of $f$. What is their product going to be?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gk_prime * gk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(gk_prime * gk) % y^(2**5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aha! When we combine the two factors together and go modulo $y^{2^k}$, we _do_ recover our original factor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, this perhaps suggest the following:\n",
    "\n",
    "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). \n",
    "\n",
    "**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)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.4",
   "language": "sage",
   "name": "sagemath-9.4"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}