eulerbooks/notebooks/problem0100.ipynb

{
 "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
}