Problem 120: Square Remainders

View on Project Euler

Project Euler Problem 120 Solution

EulerSolve provides an optimized solution for Project Euler Problem 120, Square Remainders, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer \(a\) with \(3 \le a \le 1000\), define \(r_n(a)\) to be the remainder when \((a-1)^n + (a+1)^n\) is divided by \(a^2\). The problem asks for the largest possible remainder \(r_{\max}(a)\) for each fixed \(a\), and then for the sum of these maxima over the whole range. The important point is that the exponent \(n\) is not fixed: we are free to choose the value of \(n\) that makes the remainder as large as possible. A brute-force search over many exponents would miss the main structure. Once the expression is reduced modulo \(a^2\), almost the entire binomial expansion disappears, and the problem collapses to a short parity argument together with a small modular-arithmetic case split on whether \(a\) is odd or even. Mathematical Approach The quantity to analyze is $$r_n(a) \equiv (a-1)^n + (a+1)^n \pmod{a^2}.$$ The whole solution comes from understanding which terms of the two powers can survive modulo \(a^2\). Binomial collapse modulo \(a^2\) Expand both powers around \(\pm 1\): $$ (a+1)^n = 1 + na + \binom{n}{2}a^2 + \cdots, $$ $$ (a-1)^n = (-1)^n + n(-1)^{n-1}a + \binom{n}{2}(-1)^{n-2}a^2 + \cdots. $$ Modulo \(a^2\), every term containing \(a^2\) or a higher power vanishes. Only the constant and linear terms matter, so $$ (a+1)^n \equiv 1 + na \pmod{a^2}, \qquad (a-1)^n \equiv (-1)^n + n(-1)^{n-1}a \pmod{a^2}....

Detailed mathematical approach

Problem Summary

Take a whole number \(a\) and an exponent \(n=1,2,3,\dots\), and let \(r\) be the remainder when \((a-1)^n+(a+1)^n\) is divided by \(a^2\). For a fixed \(a\) the remainder changes as \(n\) changes, and \(r_{\max}\) is the largest value it ever takes; each \(a\) may use its own best exponent. The statement's example is \(a=7\) and \(n=3\): \(6^3+8^3=728\) leaves the remainder \(42\) after division by \(49\), and no exponent does better for \(a=7\), so \(r_{\max}=42\) there. We must find the sum of \(r_{\max}\) over all \(a\) with \(3\le a\le1000\). [1]

The whole solution in one sentence: multiplying out the two powers shows that, after division by \(a^2\), the sum leaves the same remainder as \(2\) when \(n\) is even and as \(2na\) when \(n\) is odd, and the largest remainder this allows is \(a(a-1)\) for odd \(a\) and \(a(a-2)\) for even \(a\), which the program adds up.

Mathematical Approach

We start with small cases and a direct search, explain the pattern with the binomial theorem, prove the formula for \(r_{\max}\), and finally add the maxima up.

1. Small cases

Write \(r_n(a)\) for the remainder of \((a-1)^n+(a+1)^n\) after division by \(a^2\). For \(a=3\) the divisor is \(9\), and the sums \(2+4=6\), \(4+16=20=2\cdot9+2\) and \(8+64=72=8\cdot9\) for \(n=1,2,3\) leave the remainders \(6\), \(2\) and \(0\): the powers grow, but the remainders do not. For the first few \(a\) and \(n=1,2,3,\dots\) the remainders are:

  • \(a=3\), divisor \(9\): \(6,2,0,2,3,2\), then the same again; the largest is \(6\).
  • \(a=4\), divisor \(16\): \(8,2,8,2,\dots\); the largest is \(8\).
  • \(a=5\), divisor \(25\): \(10,2,5,2,0,2,20,2,15,2\), then the same again; the largest is \(20\), first reached at \(n=7\).
  • \(a=6\), divisor \(36\): \(12,2,0,2,24,2\), then the same again; the largest is \(24\).
  • \(a=7\), divisor \(49\): \(14,2,42,2,21,2,0\) for \(n=1,\dots,7\) and \(2,28,2,7,2,35,2\) for \(n=8,\dots,14\), then the same again; the largest is \(42\), as the statement says.
  • \(a=8\), divisor \(64\): \(16,2,48,2\), then the same again; the largest is \(48\).

