diff --git a/notebooks/problem0100.ipynb b/notebooks/problem0100.ipynb new file mode 100644 index 0000000..59ec40d --- /dev/null +++ b/notebooks/problem0100.ipynb @@ -0,0 +1,224 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4945116b", + "metadata": {}, + "source": [ + "# [Arranged Probability](https://projecteuler.net/problem=100)\n", + "\n", + "We can rephrase this problem as finding integer solutions to\n", + "$$\\frac{b}{t} \\times \\frac{b-1}{t-1} = \\frac{1}{2}$$\n", + "Specifically, we're looking for the first solution with $t > 10^{12}$. This is a [Diophantine problem](https://en.wikipedia.org/wiki/Diophantine_equation). " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "fa08d031", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(-1/8*(12*sqrt(2) + 17)^t*(sqrt(2) + 2) + 1/8*(-12*sqrt(2) + 17)^t*(sqrt(2) - 2) + 1/2,\n", + " -1/8*sqrt(2)*((12*sqrt(2) + 17)^t*(sqrt(2) + 2) + (-12*sqrt(2) + 17)^t*(sqrt(2) - 2)) + 1/2),\n", + " (-1/8*(12*sqrt(2) + 17)^t*(7*sqrt(2) + 10) + 1/8*(-12*sqrt(2) + 17)^t*(7*sqrt(2) - 10) + 1/2,\n", + " -1/8*sqrt(2)*((12*sqrt(2) + 17)^t*(7*sqrt(2) + 10) + (-12*sqrt(2) + 17)^t*(7*sqrt(2) - 10)) + 1/2),\n", + " (1/8*(12*sqrt(2) + 17)^t*(sqrt(2) + 2) - 1/8*(-12*sqrt(2) + 17)^t*(sqrt(2) - 2) + 1/2,\n", + " 1/8*sqrt(2)*((12*sqrt(2) + 17)^t*(sqrt(2) + 2) + (-12*sqrt(2) + 17)^t*(sqrt(2) - 2)) + 1/2),\n", + " (1/8*(12*sqrt(2) + 17)^t*(7*sqrt(2) + 10) - 1/8*(-12*sqrt(2) + 17)^t*(7*sqrt(2) - 10) + 1/2,\n", + " 1/8*sqrt(2)*((12*sqrt(2) + 17)^t*(7*sqrt(2) + 10) + (-12*sqrt(2) + 17)^t*(7*sqrt(2) - 10)) + 1/2)]" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "var('b,t')\n", + "\n", + "sols = sorted(solve_diophantine(b/t * (b-1)/(t-1) == 1/2))\n", + "sols" + ] + }, + { + "cell_type": "markdown", + "id": "59126a72", + "metadata": {}, + "source": [ + "Inspecting these parameterizations, the first two will always generate negative solutions for integer $t$, so we'll only consider the last two. We'll evaluate them at $t=0,1,2,\\ldots$ to find the first instance where the total number of discs is greater than $10^{12}$." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "4f211b1a", + "metadata": {}, + "outputs": [], + "source": [ + "from itertools import count\n", + "\n", + "candidates = []\n", + "for sol in sols[2:4]:\n", + " for n in count(0):\n", + " b, t = sol[0](t=n), sol[1](t=n)\n", + " if t > 10^12:\n", + " candidates.append((b, t))\n", + " break" + ] + }, + { + "cell_type": "markdown", + "id": "ba65bcf6", + "metadata": {}, + "source": [ + "Then we'll select the smaller $t$ of the two." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9c06a2c4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(756872327473, 1070379110497), (4411375203411, 6238626641380)]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "[(int(b), int(t)) for (b, t) in candidates]" + ] + }, + { + "cell_type": "markdown", + "id": "15f57948", + "metadata": {}, + "source": [ + "This tells us that $b=756872327473, t=1070379110497$ is the arrangement we're after.\n", + "\n", + "## Alternate method\n", + "However, we don't need SageMath's `solve_diophantine` for this problem. In fact, we can solve it using the same techniques as [problem 45](https://projecteuler.net/problem=45).\n", + "\n", + "First, rearrange:\n", + "$$2b^2 - 2b - t^2 + t = 0$$\n", + "Then complete the square:\n", + "$$(2t-1)^2 - 2(2b-1)^2 = -1$$\n", + "Then substitute $X=2t-1$ and $Y=2b-1$ to get a Pell equation:\n", + "$$X^2 - 2Y^2 = -1$$\n", + "The problem gives $b=85, t=120$ as a solution to the original equation; this corresponds to $X=239, Y=169$ for our Pell equation. To generate more solutions, we'll need a fundamental solution to this Pell equation, which we can find with SageMath, or with continued fractions, as in [problem 66](https://projecteuler.net/problem=66)." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "0349d6a5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "sqrt2 + 1" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "K. = QQ[sqrt(2)]\n", + "u, *_ = K.units()\n", + "u" + ] + }, + { + "cell_type": "markdown", + "id": "2e30da3f", + "metadata": {}, + "source": [ + "With both of these values in hand, we can easily generate more solutions to $X^2 - 2Y^2 = -1$ by computing $\\left(239+169\\sqrt{2}\\right)\\left(1+\\sqrt{2}\\right)^k$ for successive values of $k$. Then we'll check to see if they correspond to solutions to our original equation." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "989c878e", + "metadata": {}, + "outputs": [], + "source": [ + "z = 239 + 169*sqrt2\n", + "while True:\n", + " z *= u\n", + " if z.norm() != -1:\n", + " continue\n", + " \n", + " x, y = z.parts()\n", + " if (x + 1) % 2 != 0 or (y + 1) % 2 != 0:\n", + " continue\n", + " \n", + " t, b = (x + 1) // 2, (y + 1) // 2\n", + " if t > 10^12:\n", + " break" + ] + }, + { + "cell_type": "markdown", + "id": "02b30f6e", + "metadata": {}, + "source": [ + "As expected, we get the same solution as above." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "ac482fb0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(756872327473, 1070379110497)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "b, t" + ] + } + ], + "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 +}