{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e0ff1edf",
   "metadata": {},
   "source": [
    "# [Pentagon Numbers](https://projecteuler.net/problem=44)\n",
    "\n",
    "It's not difficult to test if a number is [pentagonal](https://en.wikipedia.org/wiki/Pentagonal_number)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "76d738b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "\n",
    "def is_pentagonal(x):\n",
    "    n = (sqrt(24*x+1) + 1) / 6\n",
    "    return n == int(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04ec4508",
   "metadata": {},
   "source": [
    "With this function, and given a value $j$, we can attempt to find a value $k$ such that $P_j + P_k$ and $P_j - P_k$ are both pentagonal by just iterating over all $k < j$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3b7bfb29",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pentagonal_pair(j):\n",
    "    a = polygonal_number(5, j)\n",
    "    for k in reversed(range(1, j)):\n",
    "        b = polygonal_number(5, k)\n",
    "        if is_pentagonal(a + b) and is_pentagonal(a - b):\n",
    "            return k\n",
    "        \n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b9d6594",
   "metadata": {},
   "source": [
    "A naive approach that happens to work is to just iterate over values of $j$ until this function returns a solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "182f53fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2167, 1020)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from itertools import count\n",
    "\n",
    "for j in count(2):\n",
    "    k = pentagonal_pair(j)\n",
    "    if k is not None:\n",
    "        break\n",
    "\n",
    "j, k"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac17792a",
   "metadata": {},
   "source": [
    "We can then recover the pentagonal numbers from these indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1ec74c04",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7042750, 1560090)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a, b = polygonal_number(5, j), polygonal_number(5, k)\n",
    "a, b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8756e09",
   "metadata": {},
   "source": [
    "And calculate their difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "42e1eb38",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5482660"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a - b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54b24ef2",
   "metadata": {},
   "source": [
    "This turns out to be the answer! The thing is... we sort of just got lucky that the first $j,k$ pair we found happened to minimize $|P_j - P_k|$. If you care to find out how we can *verify* that this is the minimal solution, read on.\n",
    "\n",
    "## Diophantine approach\n",
    "Let's restate the problem. We are looking for positive integers $j, k, s, d$ such that\n",
    "$$P_j + P_k = P_s$$\n",
    "$$P_j - P_k = P_d$$\n",
    "\n",
    "Suppose $d$ is some constant. Then we can solve for $j$ and $k$ in the second equation using Diophantine methods. SageMath can do this for us with [solve_diophantine](https://doc.sagemath.org/html/en/reference/calculus/sage/symbolic/expression.html), but we can do it more efficiently ourselves.\n",
    "\n",
    "$P_j - P_k = P_d$ expands to\n",
    "$$3j^2 - j - 3k^2 + k = 3d^2 - d$$\n",
    "We can factor the left hand side into\n",
    "$$(j-k)(3j + 3k - 1) = 3d^2 - d$$\n",
    "This means that any solution $j,k$ corresponds to a factorization of $3d^2 - d$. Importantly, this means there's a finite number of solutions for a fixed $d$.\n",
    "\n",
    "So if we factor $3d^2 - d = ab$, and set $a=j-k$ and $b=3j+3k-1$, we can solve for $j$ and $k$:\n",
    "$$k = \\frac{b+1-3a}{6}$$\n",
    "$$j = a + k$$\n",
    "If we iterate over all possible factorizations $ab$ of $3d^2 - d$, this will let us find all $j,k$ that solve $P_j - P_k = P_d$.\n",
    "\n",
    "Once we have found candidate pairs for $j$ and $k$, we can just check if $P_j + P_k$ is also a pentagonal number using the function from above - if so, $j$ and $k$ also satisfy $P_j + P_k = P_s$, and we've found a solution for our fixed value of $d$.\n",
    "\n",
    "Here's a function that applies all of this to find solutions (if they exist) for a given $d$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "21b7a57d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pentagonal_sums(d):\n",
    "    psums = set()\n",
    "    N = 3*d^2 - d\n",
    "    for a in divisors(N):\n",
    "        b = N / a\n",
    "        \n",
    "        k = (b + 1 - 3*a) / 6\n",
    "        if k != int(k) or k <= 0:\n",
    "            continue\n",
    "            \n",
    "        j = a + k\n",
    "            \n",
    "        p = polygonal_number(5, j) + polygonal_number(5, k)\n",
    "        if is_pentagonal(p):\n",
    "            psums.add((j, k))\n",
    "    \n",
    "    return psums"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bb82c01",
   "metadata": {},
   "source": [
    "Since we want to minimize $P_d$, we can just iterate over increasing values of $d$ and stop once this function finds a solution. This way, we know our solution is minimal since we didn't find solutions for any smaller $d$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "07d2dcc4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(2167, 1020)}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for d in count(1):\n",
    "    sols = pentagonal_sums(d)\n",
    "    if sols != set():\n",
    "        break\n",
    "            \n",
    "sols"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1af2fcc4",
   "metadata": {},
   "source": [
    "This confirms that the solution we found earlier is in fact the solution that minimizes $P_d$.\n",
    "\n",
    "## Relevant sequences\n",
    "* Pentagonal numbers: [A000326](https://oeis.org/A000326)\n",
    "* Pentagonal numbers which are the sum of two other positive pentagonal numbers: [A136117](https://oeis.org/A136117)\n",
    "\n",
    "#### 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
}