148 lines
4.2 KiB
Plaintext
148 lines
4.2 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "721eef2a",
|
|
"metadata": {},
|
|
"source": [
|
|
"# [Odd Period Square Roots](https://projecteuler.net/problem=64)\n",
|
|
"\n",
|
|
"SageMath can compute the periods of these continued fractions if you ask it nicely, but it's pretty slow since we have to construct a new quadratic field on every loop."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"id": "d9f0863f",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"1322"
|
|
]
|
|
},
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"odd_period_sqrts = set()\n",
|
|
"for N in range(1, 10001):\n",
|
|
" if is_square(N):\n",
|
|
" continue\n",
|
|
" \n",
|
|
" K.<s> = QuadraticField(N)\n",
|
|
" if continued_fraction(s).period_length() % 2 == 1:\n",
|
|
" odd_period_sqrts.add(N) \n",
|
|
" \n",
|
|
"len(odd_period_sqrts)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "a67c8448",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can speed it up if we apply an algorithm for computing the [repetend of continued fractions of square roots](https://en.wikipedia.org/wiki/Periodic_continued_fraction)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"id": "94fd3a59",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def repetend_sqrt(S):\n",
|
|
" m, d, a0 = 0, 1, isqrt(S)\n",
|
|
" \n",
|
|
" repetend = []\n",
|
|
" a = a0\n",
|
|
" while a != 2 * a0:\n",
|
|
" m = d * a - m\n",
|
|
" d = (S - m^2) / d\n",
|
|
" a = floor((a0 + m) / d)\n",
|
|
" repetend.append(a)\n",
|
|
" \n",
|
|
" return repetend"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "67b9827d",
|
|
"metadata": {},
|
|
"source": [
|
|
"To gain an intuition for this algorithm (without getting too deep into the weeds), it's helpful to consider the more general problem of finding the continued fraction of a [quadratic irrational](https://en.wikipedia.org/wiki/Quadratic_irrational_number). To do this, consider the equation\n",
|
|
"$$\\frac{m+\\sqrt{S}}{d} = a_0 + \\frac{1}{r_1}$$\n",
|
|
"The first term of the continued fraction, $a_0$, is simply the floor of the left hand side, and the following terms are just the terms of the continued fraction of $r_1$. We can then use some algebra to solve for $r_1$.\n",
|
|
"\n",
|
|
"$$r_1 = \\frac{d\\left(da_0 - m + \\sqrt{S}\\right)}{S-(m-da_0)^2}$$\n",
|
|
"\n",
|
|
"This shows that $r_1$ is *also* a quadratic irrational, so we can apply this process repeatedly to determine successive terms of the continued fraction.\n",
|
|
"\n",
|
|
"It is also known that the last term in the period of the continued fraction of a square root is twice the initial term, which provides a stopping condition."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"id": "b8490f58",
|
|
"metadata": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"1322"
|
|
]
|
|
},
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"odd_period_sqrts = {N for N in range(1, 10001) if not is_square(N) and len(repetend_sqrt(N)) % 2 == 1}\n",
|
|
"len(odd_period_sqrts)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "f6ba21d9",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Relevant sequences\n",
|
|
"* Periods of continued fractions of $\\sqrt{n}$: [A003285](https://oeis.org/A003285)\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
|
|
}
|