From 275bfae6481fc43e6148e4b90f2f8a8846354c53 Mon Sep 17 00:00:00 2001 From: filifa Date: Thu, 15 May 2025 21:16:58 -0400 Subject: [PATCH] add problem 66 --- notebooks/problem0066.ipynb | 210 ++++++++++++++++++++++++++++++++++++ 1 file changed, 210 insertions(+) create mode 100644 notebooks/problem0066.ipynb diff --git a/notebooks/problem0066.ipynb b/notebooks/problem0066.ipynb new file mode 100644 index 0000000..ab729eb --- /dev/null +++ b/notebooks/problem0066.ipynb @@ -0,0 +1,210 @@ +{ + "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 +}