Project Euler Problem 555
The McCarthy 91 function is defined as follows: We can generalize this definition by abstracting away the constants into
Solution
Answer: 208517717451208352
Let
$$M(n)=M_{m,k,s}(n)= \begin{cases} n-s & (n>m),\ M(M(n+k)) & (0\le n\le m). \end{cases}$$
We must compute
$$S(p,m)=\sum_{1\le s<k\le p} SF(m,k,s),$$
where $SF(m,k,s)$ is the sum of all fixed points of $M_{m,k,s}$.
We are asked for
$$S(10^6,10^6).$$
Step 1 — Understand the recursion
We analyze the function structurally.
For $n\le m$, if $n+k>m$, then the first recursive call immediately exits the recursion:
$$M(n+k)=n+k-s.$$
Define
$$d=k-s.$$
Then
$$M(n+k)=n+d.$$
So for $m-k<n\le m$,
$$M(n)=M(n+d).$$
Now consider what happens when repeatedly adding $d$.
Eventually we cross above $m$, and then the function returns by subtracting $s$.
The key observation is:
- If $d\nmid s$, the process cycles through residue classes and no fixed point exists.
- If $d\mid s$, then fixed points do exist.
Write
$$s=cd,\qquad k=(c+1)d.$$
Then the fixed points are exactly
$$m-s+1,; m-s+2,;\dots,; m-s+d.$$
So:
- number of fixed points = $d$,
- their sum is
$$SF(m,k,s) = \sum_{i=1}^{d}(m-s+i).$$
Compute this sum:
$$SF(m,k,s) = d(m-s)+\frac{d(d+1)}2.$$
Since $s=cd$,
$$SF = d(m-cd)+\frac{d(d+1)}2.$$
Step 2 — Re-parameterize the sum
Because
$$k=(c+1)d,$$
every valid pair $(k,s)$ with fixed points corresponds uniquely to:
- $d\ge1$,
- $c\ge1$,
- $(c+1)d\le p$.
Equivalently, define
$$t=c+1,$$
so $t\ge2$ and
$$k=td,\qquad s=(t-1)d.$$
Then
$$SF = d\left(m-(t-1)d\right)+\frac{d(d+1)}2.$$
Simplify:
$$SF = \frac d2\Bigl(2m-(2t-3)d+1\Bigr).$$
Therefore
$$S(p,m) = \sum_{d=1}^{p} \sum_{t=2}^{\lfloor p/d\rfloor} \frac d2\Bigl(2m-(2t-3)d+1\Bigr).$$
Step 3 — Closed-form inner summation
Let
$$N=\left\lfloor\frac pd\right\rfloor.$$
Then
$$S(p,m) = \sum_{d=1}^{p} \left[ \frac{d(N-1)(2m+1+3d)}2 - d^2\left(\frac{N(N+1)}2-1\right) \right].$$
This is only $O(p)$, fast enough for $p=10^6$.
Step 4 — Python implementation
def S(p, m):
total = 0
for d in range(1, p + 1):
N = p // d
cnt = N - 1
total += (
d * cnt * (2 * m + 1 + 3 * d) // 2
- d * d * (N * (N + 1) // 2 - 1)
)
return total
# Verify the examples
print(S(10, 10)) # 225
print(S(1000, 1000)) # 208724467
# Required answer
print(S(10**6, 10**6))
Step 5 — Code walkthrough
Outer loop
for d in range(1, p + 1):
Iterates over
$$d=k-s.$$
Number of admissible $t$
N = p // d
cnt = N - 1
Since
$$k=td\le p,$$
we have
$$2\le t\le N.$$
So there are $N-1$ values.
Add the closed-form contribution
total += (
d * cnt * (2 * m + 1 + 3 * d) // 2
- d * d * (N * (N + 1) // 2 - 1)
)
This is exactly the simplified inner summation derived above.
Step 6 — Verify against given examples
The program produces:
$$S(10,10)=225$$
and
$$S(1000,1000)=208724467,$$
matching the problem statement exactly.
Step 7 — Final computation
Evaluating the formula for
$$p=m=10^6$$
gives
$$208517717451208352.$$
The magnitude is reasonable:
- roughly quadratic/cubic growth in $p$,
- positive,
- consistent with the $1000$-case scale.
Answer: 208517717451208352