Three things stand out. Every even \(n\) gives the remainder \(2\). Every odd \(n\) gives a multiple of \(a\), such as \(14=2\cdot7\), \(42=6\cdot7\) and \(21=3\cdot7\). And the remainders repeat, so the maximum over infinitely many exponents is really a maximum over finitely many. The maxima \(6,8,20,24,42,48\) are \(3\cdot2\), \(4\cdot2\), \(5\cdot4\), \(6\cdot4\), \(7\cdot6\) and \(8\cdot6\), which suggests the guess \(r_{\max}=a(a-1)\) for odd \(a\) and \(r_{\max}=a(a-2)\) for even \(a\). The next sections explain the three observations and prove the guess.

2. Remainders and congruences

Dividing \(728\) by \(49\) gives \(728=14\cdot49+42\): the quotient is \(14\) and the remainder is \(42\). For a divisor \(m\ge1\), every integer \(x\) can be written in exactly one way as \(x=qm+s\) with an integer \(q\) and a remainder \(s\) between \(0\) and \(m-1\); we write \(s=x\bmod m\). Two numbers that leave the same remainder are called congruent modulo \(m\), written \(x\equiv y\pmod m\), and this happens exactly when \(m\) divides the difference \(x-y\). For example \(728\equiv42\pmod{49}\), because \(728-42=686=14\cdot49\).

Congruences can be added and multiplied like equations: if \(x\equiv x'\) and \(y\equiv y'\) modulo \(m\), then \(x+y\equiv x'+y'\) and \(xy\equiv x'y'\), and therefore also \(x^k\equiv x'^k\) for every whole-number exponent \(k\ge0\). [2] The reason fits in one line: \(xy-x'y'=x\,(y-y')+(x-x')\,y'\) is a sum of two multiples of \(m\). So we may replace any number by its remainder at any moment and drop every multiple of \(m\).

3. A direct search, and why it may stop

A direct search computes the two powers one exponent at a time. It keeps \(x=(a-1)^n\bmod a^2\) and \(y=(a+1)^n\bmod a^2\), starts from \(x=y=1\) for \(n=0\), and for each new exponent multiplies \(x\) by \(a-1\) and \(y\) by \(a+1\), reducing modulo \(a^2\) after every multiplication; by Section 2 this gives the same remainders as the exact powers. This is the memory-efficient way of computing powers modulo a number, [3] and it matters: \(1001^{2000}\) has \(6001\) digits, while the reduced numbers stay below \(a^2\).

When may the search stop? The pair \((x,y)\) takes at most \(a^4\) different values, so by the pigeonhole principle, which says that more objects than boxes force two objects into one box, some pair appears twice. [4] In fact the first pair to come back is the start \((1,1)\), because every step can be undone. A common divisor of \(a\) and \(a\pm1\) divides their difference \(1\), so \(a\pm1\) shares no factor with \(a^2\) and has an inverse modulo \(a^2\), a partner whose product with it leaves the remainder \(1\). [5] For \(a=7\), \(8\cdot43=344=7\cdot49+1\), so multiplying by \(43\) undoes a multiplication by \(8\) modulo \(49\). If the first repetition were the pair of a step \(i\ge1\) coming back at step \(j\), undoing \(i\) steps would bring \((1,1)\) back already at step \(j-i\lt j\). So the search may stop at the first \(n\ge1\) with \(x=y=1\); from there on, everything repeats.

For \(a=7\) this happens at \(n=14\), and for \(a=8\) at \(n=8\). The first exponent at which the powers of one number return to \(1\) is called its multiplicative order; [6] Section 6 shows that the pair returns at \(n=a\) for even \(a\) and at \(n=2a\) for odd \(a\). The search over all \(3\le a\le1000\) then takes \(750496\) steps; it solves the problem on its own and checks everything that follows.

4. Half of the binomial expansion cancels

To see why the pattern holds we multiply the powers out. The binomial theorem says

$$(x+y)^n=\sum_{j=0}^{n}\binom nj x^{n-j}y^{j},$$

