filifa
/
eulerbooks
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity
eulerbooks/notebooks/problem0045.ipynb

370 lines
16 KiB
Plaintext
Raw Permalink Blame History

{
"cells": [
{
"cell_type": "markdown",
"id": "4f36b586",
"metadata": {},
"source": [
"# [Triangular, Pentagonal, and Hexagonal](https://projecteuler.net/problem=45)\n",
"\n",
"[All hexagonal numbers are triangular numbers](https://en.wikipedia.org/wiki/Polygonal_number). Therefore this problem reduces to finding numbers that are both pentagonal and hexagonal.\n",
"\n",
"One way to frame this problem is that we are finding integer solutions to\n",
"$$\\frac{x(3x-1)}{2} = y(2y-1)$$\n",
"In other words, we are solving a [Diophantine equation](https://en.wikipedia.org/wiki/Diophantine_equation), which is no easy task, in general. Fortunately, there are methods to tackle equations of this form, and SageMath knows them."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "bae28010",
"metadata": {},
"outputs": [],
"source": [
"var('x,y')\n",
"\n",
"P(n) = n * (3*n - 1) / 2\n",
"H(n) = n * (2*n - 1)\n",
"\n",
"params = solve_diophantine(P(x) == H(y))"
]
},
{
"cell_type": "markdown",
"id": "f87cb6b3",
"metadata": {},
"source": [
"Here, `params` ends up being a whole bunch of ugly parameterizations in $t$ for $x$ and $y$ satisfying $P(x) = H(y)$. Here's one, just as an example:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "2c65dbea",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-1/12*(2107560*sqrt(3) + 3650401)^t*(153*sqrt(3) + 265) + 1/12*(-2107560*sqrt(3) + 3650401)^t*(153*sqrt(3) - 265) + 1/6,\n",
" -1/24*sqrt(3)*((2107560*sqrt(3) + 3650401)^t*(153*sqrt(3) + 265) + (-2107560*sqrt(3) + 3650401)^t*(153*sqrt(3) - 265)) + 1/4)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"params[0]"
]
},
{
"cell_type": "markdown",
"id": "d4ef92eb",
"metadata": {},
"source": [
"We want to check all of these parameterizations, because one of them might generate the next pentagonal hexagonal number. We'll guess a small range of values to try for $t$ (which turns out to be all we need). We're also only interested in positive indices, but some of these parameterizations generate negative values, so we'll filter those out."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "7746c170",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[1,\n",
" 40755,\n",
" 1533776805,\n",
" 57722156241751,\n",
" 2172315626468283465,\n",
" 81752926228785223683195,\n",
" 3076689623521787481625080301,\n",
" 115788137209866023854693048367775,\n",
" 4357570752679408318225730700647767185,\n",
" 163992817590548715438241125333485021875651,\n",
" 6171705692845139604123358192574644612620485685,\n",
" 232265971880541166271029746781113050017874336396775,\n",
" 8741097579580580558598793886237050331798038163335747801,\n",
" 328962466077669596861765842843615405405774318221103196349195,\n",
" 12380173439625920028715115170977828280803860360134959528069859965]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pent_hexes = dict()\n",
"for m, n in params:\n",
" for t in range(-5, 5):\n",
" u, v = m(t=t), n(t=t)\n",
" if u <= 0 or v <= 0:\n",
" continue\n",
" \n",
" pent_hexes[(u, v)] = P(u).simplify_full()\n",
"\n",
"sorted(pent_hexes.values())"
]
},
{
"cell_type": "markdown",
"id": "3f115b65",
"metadata": {},
"source": [
"This shows the next number in the sequence as 1533776805. (This number also shows up rather quickly with a few web searches about the pentagonal hexagonal numbers, but isn't it more satisfying to compute it yourself?)\n",
"\n",
"Of course, there's a lot of interesting math we're glossing over by having SageMath solve the hardest part, so let's look at how we find Diophantine solutions.\n",
"\n",
"## Transformation to Pell equation\n",
"It's convenient to rearrange our equation as\n",
"$$3x^2 - x - 4y^2 + 2y = 0$$\n",
"Then, by [completing the square](https://en.wikipedia.org/wiki/Completing_the_square) on both $x$ and $y$, we can transform this to\n",
"$$(6x-1)^2 - 3(4y-1)^2 = -2$$\n",
"If we make the substitutions $X = 6x - 1$ and $Y = 4y - 1$,\n",
"$$X^2 - 3Y^2 = -2$$\n",
"then we have a generalized [Pell equation](https://en.wikipedia.org/wiki/Pell%27s_equation), which is a special category of Diophantine equations that has known solution methods.\n",
"\n",
"It's important to note the implications of making these substitutions. If we think of $X$ and $Y$ as functions of $x$ and $y$, respectively, then they are [injective](https://en.wikipedia.org/wiki/Injective_function). This means that every distinct solution to our original equation maps to a distinct solution to our Pell equation, as well, so we will not \"lose\" any solutions by making these substitutions. However, these functions are not [surjective](https://en.wikipedia.org/wiki/Surjective_function), since some integers can not be expressed as $6x-1$ or $4y-1$ with integer $x,y$. This means a solution to our general Pell equation does not necessarily correspond to a solution to our original equation. For example, $X=3691, Y=2131$ is a solution for $X^2 - 3Y^2 = -2$, but these values correspond to $x = \\frac{1846}{3}, y=533$, which are not solutions to our original equation since they aren't both integers.\n",
"\n",
"The main takeaway from this should be that when we find the solutions to $X^2 - 3Y^2 = -2$, we should check that they actually correspond to integer solutions to $3x^2 - x - 4y^2 + 2y = 0$, as some will be extraneous.\n",
"\n",
"## But first, we need to talk about algebraic number theory\n",
"[And if you thought my other tangents were complicated...](https://www.youtube.com/watch?v=kpk2tdsPh0A)\n",
"\n",
"We could jump right in from here with the methods to solve Pell equations, but to understand *why* these methods work, it's helpful to establish a few facts.\n",
"\n",
"Consider that the equation $x^2 - dy^2 = n$ can be factored into\n",
"$$\\left(x + y\\sqrt{d}\\right)\\left(x - y\\sqrt{d}\\right) = n$$\n",
"Also recall that we are looking for integer values for $x$ and $y$ when solving Pell equations. Numbers of the form $x + y\\sqrt{d}$ where $x$ and $y$ are integers (and $d$ is a non-square integer) are [quadratic integers](https://en.wikipedia.org/wiki/Quadratic_integer) (there are other quadratic integers, too, but they have non-integer $x$ and $y$, so they don't correspond to solutions to Pell equations).\n",
"\n",
"Numbers of the form $x + y\\sqrt{d}$ where $x,y$ are rational (instead of just integers) form a [quadratic field](https://en.wikipedia.org/wiki/Quadratic_field), denoted $\\mathbb{Q}\\left(\\sqrt{d}\\right)$.\n",
"\n",
"The *norm* of a quadratic integer $x + y\\sqrt{d}$ is given by $x^2 - dy^2$. Therefore, in the language of [algebraic number theory](https://en.wikipedia.org/wiki/Algebraic_number_theory), finding integer solutions to $x^2 - dy^2 = n$ can be rephrased as finding quadratic integers in $\\mathbb{Q}\\left(\\sqrt{d}\\right)$ with norm $n$. Because of this, a solution $(x, y)$ to a Pell equation is sometimes written as $x + y\\sqrt{d}$.\n",
"\n",
"Furthermore, let $a = x_1 + y_1\\sqrt{d}$ and let $b = x_2 + y_2\\sqrt{d}$. Through some tedious algebra, you can show that the norm of $ab$ is equal to the product of the norm of $a$ and the norm of $b$, i.e. the norm is [completely multiplicative](https://en.wikipedia.org/wiki/Completely_multiplicative_function)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "b97af696",
"metadata": {},
"outputs": [],
"source": [
"var('x1,y1,x2,y2,d')\n",
"\n",
"a = x1 + y1*sqrt(d)\n",
"b = x2 + y2*sqrt(d)\n",
"\n",
"x = x1*x2 + y1*y2*d\n",
"y = x1*y2 + x2*y1\n",
"\n",
"assert a * b == x + y*sqrt(d)\n",
"\n",
"a_norm = x1^2 - d*y1^2\n",
"b_norm = x2^2 - d*y2^2\n",
"\n",
"ab_norm = x^2 - d*y^2\n",
"\n",
"assert a_norm * b_norm == ab_norm"
]
},
{
"cell_type": "markdown",
"id": "ff78b3ea",
"metadata": {},
"source": [
"Finally, note that a quadratic integer with a norm of 1 or -1 is called a *unit*. [Dirichlet's unit theorem](https://en.wikipedia.org/wiki/Dirichlet%27s_unit_theorem) leads to a powerful result about units of real quadratic integers: there exists a [fundamental unit](https://en.wikipedia.org/wiki/Fundamental_unit_(number_theory)) that generates all other units. In other words, there exists a unit $x + y\\sqrt{d}$ in $\\mathbb{Q}\\left(\\sqrt{d}\\right)$ such that every other unit can be expressed as $\\left(x + y\\sqrt{d}\\right)^k$ for some integer $k$.\n",
"\n",
"As an example, we'll have SageMath compute a fundamental unit of $\\mathbb{Q}\\left(\\sqrt{3}\\right)$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "d25dc6a1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"sqrt3 - 2"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"K = QQ[sqrt(3)]\n",
"u, *_ = K.units()\n",
"u"
]
},
{
"cell_type": "markdown",
"id": "730cb762",
"metadata": {},
"source": [
"$u$ has norm 1, and as we expect, powers of $u$ also have norm 1."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "927c1701",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{-4*sqrt3 + 7: 1, 15*sqrt3 - 26: 1, -56*sqrt3 + 97: 1, 209*sqrt3 - 362: 1}"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"{u^k: (u^k).norm() for k in range(2, 6)}"
]
},
{
"cell_type": "markdown",
"id": "5fa25571",
"metadata": {},
"source": [
"We can apply these properties of norms and units to help solve Pell equations.\n",
"\n",
"## Finding solutions to the general Pell equation\n",
"There are methods for finding solutions to a general Pell equation given only the equation, such as the Lagrange-Matthews-Mollin algorithm. However, the problem statement gives us that $P_{165} = H_{143}$. Therefore $x=165, y=143$ is a solution to our original equation; this corresponds to $989 + 571\\sqrt{3}$ for our general Pell equation $X^2 - 3Y^2 = -2$.\n",
"\n",
"Knowing this information lets us take a shortcut: if $X + Y\\sqrt{d}$ is a solution to the *general* Pell equation $x^2 - dy^2 = n$ (i.e. $X + Y\\sqrt{d}$ has norm $n$), and $x_0 + y_0\\sqrt{d}$ is a solution to the *standard* Pell equation $x^2 - dy^2 = 1$ (i.e. $x_0 + y_0\\sqrt{d}$ has norm 1), then\n",
"$$\\left(X + Y\\sqrt{d}\\right)\\left(x_0 + y_0\\sqrt{d}\\right)$$\n",
"gives *another* solution to $x^2 - dy^2 = n$ (this follows from the norm of quadratic integers being completely multiplicative). Since we have $989 + 571\\sqrt{3}$ as a solution to $x^2 - 3y^2 = -2$, we can multiply it by any unit of $\\mathbb{Q}\\left(\\sqrt{3}\\right)$ with norm 1 to get another solution. Furthermore, if we can find a *fundamental unit* $u$ of $\\mathbb{Q}\\left(\\sqrt{3}\\right)$ with norm 1, we can generate *every successive solution* to $x^2 - 3y^2 = -2$ by evaluating $\\left(989 + 571\\sqrt{3}\\right)u^k$ for successive values of $k$.\n",
"\n",
"We've now reduced the problem to finding a fundamental solution to the standard Pell equation $x^2 - 3y^2 = 1$, which we happened to already find above: $-2 + \\sqrt{3}$ (see [problem 66](https://projecteuler.net/problem=66) for *even more discussion* on finding solutions to standard Pell equations).\n",
"\n",
"## Solution to the original problem\n",
"We now have everything we need to solve the original problem:\n",
"* an initial solution to $x^2 - 3y^2 = -2$ ($989 + 571\\sqrt{3}$)\n",
"* a fundamental solution to $x^2 - 3y^2 = 1$ ($-2 + \\sqrt{3}$)\n",
"\n",
"To generate more solutions to $x^2 - 3y^2 = -2$, we can just repeatedly multiply $989 + 571\\sqrt{3}$ by $-2 + \\sqrt{3}$. However, since we want the next largest pentagonal hexagonal number, we will instead multiply by $2 + \\sqrt{3}$, which will also generate solutions.\n",
"\n",
"As a reminder, for each one of these solutions $X + Y\\sqrt{3}$, we'll need to check if $X = 6x-1$ and $Y = 4y-1$ for any integers $x$ and $y$. For example, $\\left(989 + 571\\sqrt{3}\\right)\\left(2 + \\sqrt{3}\\right) = 3691 + 2131\\sqrt{3}$. As we noted above, $(3691, 2131)$ is a solution to $X^2 - 3Y^2 = -2$, but does *not* correspond to a solution to $3x^2 - x - 4y^2 + 2y = 0$.\n",
"\n",
"Iterating this way, we eventually get $X=191861, Y=110771$, which corresponds to the integer solutions $x=31977, y=27693$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "6257dd3e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(31977, 27693)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x0, y0 = 2, 1 \n",
"X, Y = 989, 571\n",
"while True:\n",
" X, Y = x0*X + 3*y0*Y, x0*Y + y0*X\n",
" if (X + 1) % 6 == 0 and (Y + 1) % 4 == 0:\n",
" x, y = (X + 1) // 6, (Y + 1) // 4\n",
" break\n",
"\n",
"(x, y)"
]
},
{
"cell_type": "markdown",
"id": "993746fd",
"metadata": {},
"source": [
"As expected, $P_x = H_y = 1533776805$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "b9342827",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1533776805"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P(x)"
]
},
{
"cell_type": "markdown",
"id": "74cb1521",
"metadata": {},
"source": [
"## Further reading\n",
"* [Solving the Pell Equation](https://www.ams.org/notices/200202/fea-lenstra.pdf)\n",
"* [The Diophantine Equation $x^2 - Dy^2 = N, D>0$](http://www.numbertheory.org/PDFS/patz5.pdf)\n",
"* [Solving the Diophantine Equation $ax^2 + bxy + cy^2 + dx + ey + f = 0$](http://www.numbertheory.org/PDFS/general_quadratic_solution.pdf)\n",
"* [Solving the generalized Pell equation $x^2 - Dy^2 = N$](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=5ac34a344ee346855184ff949eeaed18685b155c)\n",
"* [Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0](https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf)\n",
"* [Pell's Equation, I](https://kconrad.math.uconn.edu/blurbs/ugradnumthy/pelleqn1.pdf)\n",
"* [Pell's Equation, II](https://kconrad.math.uconn.edu/blurbs/ugradnumthy/pelleqn2.pdf)\n",
"\n",
"## Relevant sequences\n",
"* Hexagonal pentagonal numbers: [A046180](https://oeis.org/A046180)\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
}
Reference in New Issue View Git Blame Copy Permalink