eulerbooks/notebooks/problem0066.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ca2d9c5a",
   "metadata": {},
   "source": [
    "# [Diophantine Equation](https://projecteuler.net/problem=66)\n",
    "\n",
    "SageMath can solve these [Diophantine equations](https://en.wikipedia.org/wiki/Diophantine_equation) for us. For equations of the form $x^2 - Dy^2 = 1$ (which are called [Pell equations](https://en.wikipedia.org/wiki/Pell%27s_equation)), it will give parameterizations in $t$ for both $x$ and $y$. Here's $D=2$ as an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d6c3ebd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-sqrt(2)*(2*sqrt(2) + 3)^t + sqrt(2)*(-2*sqrt(2) + 3)^t - 3/2*(2*sqrt(2) + 3)^t - 3/2*(-2*sqrt(2) + 3)^t,\n",
       "  3/4*sqrt(2)*(2*sqrt(2) + 3)^t - 3/4*sqrt(2)*(-2*sqrt(2) + 3)^t + (2*sqrt(2) + 3)^t + (-2*sqrt(2) + 3)^t),\n",
       " (sqrt(2)*(2*sqrt(2) + 3)^t - sqrt(2)*(-2*sqrt(2) + 3)^t + 3/2*(2*sqrt(2) + 3)^t + 3/2*(-2*sqrt(2) + 3)^t,\n",
       "  -3/4*sqrt(2)*(2*sqrt(2) + 3)^t + 3/4*sqrt(2)*(-2*sqrt(2) + 3)^t - (2*sqrt(2) + 3)^t - (-2*sqrt(2) + 3)^t)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('x,y')\n",
    "solve_diophantine(x^2 - 2*y^2 == 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "519c091f",
   "metadata": {},
   "source": [
    "$t=0$ seems to correspond to the minimal solution we are after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4b9b6240",
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimal_solution(d):\n",
    "    var('x,y')\n",
    "    \n",
    "    sols = solve_diophantine(x^2 - d*y^2 == 1)\n",
    "    u, v = sols[0]\n",
    "    return (abs(u(t=0).simplify_full()), abs(v(t=0).simplify_full()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42ba15f7",
   "metadata": {},
   "source": [
    "Now we can just iterate (although it is pretty slow)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "300b1afb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "661"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max((d for d in range(1, 1001) if not is_square(d)), key=lambda k: minimal_solution(k)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b16636f7",
   "metadata": {},
   "source": [
    "Let's dig a little deeper so we can optimize.\n",
    "\n",
    "## Solving Pell equations\n",
    "Lagrange proved that if $(x_0, y_0)$ is a solution to\n",
    "$$x^2 - dy^2 = 1$$\n",
    "then $\\frac{x_0}{y_0}$ is a [convergent of the continued fraction](https://en.wikipedia.org/wiki/Simple_continued_fraction) of $\\sqrt{d}$. This is great for us, since we can write generators for computing convergents of square roots (FYI, SageMath can do this with built-in methods: `continued_fraction(sqrt(d)).convergents()`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c51f9b7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def continued_fraction_sqrt(d):\n",
    "    x = sqrt(d)\n",
    "    while True:\n",
    "        b = floor(x)\n",
    "        yield b\n",
    "        x = (x - b)^-1\n",
    "        \n",
    "        \n",
    "def convergents(partial_denoms):\n",
    "    h, hprev = 1, 0\n",
    "    k, kprev = 0, 1\n",
    "    for b in partial_denoms:\n",
    "        h, hprev = b * h + hprev, h\n",
    "        k, kprev = b * k + kprev, k\n",
    "        yield h/k"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5cbf5eb",
   "metadata": {},
   "source": [
    "Then we can just iterate over each convergent to see if its numerator and denominator are a solution to the Pell equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d95125c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pell_fundamental_solution(d):\n",
    "    partial_denoms = continued_fraction_sqrt(d)\n",
    "    for f in convergents(partial_denoms):\n",
    "        x, y = f.as_integer_ratio()\n",
    "        if x^2 - d*y^2 == 1:\n",
    "            return (x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0bec674",
   "metadata": {},
   "source": [
    "Now we'll just iterate with this function instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "03d7ba9d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "661"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max((d for d in range(1, 1001) if not is_square(d)), key=lambda x: pell_fundamental_solution(x)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfb42d20",
   "metadata": {},
   "source": [
    "## Relevant sequences\n",
    "* Minimal values of $x$ for solutions to the Pell equation: [A002350](https://oeis.org/A002350)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.5",
   "language": "sage",
   "name": "sagemath"
  },
  "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}