Problem Summary

Define

$$T(a,n)=\left\lfloor\left(\lceil\sqrt{a}\rceil+\sqrt{a}\right)^n\right\rfloor,\qquad M=999999937.$$

The task is to evaluate

$$S(N)=\sum_{a=1}^{N} T(a,a^2)\pmod{M}$$

for \(N=5{,}000{,}000\). The exponent is already quadratic in \(a\), so direct floating-point evaluation is completely impractical; the solution instead converts the irrational power into an exact integer recurrence that can be computed modulo \(M\).

Mathematical Approach

The core observation is that the expression becomes easy once we pair it with its algebraic conjugate.

Step 1: Introduce the Conjugate Pair

For a fixed integer \(a\), let

$$m=\lceil\sqrt{a}\rceil,\qquad \alpha=m+\sqrt{a},\qquad \beta=m-\sqrt{a}.$$

Then \(\alpha\) and \(\beta\) are conjugates, and they satisfy

$$\alpha+\beta=2m,\qquad \alpha\beta=m^2-a.$$

These two symmetric quantities are integers, which is what makes the later recurrence integral.

Step 2: Replace the Floor by an Exact Integer Formula

If \(a\) is not a perfect square, then \(m-1<\sqrt{a}<m\), so

$$0<\beta=m-\sqrt{a}<1.$$

For every positive \(n\), this gives

$$0<\beta^n<1.$$

Now define

$$U_n=\alpha^n+\beta^n.$$

Since \(U_n-\alpha^n=\beta^n\) lies strictly between \(0\) and \(1\), we obtain

$$\left\lfloor\alpha^n\right\rfloor=U_n-1\qquad\text{for non-square }a.$$

If \(a\) is a perfect square, then \(m=\sqrt{a}\) and \(\beta=0\), so instead

$$\left\lfloor\alpha^n\right\rfloor=U_n=(2m)^n.$$

Therefore the entire problem reduces to computing \(U_n\) efficiently.

Step 3: Derive an Integer Recurrence for \(U_n\)

The numbers \(\alpha\) and \(\beta\) are the two roots of

$$x^2-2mx+(m^2-a)=0.$$

Therefore the sequence \(U_n=\alpha^n+\beta^n\) satisfies the standard second-order linear recurrence

$$U_0=2,\qquad U_1=2m,$$

$$U_n=2m\,U_{n-1}-(m^2-a)\,U_{n-2}\qquad(n\ge 2).$$

Because the coefficients and initial values are integers, every \(U_n\) is an integer. This is the exact quantity hidden behind the irrational-looking power.

Step 4: Compute \(U_n\) by Fast Exponentiation in a Quadratic Ring

Instead of iterating the recurrence all the way up to \(n=a^2\), the implementation exponentiates the base element \(m+\sqrt{a}\) directly. Represent

$$x+y\sqrt{a}$$

by the pair \((x,y)\). Then multiplication becomes

$$\left(x_1,y_1\right)\left(x_2,y_2\right)=\left(x_1x_2+a y_1y_2,\ x_1y_2+x_2y_1\right).$$

If

$$\alpha^n=X_n+Y_n\sqrt{a},$$

then its conjugate is

$$\beta^n=X_n-Y_n\sqrt{a},$$

so

$$U_n=\alpha^n+\beta^n=2X_n.$$

Binary exponentiation computes \(X_n\) in \(O(\log n)\) ring multiplications, which is essential because \(n=a^2\) can be enormous.

Step 5: Assemble the Final Summation

For each \(a\), the contribution is

$$T(a,a^2)=\begin{cases} U_{a^2}-1, & \text{if } a \text{ is not a perfect square},\\ U_{a^2}, & \text{if } a \text{ is a perfect square}. \end{cases}$$

Hence

$$S(N)=\sum_{a=1}^{N} T(a,a^2)\pmod{M}.$$

The only case distinction is whether \(a\) is square; the same ring exponentiation handles both branches.

Worked Example: \(a=5\)

Here \(m=\lceil\sqrt{5}\rceil=3\), so

$$\alpha=3+\sqrt{5},\qquad \beta=3-\sqrt{5},\qquad \alpha\beta=4.$$

The recurrence becomes

$$U_0=2,\qquad U_1=6,\qquad U_n=6U_{n-1}-4U_{n-2}.$$

Then

$$U_2=6\cdot 6-4\cdot 2=28,$$

so

$$T(5,2)=\left\lfloor(3+\sqrt{5})^2\right\rfloor=28-1=27.$$

Continuing,

$$U_3=144,\qquad U_4=752,\qquad U_5=3936,$$

and since \(0<\beta<1\),

$$T(5,5)=3936-1=3935.$$

This is exactly the kind of identity the implementations use, but computed modulo \(M\) and at much larger exponents.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They scan \(a\) over a contiguous range, keep track of the current value of \(\lceil\sqrt{a}\rceil\), and update it only when \(a\) crosses a perfect square. For each \(a\), they set \(n=a^2\), exponentiate \(m+\sqrt{a}\) by binary exponentiation in the pair representation above, double the real component to recover \(U_n\), subtract \(1\) when \(a\) is not square, and add the result into a running modular sum.

To speed up the full computation, the interval \([1,N]\) is partitioned into independent chunks. Each worker computes its own partial sum modulo \(M\), and the partial results are combined at the end. The mathematics is identical in all three languages; only the concurrency mechanism differs.

Complexity Analysis

For one value of \(a\), the cost is \(O(\log(a^2))=O(\log a)\) pair multiplications. Summing over all \(1\le a\le N\) gives

$$\sum_{a=1}^{N} O(\log a)=O(N\log N).$$

The working memory inside one worker is \(O(1)\). With parallel execution, total auxiliary memory is linear in the number of workers, but still constant per worker range.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=721
2. Lucas sequence: Wikipedia — Lucas sequence
3. Binary exponentiation: Wikipedia — Exponentiation by squaring
4. Integer square root: Wikipedia — Integer square root
5. Quadratic field: Wikipedia — Quadratic field