"# [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:"
"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."
"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",
"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."
"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",
"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",
"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)."
