Project Euler Problem 555

The McCarthy 91 function is defined as follows: We can generalize this definition by abstracting away the constants into

Project Euler Problem 555

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