where the binomial coefficient \(\binom nj=\frac{n!}{j!\,(n-j)!}\), with \(t!=1\cdot2\cdots t\) and \(0!=1\), counts the ways to choose the \(j\) factors \(x+y\), out of \(n\), from which \(y\) is taken; these numbers form Pascal's triangle. [7] For \(n=3\), \(y=a\) and \(x=\pm1\):

$$\begin{aligned}(a+1)^3&=a^3+3a^2+3a+1,\\(a-1)^3&=a^3-3a^2+3a-1.\end{aligned}$$

When the two lines are added, \(3a^2\) and \(1\) cancel and the other terms double: \((a+1)^3+(a-1)^3=2a^3+6a\). For \(a=7\) this is \(686+42=728\), and since \(686=14\cdot49\), the remainder \(42\) of the statement appears without any division. In the same way \(n=2\) gives \(2a^2+2\), \(n=4\) gives \(2a^4+12a^2+2\), and \(n=5\) gives \(2a^5+20a^3+10a\).

In general the term with \(a^j\) has the sign \(+1\) in \((a+1)^n\) and \((-1)^{n-j}\) in \((a-1)^n\), so in the sum it gets the factor \(1+(-1)^{n-j}\), which is \(2\) when \(n-j\) is even and \(0\) when \(n-j\) is odd. Hence

$$(a+1)^n+(a-1)^n=2\sum_{j}\binom nj a^j,$$

where \(j\) runs over the numbers from \(0\) to \(n\) that have the same parity as \(n\): an even \(n\) keeps only the even powers of \(a\), and an odd \(n\) only the odd powers.

5. The remainder for every exponent

Modulo \(a^2\) every term with \(j\ge2\) contains the factor \(a^2\) and can be dropped. For even \(n\) only the term \(j=0\) is left, which is \(2\binom n0=2\); for odd \(n\) only the term \(j=1\), which is \(2\binom n1a=2na\). So the sum is congruent to \(2\) for even \(n\) and to \(2na\) for odd \(n\), modulo \(a^2\).

These are congruences, but the remainders follow at once. For even \(n\) the number \(2\) already lies between \(0\) and \(a^2-1\), so \(r_n(a)=2\). For odd \(n\), divide \(2n\) by \(a\), say \(2n=qa+s\) with \(0\le s\le a-1\). Then \(2na=q\,a^2+s\,a\) with \(0\le sa\lt a^2\), so for odd \(n\)

$$r_n(a)=a\cdot(2n\bmod a).$$

For example, \(a=7\) and \(n=5\) give \(2n=10=7+3\) and the prediction \(7\cdot3=21\); directly, \(6^5+8^5=7776+32768\)\({}=40544=827\cdot49+21\). This proves the first two observations of Section 1.

6. When the search returns to the start

The argument of Section 5, applied to one power at a time, gives

$$\begin{aligned}(a+1)^n&\equiv1+na,\\(a-1)^n&\equiv(-1)^n(1-na)\end{aligned}$$

modulo \(a^2\). The first congruence also has a proof by induction, which checks a statement for \(n=0\) and then shows that each case implies the next: \((a+1)^{n+1}\equiv(1+na)(a+1)\)\({}=1+(n+1)a+na^2\). For \(a=7\) and \(n=3\): \(8^3=512=10\cdot49+22\) and \(1+3\cdot7=22\), while \(6^3=216=4\cdot49+20\) and \(-(1-21)=20\).

Now we can see when the search of Section 3 returns to \((1,1)\). The number \(1+na\) leaves the remainder \(1\) exactly when \(a^2\) divides \(na\), that is, when \(a\) divides \(n\); so the order of \(a+1\) is \(a\). For \(a-1\) and odd \(n\) we would need \(na-1\equiv1\), so \(a^2\) would divide \(na-2\) and \(a\) would divide \(2\), impossible for \(a\ge3\); for even \(n\) we need, again, that \(a\) divides \(n\). So \(x\) and \(y\) are both \(1\) exactly when \(n\) is even and a multiple of \(a\), and the first such \(n\) is \(a\) for even \(a\) and \(2a\) for odd \(a\). This is the period of the pair \((x,y)\); the remainders themselves can repeat sooner: for \(a=8\) the pair needs \(8\) steps, but the list of Section 1 repeats after \(4\).

