eulerbooks/notebooks/problem0074.ipynb

236 lines
7.4 KiB
Plaintext

{
"cells": [
{
"cell_type": "markdown",
"id": "c5bcd383",
"metadata": {},
"source": [
"# [Digit Factorial Chains](https://projecteuler.net/problem=74)\n",
"\n",
"We previously looked at the sum of digit factorials in [problem 34](https://projecteuler.net/problem=34). This time's a little different, since we're *repeatedly* summing digit factorials to find cycles. We can reuse our function implementation from that problem, though."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "9911ad22",
"metadata": {},
"outputs": [],
"source": [
"from functools import cache\n",
"\n",
"@cache\n",
"def sfd(n):\n",
" q = n // 10\n",
" if q == 0:\n",
" return factorial(n)\n",
" \n",
" return factorial(n % 10) + sfd(q)"
]
},
{
"cell_type": "markdown",
"id": "a2fa0cb6",
"metadata": {},
"source": [
"A natural way to find the chain lengths is with a memoized recursive function, since we know every starting value will eventually reach one of these cases."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "cbd97ddc",
"metadata": {},
"outputs": [],
"source": [
"@cache\n",
"def chain_length(n):\n",
" if n in {1, 2, 145, 40585}:\n",
" return 1\n",
" elif n in {169, 363601, 1454}:\n",
" return 3\n",
" elif n in {871, 45361, 872, 45362}:\n",
" return 2\n",
" \n",
" return 1 + chain_length(sfd(n))"
]
},
{
"cell_type": "markdown",
"id": "21d902d1",
"metadata": {},
"source": [
"Then we just iterate through every starting value below 1000000."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "084ed910",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"402"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"long_chains = []\n",
"for n in range(0, 1000000):\n",
" if chain_length(n) == 60:\n",
" long_chains.append(n)\n",
" \n",
"len(long_chains)"
]
},
{
"cell_type": "markdown",
"id": "b01d0dd6",
"metadata": {},
"source": [
"This is decently fast, but with some clever combinatorics, we can get the answer even quicker.\n",
"\n",
"## Multisets and multinomials\n",
"\n",
"Suppose a number $n$ is composed solely of digits 1-9 (we'll get to including 0 in a bit). Every permutation of those digits will have the same digit factorial sum as $n$, and therefore the same chain length. This means we only have to check each multiset of digits 1-9 once, and if the chain length is 60, we'll count every distinct permutation of those digits using the multinomial coefficient. (See [problem 92](https://projecteuler.net/problem=92) for a little more explanation of these concepts.)\n",
"\n",
"There are quirks with handling 0 here. Since $0! = 1$, it might seem that we can replace any instance of 1 in a number with 0 and get another number with the same chain length. However, this fails to handle the case that 1 is the first digit. When that happens, since the 0 is now a leading 0, we don't include it in the sum of digit factorials calculation, which changes the chain length. For example, the number 1218 has a digit factorial sum of 40324, as does 1208, but 0218 is just 218, and its digit factorial sum is 40323.\n",
"\n",
"For that reason, if we simply count the multinomial coefficient twice for every number with 1 as a digit, we'll overcount. Fortunately, we can correct for this using the [inclusion-exclusion principle](https://en.wikipedia.org/wiki/Inclusion%E2%80%93exclusion_principle).\n",
"\n",
"Counting with combinatorics reduces the search space significantly - from a million numbers to check to only thousands. We also end up finding that there aren't that many combinations of digits that result in a chain length of 60."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "6409c76a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{(1, 4, 7, 9), (2, 2, 3, 4, 7, 9)}"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from itertools import combinations_with_replacement\n",
"\n",
"cs = set()\n",
"for i in range(1, 7):\n",
" for digits in combinations_with_replacement(range(1, 10), i):\n",
" n = sum(10^k * d for (k, d) in enumerate(reversed(digits)))\n",
" if chain_length(n) == 60:\n",
" cs.add(digits)\n",
" \n",
"cs"
]
},
{
"cell_type": "markdown",
"id": "34e972e0",
"metadata": {},
"source": [
"As explained above, any permutation of these two multisets of digits will result in a number with a chain length of 60. Additionally, if we replace 1 in `(1, 4, 7, 9)` with 0, any permutation where 0 is not the most significant digit will also have a chain length of 60.\n",
"\n",
"You could just work this by hand at this point:\n",
"$$\\frac{4!}{(1!)(1!)(1!)(1!)} + 3 \\times \\frac{3!}{(1!)(1!)(1!)} + \\frac{6!}{(2!)(1!)(1!)(1!)(1!)(1!)} = 402$$\n",
"Or we could have the computer count for us."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "3141e8dc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"402"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from collections import Counter\n",
"from itertools import count\n",
"\n",
"def include_exclude(multiset):\n",
" total = 0\n",
" for i in count(0):\n",
" total += (-1)^i * multinomial(multiset.values())\n",
" multiset[1] -= 1\n",
" if multiset[1] == 0:\n",
" del multiset[1]\n",
" total += (-1)^(i+1) * multinomial(multiset.values())\n",
" break\n",
" \n",
" return total\n",
"\n",
"\n",
"def sfd_permutations(digits):\n",
" multiset = Counter(digits)\n",
" total = multinomial(multiset.values())\n",
"\n",
" if multiset[1] >= 1:\n",
" total += include_exclude(multiset)\n",
" \n",
" return total\n",
"\n",
"\n",
"sum(sfd_permutations(digits) for digits in cs)"
]
},
{
"cell_type": "markdown",
"id": "a0f338f0",
"metadata": {},
"source": [
"## Relevant sequences\n",
"* Numbers that eventually cycle when summing digit factorials: [A188284](https://oeis.org/A188284)\n",
"* Digit factorial chain lengths: [A303935](https://oeis.org/A303935)\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
}