235 lines
6.5 KiB
Plaintext
235 lines
6.5 KiB
Plaintext
{
|
|
"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*(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*(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.<sqrt2> = 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"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "10873881",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### 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
|
|
}
|