7. Which multiples of a can appear

By Section 5, \(r_{\max}\) is the larger of \(2\) and \(a\) times the largest value of \(s=2n\bmod a\) over odd \(n\). With \(n=2k+1\) we get \(2n=4k+2\): from one odd \(n\) to the next, \(s\) moves four places forward on a clock with \(a\) hours. Let \(d=\gcd(4,a)\) be the greatest common divisor of \(4\) and \(a\), the largest number that divides both: \(d=1\) for odd \(a\), \(d=2\) when \(a/2\) is odd, and \(d=4\) when \(4\) divides \(a\). Since \(d\) divides \(4\) and \(a\), every \(s=4k+2-qa\) leaves the same remainder as \(2\) after division by \(d\). The clock is back at its start after \(L\) steps exactly when \(a\) divides \(4L\), that is, when \(a/d\) divides \(L\), because dividing \(a\) and \(4\) by their greatest common divisor leaves no common factor. [8] So the first \(a/d\) values of \(s\) are all different, and as exactly \(a/d\) numbers from \(0\) to \(a-1\) leave the same remainder as \(2\) after division by \(d\), the odd exponents reach every one of them:

  • for odd \(a\), every \(s\) from \(0\) to \(a-1\), in the order \(2,6,3,0,4,1,5\) for \(a=7\);
  • when \(a/2\) is odd, every even \(s\), for example \(2,0,4\) for \(a=6\);
  • when \(4\) divides \(a\), only \(s=2,6,\dots,a-2\), for example \(2,6\) for \(a=8\) and \(2,6,10\) for \(a=12\).

8. The largest remainder

Now the guess of Section 1 follows. By Section 7 the largest value of \(s\) is \(a-1\) for odd \(a\) and \(a-2\) for even \(a\); for even \(a\), for example, the odd exponent \(n=a-1\) gives \(2n=2a-2=a+(a-2)\). The even exponents give only \(2\), less than the remainder \(2a\) of \(n=1\), although not every odd exponent does better: for \(a=3\), \(n=3\) gives \(0\). Hence \(r_{\max}=a(a-1)\) for every odd \(a\ge3\) and \(r_{\max}=a(a-2)\) for every even \(a\ge4\). In both cases the best \(s\) is the largest even number below \(a\), which is \(2\lfloor(a-1)/2\rfloor\), where \(\lfloor z\rfloor\) means \(z\) rounded down; so the two cases are one formula:

$$r_{\max}=2a\left\lfloor\frac{a-1}{2}\right\rfloor.$$

The On-Line Encyclopedia of Integer Sequences lists these maxima, \(6,8,20,24,42,48,72,80,110,120\) for \(a=3,\dots,12\) and so on, as sequence A159469, with the same formula and a link to this problem. [9]

Which exponent reaches the maximum first? For odd \(a\) we need \(2n\equiv a-1\), that is, \(n\equiv(a-1)/2\pmod a\), with \(n\) odd: this gives \(n=(a-1)/2\) when \(a\) leaves the remainder \(3\) after division by \(4\), as \(n=3\) for \(a=7\), and \(n=(a-1)/2+a=(3a-1)/2\) when it leaves the remainder \(1\), as \(n=7\) for \(a=5\). For even \(a\) we need \(n\equiv a/2-1\) modulo \(a/2\), with \(n\) odd: this gives \(n=a/2-1\) when \(4\) divides \(a\), as \(n=3\) for \(a=8\), and \(n=a-1\) otherwise, as \(n=5\) for \(a=6\).

9. A complete example

For \(a=6\), the odd exponents \(n=1,3,5\) give \(2n=2,6,10\), which leave the remainders \(s=2,0,4\) after division by \(6\), so Section 5 predicts the remainders \(12\), \(0\) and \(24\); the even exponents give \(2\). The maximum is \(24=6\cdot4\), at \(n=5\), and directly \(5^5+7^5=3125+16807=19932\)\({}=553\cdot36+24\). Adding the maxima one by one from \(a=3\) gives the running totals \(6,14,34,58,100\) up to \(a=7\), and then \(148,220,300\) up to \(a=10\).

