217 lines
6.0 KiB
Plaintext
217 lines
6.0 KiB
Plaintext
{
|
|
"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",
|
|
" sols = solve_diophantine(x^2 - d*y^2 == 1)\n",
|
|
" u, v = sols[0]\n",
|
|
" return (abs(u(t=0)), abs(v(t=0)))"
|
|
]
|
|
},
|
|
{
|
|
"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}$. Specifically, if $p$ is the period of the continued fraction, then the first solution will be the $(p-1)$th convergent if $p$ is even, and the $(2p-1)$th convergent if $p$ is odd.\n",
|
|
"\n",
|
|
"This is great for us, since there are algorithms to compute these convergents. We'll use SageMath here; see [problem 64](https://projecteuler.net/problem=64) for how to compute the partial denominators of the continued fraction of a square root, and see [problem 65](https://projecteuler.net/problem=65) for an algorithm that uses partial denominators to compute convergents of continued fractions (SageMath's constructions make this implementation a little slow, but it makes the code easier to read - and it's still considerably faster than using `solve_diophantine`)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"id": "5d95125c",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def pell_fundamental_solution(d):\n",
|
|
" K.<s> = QuadraticField(d)\n",
|
|
" f = continued_fraction(s)\n",
|
|
" p = f.period_length()\n",
|
|
" if p % 2 == 0:\n",
|
|
" n = p - 1\n",
|
|
" else:\n",
|
|
" n = 2 * p - 1\n",
|
|
" \n",
|
|
" x, y = f.convergent(n).as_integer_ratio()\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": 5,
|
|
"id": "03d7ba9d",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"661"
|
|
]
|
|
},
|
|
"execution_count": 5,
|
|
"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": "3f573e65",
|
|
"metadata": {},
|
|
"source": [
|
|
"And in case you want to know the minimal $x$ for $x^2 - 661y^2 = 1$:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"id": "b4b119a0",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(16421658242965910275055840472270471049, 638728478116949861246791167518480580)"
|
|
]
|
|
},
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"pell_fundamental_solution(661)"
|
|
]
|
|
},
|
|
{
|
|
"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)\n",
|
|
"\n",
|
|
"#### Copyright (C) 2025 filifa\n",
|
|
"\n",
|
|
"This work is licensed under the [Creative Commons Attribution-ShareAlike 4.0 International license](https://creativecommons.org/licenses/by-sa/4.0/) and the [BSD Zero Clause license](https://spdx.org/licenses/0BSD.html)."
|
|
]
|
|
}
|
|
],
|
|
"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
|
|
}
|