From 4fa8149aa4dc8fe7fdcda652da87227baacdf36d Mon Sep 17 00:00:00 2001 From: filifa Date: Wed, 9 Jul 2025 21:38:53 -0400 Subject: [PATCH] add problem 45 --- notebooks/problem0045.ipynb | 365 ++++++++++++++++++++++++++++++++++++ 1 file changed, 365 insertions(+) create mode 100644 notebooks/problem0045.ipynb diff --git a/notebooks/problem0045.ipynb b/notebooks/problem0045.ipynb new file mode 100644 index 0000000..b4012ac --- /dev/null +++ b/notebooks/problem0045.ipynb @@ -0,0 +1,365 @@ +{ + "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*(3*sqrt(3) + 5) - 1/12*(-2107560*sqrt(3) + 3650401)^t*(3*sqrt(3) - 5) + 1/6,\n", + " 1/24*sqrt(3)*((2107560*sqrt(3) + 3650401)^t*(3*sqrt(3) + 5) + (-2107560*sqrt(3) + 3650401)^t*(3*sqrt(3) - 5)) + 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)" + ] + } + ], + "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 +}