10. Pitfalls

  • Small \(a\). For \(a=2\) the divisor is \(4\), the odd exponents give \(2\cdot(2n\bmod2)=0\) and the even ones give \(2\), so \(r_{\max}=2\) comes from the even exponents and the formula \(a(a-2)=0\) fails; the sum, which starts at \(a=3\), avoids this case. The exponent \(n=0\), if it were allowed, would only add the remainder \(2\) and change nothing for \(a\ge3\).
  • A search that stops too early. Trying only \(n\le a\) misses the maximum whenever \(a\) leaves the remainder \(1\) after division by \(4\), since the first maximising exponent \((3a-1)/2\) exceeds \(a\): for \(a=5\) such a search finds \(10\) instead of \(20\), and for \(3\le a\le1000\) it gets \(249\) of the \(998\) maxima wrong. Trying all \(n\le2a-1\) is always enough.
  • Floating-point powers. A double-precision number keeps only \(53\) binary digits, about \(16\) decimal digits, so above \(2^{53}\) not every whole number is stored exactly. [10] So remainders of large powers computed in floating point cannot be trusted: for \(a=1000\), the double-precision numbers nearest to \(999^6+1001^6\) and \(999^7+1001^7\) leave the remainders \(128\) and \(984640\) after division by \(10^6\), instead of \(2\) and \(14000\), and building the seventh powers by repeated floating-point multiplication gives yet another wrong value, \(246784\). Exact whole numbers, reduced after every multiplication, avoid this.

11. A larger divisor

The same algebra handles the divisor \(a^3\), but now only the terms with \(j\ge3\) drop out: for \(a=3\) and \(n=2\), \(2^2+4^2=20\) leaves the remainder \(2\) after division by \(9\), but \(20\) after division by \(27\). For odd \(n\) the identity of Section 4 has no term with \(a^2\), so the remainder is \(a\cdot(2n\bmod a^2)\). For even \(n\) the sum is congruent to \(2+n(n-1)a^2\), and its remainder \(2+a^2\cdot\bigl(n(n-1)\bmod a\bigr)\) is at most \(a^3-a^2+2\). As in Sections 7 and 8, \(2n\bmod a^2\) for odd \(n\) reaches \(a^2-1\) for odd \(a\) and \(a^2-2\) for even \(a\), so the largest remainder is \(a^3-a\) for odd \(a\) and \(a^3-2a\) for even \(a\); the even exponents fall short, because \(a^3-a^2+2\lt a^3-2a\) for \(a\ge3\). For \(a=3\), \(4\) and \(7\) the maxima are \(24\), \(56\) and \(336\), first reached at \(n=13\), \(7\) and \(73\), and a search over all \(n\le2a^2\), a full period, confirms the formula for \(3\le a\le80\).

12. Adding up the maxima

The requested number is \(S(A)=\sum_{a=3}^{A}r_{\max}\) for \(A=1000\). The program adds the \(998\) terms in a loop, but a closed formula exists as well. Since \(r_{\max}=a(a-1)\) for odd \(a\) and \(a(a-1)-a\) for even \(a\), \(S(A)\) is the sum of \(a(a-1)\) over \(3\le a\le A\) minus the sum of the even \(a\) from \(4\) to \(A\). Now \(a(a-1)=2\binom a2\), and the hockey-stick identity \(\sum_{a=2}^{A}\binom a2=\binom{A+1}3\) says that a diagonal of Pascal's triangle adds up to the number just below its end, as in \(1+3+6+10+15=35\). [11] Without the term \(a=2\), the first sum is therefore \((A+1)A(A-1)/3-2\). With \(m=\lfloor A/2\rfloor\), the even numbers \(4,6,\dots,2m\) are twice \(2,3,\dots,m\), so they add up to \(m(m+1)-2\). The two corrections \(-2\) cancel, and

$$\begin{aligned}S(A)&=\frac{(A+1)A(A-1)}{3}\\&\quad{}-m(m+1).\end{aligned}$$

For example, \(S(7)=112-12=100\) and \(S(10)=330-30=300\), as in Section 9, and the formula agrees with the loop for every \(A\) from \(3\) to \(10000\). In particular \(S(A)\) grows like \(A^3/3\).

How the Code Works

None of the programs computes a power; what is left is the loop of Section 12. The C++ program [12] keeps its settings in the structure Options: the limit max_a = 1000 and the switch run_checkpoints = true. parse_arguments accepts --skip-checkpoints, which turns the checks off, and --max-a= followed by digits, which parse_int_after_prefix converts one digit at a time into a new limit; any other argument prints "Unknown argument" and ends the program with exit code \(1\), and a limit below \(3\) ends it with the same code, silently.

r_max is Section 8 in two lines: (a & 1) == 0 tests the last binary digit of a, which is \(0\) exactly for even numbers, and the product a * (a - 2) or a * (a - 1) is formed after both factors are converted to the 64-bit type i64, that is, std::int64_t. solve adds r_max(a) for a from \(3\) to max_a in an i64 total. Unless the checks are switched off, main first calls run_checkpoints, which requires \(r_{\max}=42\) for \(a=7\), the example of the statement, and the sum \(100\) for \(3\le a\le7\) from Section 9 and otherwise returns false, so that main ends with exit code \(2\); then main prints the total for the chosen limit. Neither size limit is checked: the total is exact up to max_a \(=3024616\), beyond which it no longer fits into a signed 64-bit integer, and the digits after --max-a= are not checked against the range of int. Such an overflow of a signed integer is undefined behaviour in C++; [13] in a test, --max-a=4294967299 was read as \(3\).

The Python version is one function solve with the range fixed by range(3, 1001), which stops just before \(1001\), and the test a % 2 == 0, without options or checkpoints; its integers have unlimited precision. The Java version puts the same loop on one line of main, with an int loop variable: in (long)a*(a-2) the cast turns a into a long before the multiplication, so the product and total use 64-bit arithmetic. [13] Even 32-bit int would do here, since \(S(A)\lt2^{31}\) for every \(A\le1860\).

Complexity and Verification

Counting each arithmetic operation on machine-sized numbers as one step, the program does a fixed amount of work for each \(a\), a parity test, a multiplication and an addition, so its time grows in proportion to \(A\), with \(998\) terms for \(A=1000\), and its memory is a few numbers. The closed formula of Section 12 needs a fixed number of operations for any \(A\). The direct search of Section 3 makes \(a\) or \(2a\) steps for each \(a\), about \(3A^2/4\) in all, each with two multiplications and reductions of numbers below \(a^2(a+1)\), so its time grows like \(A^2\).

The C++ program checks itself with its two checkpoints, and all three programs print the same number. Separate computations, not part of the published solutions, confirmed the following. The remainder formula of Section 5 agrees with repeated multiplication and with independently computed modular powers for all \(3\le a\le1000\) and \(1\le n\le2a\), \(1000994\) pairs, and the congruences of Section 6 hold for \(2\le a\le300\) and \(0\le n\le1000\). For \(3\le a\le1000\) the search of Section 3 first returns at the exponent of Section 6, and for \(3\le a\le2000\) the values of \(s\) are those of Section 7. A full-period search gives the maxima and the first maximising exponents of Section 8 for every \(a\) from \(3\) to \(2000\). The number printed by the programs agrees with the full-period search over \(3\le a\le1000\) and with the closed formula; the C++ program prints \(6\), \(14\), \(100\) and \(300\) for max_a \(=3\), \(4\), \(7\) and \(10\), and a build with the undefined-behaviour sanitizer finds nothing on the default range but flags the oversized option.

Footnotes and References

The references below supply the background facts; the remainders, the formula for \(r_{\max}\), the divisor \(a^3\) and the sum formula are derived in this article.

  1. Project Euler 120 — Square Remainders. The official statement: r is the remainder when (a − 1)ⁿ + (a + 1)ⁿ is divided by a²; for a = 7 and n = 3, 6³ + 8³ = 728 ≡ 42 (mod 49), and as n varies, the largest remainder for a = 7 is r_max = 42. The task is to find the sum of r_max for 3 ≤ a ≤ 1000.
  2. Wikipedia — Modular arithmetic. Two integers are congruent modulo m when their difference is a multiple of m; congruence is compatible with addition, subtraction, multiplication and exponentiation.
  3. Wikipedia — Modular exponentiation. The memory-efficient method reduces modulo m after every multiplication, using (a · b) mod m = [(a mod m) · (b mod m)] mod m, which keeps the numbers small.
  4. Wikipedia — Pigeonhole principle. If n items are put into m containers with n > m, at least one container must contain more than one item.
  5. Wikipedia — Modular multiplicative inverse. An inverse of a modulo m is an integer x with ax ≡ 1 (mod m), and it exists if and only if gcd(a, m) = 1.
  6. Wikipedia — Multiplicative order. For a coprime to n, the multiplicative order of a modulo n is the smallest positive integer k with aᵏ ≡ 1 (mod n).
  7. Wikipedia — Binomial theorem. (x + y)ⁿ is the sum of C(n, k) xⁿ⁻ᵏ yᵏ for k = 0, …, n, where C(n, k) = n!/(k!(n − k)!) is an entry of Pascal's triangle.
  8. Wikipedia — Greatest common divisor. The greatest common divisor of two integers is the largest positive integer dividing both; if gcd(a, b) = d, then gcd(a/d, b/d) = 1, and if a divides b · c, then a/d divides c.
  9. OEIS A159469. "Maximum remainder when (k + 1)ⁿ + (k − 1)ⁿ is divided by k² for variable n and k > 2": 6, 8, 20, 24, 42, 48, … from k = 3, with the formula k² − 2k for even k and k² − k for odd k, and a link to Project Euler Problem 120.
  10. Wikipedia — Double-precision floating-point format. The 53-bit significand gives 15 to 17 significant decimal digits; integers up to 2⁵³ in size are exact, larger ones are rounded to multiples of 2, 4 and so on.
  11. Wikipedia — Hockey-stick identity. For integers n ≥ r ≥ 0, the sum of C(i, r) for i = r, …, n equals C(n + 1, r + 1); for r = 2 and n = 6 the page gives 1 + 3 + 6 + 10 + 15 = 35.
  12. C++ — Euler120.cpp. The linked immutable revision contains r_max, solve, run_checkpoints and the two options; it is the implementation record for this article.
  13. C++ — Arithmetic operators: when signed integer arithmetic overflows, the behaviour is undefined. Java Language Specification — Types, Values, and Variables: int and long are 32-bit and 64-bit two's-complement integers, and an integer operator other than a shift that has a long operand is carried out in 64-bit precision, after the other operand is widened to long. Python — numeric types: integers have unlimited precision.

Problem 120 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

using i64 = std::int64_t;

struct Options {
    int max_a = 1000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--max-a=", options.max_a)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.max_a >= 3;
}

i64 r_max(const int a) {
    if ((a & 1) == 0) {
        return static_cast<i64>(a) * static_cast<i64>(a - 2);
    }
    return static_cast<i64>(a) * static_cast<i64>(a - 1);
}

i64 solve(const int max_a) {
    i64 total = 0;
    for (int a = 3; a <= max_a; ++a) {
        total += r_max(a);
    }
    return total;
}

bool run_checkpoints() {
    if (r_max(7) != 42) {
        std::cerr << "Checkpoint failed for a=7" << '\n';
        return false;
    }
    if (solve(7) != 100) {
        std::cerr << "Checkpoint failed for range 3..7" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.max_a) << '\n';
    return 0;
}

Python

# Problem 120: Square remainders
# For 3 <= a <= 1000, find max r = ((a-1)^n + (a+1)^n) mod a^2.
# r_max(a) = a*(a-2) if a even, a*(a-1) if a odd.

def solve():
    total = 0
    for a in range(3, 1001):
        if a % 2 == 0:
            total += a * (a - 2)
        else:
            total += a * (a - 1)
    print(total)

solve()

Java

public class Euler120 {
    public static void main(String[] args) {
        long total=0;
        for(int a=3;a<=1000;a++)total+=(a%2==0)?(long)a*(a-2):(long)a*(a-1);
        System.out.println(total);
    }
}