Problem 403: Lattice Points Enclosed by Parabola and Line

Project Euler Problem 403 Solution

EulerSolve provides an optimized solution for Project Euler Problem 403, Lattice Points Enclosed by Parabola and Line, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary This problem comes from a lattice-point configuration involving a parabola and a line. After the geometric reduction used by the implementation, the remaining task is to evaluate an exact arithmetic counting function $S(N)$. The Project Euler target is the residue of $S(10^{12})$ modulo $10^8$. The important point is that the implementation never enumerates lattice points one by one. Instead, it uses a closed formula made of a fixed polynomial part and a quotient-dependent summation whose values are grouped into floor-division blocks. Mathematical Approach Step 1: Polynomial Prefix Functions The reduced counting formula is built from three polynomial prefix functions: $$P_1(t)=\frac{(t^2-t+6)(t+1)}{6},$$ $$P_2(t)=\frac{(t^2-t+12)(t+1)(t+2)}{24},$$ $$P_3(t)=\frac{(t^2-t+20)(t+1)(t+2)(t+3)}{120}.$$ Expanding them shows that they are cubic, quartic, and quintic polynomials: $$P_1(t)=\frac{t^3+5t+6}{6},\qquad P_2(t)=\frac{t^4+2t^3+11t^2+34t+24}{24},$$ $$P_3(t)=\frac{t^5+5t^4+25t^3+115t^2+214t+120}{120}.$$ What matters computationally is their telescoping behavior: $$P_1(t)-P_1(t-1)=\frac{t(t-1)}{2}+1,$$ $$P_2(t)-P_2(t-1)=P_1(t),\qquad P_3(t)-P_3(t-1)=P_2(t).$$ So $P_2$ is a prefix sum of $P_1$, and $P_3$ is a prefix sum of $P_2$. This is exactly why long ranges can later be collapsed to differences of polynomial values....

Detailed mathematical approach

Problem Summary

For integers $a$ and $b$, let $D(a,b)$ be the domain enclosed by the parabola $y=x^2$ and the line $y=ax+b$,

$$D(a,b)=\{(x,y)\mid x^2\le y\le ax+b\},$$

and let $L(a,b)$ be the number of lattice points it contains, that is, points with integer coordinates, boundary included. The statement gives $L(1,2)=8$ and $L(2,-1)=1$. For a bound $N$, the quantity $S(N)$ adds up $L(a,b)$ over every pair $(a,b)$ with $|a|\le N$ and $|b|\le N$ whose domain has a rational area, and we must find $S(10^{12})$ modulo $10^8$. Two checkpoints are supplied: $S(5)=344$ and $S(100)=26709528$. ^[1]

The whole solution in one sentence: leaving aside the empty domains, which have area $0$ and no lattice points, the area is rational exactly when the line meets the parabola at two integer abscissas $p\le q$, possibly equal; the lattice count then depends only on $q-p$ through a single cubic polynomial, and summing that polynomial over all admissible pairs $(p,q)$ telescopes into $O(\sqrt N)$ blocks on which $\lfloor N/m\rfloor$ is constant.

Mathematical Approach

1. When the area is rational

The line meets the parabola where $x^2=ax+b$, that is, at the roots of $x^2-ax-b$. Write $\Delta=a^2+4b$; the roots are $\bigl(a\pm\sqrt\Delta\bigr)/2$.

If $\Delta\lt0$ the line stays strictly below the parabola, $D(a,b)$ is empty, its area is $0$ and $L(a,b)=0$; such pairs satisfy the rationality condition but add nothing to $S(N)$. If $\Delta\ge0$, call the roots $p\le q$. Then

$$ax+b-x^2=(x-p)(q-x),$$

which is non-negative exactly for $p\le x\le q$, so the domain lies between the two curves over that interval and its area is

$$\begin{aligned}\int_p^q(x-p)(q-x)\,dx&=\frac{(q-p)^3}{6}\\&=\frac{\Delta^{3/2}}{6}.\end{aligned}$$

This matches the classical fact that a parabolic segment fills two thirds of the parallelogram built on its chord and the parallel tangent. ^[2] That parallelogram has area $(q-p)^3/4$, since the tangent parallel to the chord touches at $x=(p+q)/2$, where the vertical gap is $(q-p)^2/4$.

The area is therefore $\Delta\sqrt\Delta/6$, and it is rational precisely when $\sqrt\Delta$ is. But $\sqrt\Delta$ is a root of the monic integer polynomial $t^2-\Delta$, and every rational root of a monic integer polynomial is an integer. ^[3] So the area is rational if and only if $\Delta$ is a perfect square, and only those pairs matter.

2. From the pair $(a,b)$ to two integer abscissas

Let $\Delta=d^2$ with $d\ge0$. From $4b=d^2-a^2$ we get $d^2\equiv a^2\pmod4$, so $a$ and $d$ have the same parity: if $a$ is even then $d^2\equiv0\pmod4$ forces $d$ even, and if $a$ is odd then $d^2\equiv1\pmod4$ forces $d$ odd. Hence

$$p=\frac{a-d}{2},\qquad q=\frac{a+d}{2}$$

are integers, and $a=p+q$, $b=-pq$, $d=q-p$. Conversely, every integer pair $p\le q$ gives a pair $(a,b)=(p+q,-pq)$ whose discriminant $(q-p)^2$ is a square. The two descriptions are therefore in bijection, and the box $|a|\le N$, $|b|\le N$ turns into

$$|p+q|\le N,\qquad |pq|\le N.$$

3. How many lattice points the domain holds

Fix such a pair and put $n=q-p\ge0$. For an integer $x$, the admissible values of $y$ are the integers of the interval from $x^2$ to $ax+b$, and both ends are integers, so there are exactly $ax+b-x^2+1=(x-p)(q-x)+1$ of them when that count is positive, and none otherwise. Summing over $x=p,\dots,q$ and substituting $t=x-p$,

$$L(a,b)=(n+1)+\sum_{t=0}^{n}t\,(n-t).$$

The last sum equals $\binom{n+1}{3}$. To see it, count the triples $i\lt t\lt j$ drawn from $\{0,1,\dots,n\}$ by first choosing the middle element $t$, which leaves $t$ choices for $i$ and $n-t$ for $j$. So the lattice count depends on $n$ alone:

$$F_1(n)=(n+1)+\binom{n+1}{3}=\frac{(n+1)(n^2-n+6)}{6},$$

equivalently $F_1(n)=(n^3+5n+6)/6$. Its first values, for $n=0,\dots,7$, are $1,2,4,8,15,26,42,64$, the cake numbers. ^[4] Both examples of the statement check out: $a=1$, $b=2$ give $\Delta=9$ and $(p,q)=(-1,2)$, so $L=F_1(3)=8$; and $a=2$, $b=-1$ give $\Delta=0$ and $(p,q)=(1,1)$, so $L=F_1(0)=1$, the single point $(1,1)$.

4. The smallest case in full

Take $N=1$. Of the nine pairs with $|a|,|b|\le1$, three have $\Delta\lt0$ and an empty domain, two have $\Delta=5$ and an irrational area, and four remain:

$(a,b)=(-1,0)$: $(p,q)=(-1,0)$, $n=1$, $L=2$;
$(a,b)=(0,0)$: $(p,q)=(0,0)$, $n=0$, $L=1$;
$(a,b)=(0,1)$: $(p,q)=(-1,1)$, $n=2$, $L=4$;
$(a,b)=(1,0)$: $(p,q)=(0,1)$, $n=1$, $L=2$.

Hence $S(1)=2+1+4+2=9$. For instance $D(0,1)$ is the sliver under the line $y=1$, with the four lattice points $(-1,1)$, $(0,0)$, $(0,1)$ and $(1,1)$.

5. The sum as a region of the plane

Summing $F_1(q-p)$ over the pairs of Section 2 becomes symmetric after the substitution $u=-p$, $v=q$. Then $p\le q$ reads $u+v\ge0$, while $a=v-u$, $b=uv$ and $q-p=u+v$, so

$$S(N)=\sum F_1(u+v),$$

the sum running over all integer pairs $(u,v)$ that satisfy

$$u+v\ge0,\quad |u-v|\le N,\quad |uv|\le N.$$

The three conditions are exactly $p\le q$, $|a|\le N$ and $|b|\le N$. The region they describe lies above the anti-diagonal $u+v=0$, inside the diagonal strip $|u-v|\le N$ and between the branches of the hyperbolas $uv=N$ and $uv=-N$ ^[6]; it is symmetric in $u$ and $v$. Two features make it tractable: the summand depends on $u+v$ only, and for each fixed $u$ the admissible $v$ form one unbroken interval. A vertical line of the region therefore contributes a run of consecutive values of $F_1$ — what prefix sums are for.

6. Three prefix polynomials

For $n\ge0$ put $F_2(n)=\sum_{j=0}^{n}F_1(j)$ and $F_3(n)=\sum_{j=0}^{n}F_2(j)$ — running totals, so $F_2(2)=1+2+4=7$ and $F_3(2)=1+3+7=11$. Wherever a negative argument appears below it means the value of the closed form, not an empty sum. Writing $F_1(n)=\binom{n+1}{3}+\binom{n+1}{1}$ and applying the hockey-stick identity $\sum_{i=r}^{m}\binom{i}{r}=\binom{m+1}{r+1}$ ^[5] to each of the two terms gives closed forms:

$$\begin{aligned}F_2(n)&=\binom{n+2}{4}+\binom{n+2}{2}\\&=\frac{(n^2-n+12)(n+1)(n+2)}{24},\end{aligned}$$

$$\begin{aligned}F_3(n)&=\binom{n+3}{5}+\binom{n+3}{3}\\&=\frac{(n^2-n+20)(n+1)(n+2)(n+3)}{120}.\end{aligned}$$

Thus $F_2$ starts $1,3,7,15,30,56$ and $F_3$ starts $1,4,11,26,56,112$. Two boundary properties matter later. A binomial coefficient $\binom{M}{r}$ vanishes for $0\le M\lt r$, so

$$F_2(-1)=F_2(-2)=0,\quad F_3(-1)=F_3(-2)=F_3(-3)=0,$$

as the explicit factors $(n+1)(n+2)$ and $(n+1)(n+2)(n+3)$ also show. These zeros let the edge cases below be written with no special case. One step further the polynomials stop vanishing — $F_3(-4)=-2$ — so we must check that no argument drops below $-3$.

7. One vertical line at a time

Let $T_N(u)$ be the contribution of the vertical line at $u$, and write $s=u+v$ for the argument of $F_1$. For $u\ne0$ the condition $|uv|\le N$ says $|v|\le\lfloor N/|u|\rfloor$, the condition $u+v\ge0$ says $v\ge-u$, and $|u-v|\le N$ says $u-N\le v\le u+N$. Assume $N\ge2$, and in the last two items write $c=\lfloor N/m\rfloor$.

$u=0$: the hyperbola condition is vacuous, and $0\le v\le N$. So $0\le s\le N$ and $T_N(0)=F_2(N)$.
$u=1$: here $|v|\le N$, $v\ge-1$ and $v\ge1-N$; the bound $v\ge-1$ wins because $N\ge2$. So $0\le s\le N+1$ and $T_N(1)=F_2(N+1)$.
$u=-1$: here $|v|\le N$, $v\ge1$, and $|{-1}-v|\le N$ gives $v\le N-1$. So $0\le s\le N-2$ and $T_N(-1)=F_2(N-2)$.
$u=N$: now $\lfloor N/N\rfloor=1$, so $|v|\le1$, while $|N-v|\le N$ forces $v\ge0$. Only $v=0$ and $v=1$ survive, so $T_N(N)=F_1(N)+F_1(N+1)$.
$u=-N$: the bounds $|v|\le1$ and $v\ge N\ge2$ are incompatible, so $T_N(-N)=0$. And if $|u|\gt N$ then $|uv|\le N$ forces $v=0$, which violates $|u-v|\le N$: nothing lies outside $|u|\le N$.
$u=m$ with $2\le m\le N-1$: the upper bound is $v\le c$, and the lower bound is the largest of $-c$, $-m$ and $m-N$. If $m\le c$, then $-m\ge-c$, and $2m\le m^2\le mc\le N$ gives $m-N\le-m$, so the lower bound is $-m$: thus $0\le s\le m+c$ and $T_N(m)=F_2(m+c)$. If $m\gt c$, then $-c\gt-m$, and $m+c\le N$ holds — for $c=1$ because $m\le N-1$, and for $c\ge2$ because $(m-1)(c-1)\ge1$ gives $m+c\le mc\le N$ — so the lower bound is $-c$: thus $m-c\le s\le m+c$, and the line contributes $F_2(m+c)-F_2(m-c-1)$.
$u=-m$ with $2\le m\le N-1$: now $v\ge m$ and $|v|\le c$, so the line is empty unless $m\le c$; and then $m+c\le mc\le N$ shows that $|{-m}-v|\le N$ is implied. Thus $0\le s\le c-m$ and $T_N(-m)=F_2(c-m)$ when $m\le c$, and $T_N(-m)=0$ otherwise.

8. The master formula

Adding the lines of Section 7, with $c=\lfloor N/m\rfloor$ inside the sum,

$$\begin{aligned}S(N)&=F_2(N)+F_2(N+1)+F_2(N-2)\\&\quad+F_1(N)+F_1(N+1)+\sum_{m=2}^{N-1}H_N(m),\end{aligned}$$

where the two cases of the last two items combine into

$$H_N(m)=\begin{cases}F_2(m+c)+F_2(c-m),&m\le c,\\F_2(m+c)-F_2(m-c-1),&m\gt c.\end{cases}$$

The first branch needs no emptiness test: for $m\le c$ the term $F_2(c-m)$ is the contribution of the line $u=-m$, and the second branch drops it because that line is then empty. The formula holds for every $N\ge2$. It fails at $N=1$, where the lines $u=1$ and $u=N$ are one and the same line, and where the derivation of $T_N(1)$ used $N\ge2$ as well; the verification section works the resulting overcount out in full.

9. Worked example: $S(5)=344$

Four lines sit in the fixed part: $T_5(0)=F_2(5)=56$, $T_5(1)=F_2(6)=98$, $T_5(-1)=F_2(3)=15$ and $T_5(5)=F_1(5)+F_1(6)=26+42=68$, together $237$. The sum over $m$ has three terms:

$m=2$, $c=2$, so $m\le c$: $F_2(4)+F_2(0)=30+1=31$;
$m=3$, $c=1$, so $m\gt c$: $F_2(4)-F_2(1)=30-3=27$;
$m=4$, $c=1$, again $m\gt c$: $F_2(5)-F_2(2)=56-7=49$.

Hence $S(5)=237+31+27+49=344$, as the statement says. Line by line instead, the same computation gives the contributions $1,15,56,98,30,27,49,68$ for $u=-2,-1,0,1,2,3,4,5$, and every other line is empty.

10. Collapsing the quotient into blocks

The remaining sum still has $N-2$ terms, but it depends on $m$ only through $m$ and $c=\lfloor N/m\rfloor$, and $c$ is constant on long blocks: if a block starts at $m=x$ and $c=\lfloor N/x\rfloor$, then every $m$ with $x\le m\le\lfloor N/c\rfloor$ has that same quotient. The blocks are few: put $r=\lfloor\sqrt N\rfloor$; there are at most $r$ indices with $m\le r$, while for $m\gt r$ we have $N/m\lt\sqrt N$, so $\lfloor N/m\rfloor$ is one of at most $r$ values and that range holds at most $r$ blocks. Hence at most $2r$ blocks in all, the square-root scale familiar from divisor sums. ^[6]

Since $F_3$ is the prefix sum of $F_2$, a whole block collapses. The master sum stops at $m=N-1$, so a block starting at $x$ ends at $e=\min(\lfloor N/c\rfloor,N-1)$; the cap bites only in the last block, where $c=1$. Without it the term $m=N$ would slip in: at $N=5$ it would add $F_2(6)-F_2(3)=83$ and give $427$ instead of $344$. With that $e$,

$$\sum_{m=x}^{e}F_2(m+c)=F_3(e+c)-F_3(x+c-1),$$

$$\sum_{m=x}^{e}F_2(c-m)=F_3(c-x)-F_3(c-e-1),$$

$$\sum_{m=x}^{e}F_2(m-c-1)=F_3(e-c-1)-F_3(x-c-2).$$

Only one branch of $H_N$ is needed per block, because no block straddles the diagonal $m=c$: if $x\le c$ then $e\le\lfloor N/c\rfloor\le\lfloor N/x\rfloor=c$, so the whole block satisfies $m\le c$; and if $x\gt c$ then so does every $m\ge x$. The smallest arguments reached are $c-e-1\ge-1$ in the first case, because $e\le c$, and $x-c-2\ge-1$ in the second, because $x\ge c+1$. As $F_3(-1)=0$, the collapsed formulas need no special case and never reach the region where $F_3$ stops vanishing.

Take $N=5$ once more as a check. The indices $m=3,4$ form one block, with $c=1$, $x=3$ and $e=4$. The collapsed formulas give $F_3(5)-F_3(3)=86$ for the first sum and $F_3(2)-F_3(0)=10$ for the correction, and $86-10=76$ is exactly the $27+49$ of Section 9. The whole computation is now a walk over $O(\sqrt N)$ blocks.

How the Code Works

The C++ program ^[9] is a direct transcription of Sections 8 and 10. F1, F2 and F3 evaluate the three closed forms of Sections 3 and 6 on boost::multiprecision::cpp_int, an arbitrary-precision signed integer type ^[7]; the divisions by $6$, $24$ and $120$ are exact, because each function is an integer-valued prefix sum.

S_exact first adds the four fixed lines as F2(N) + F2(N+1) + F2(N-2) + F1(N) + F1(N+1), then walks the blocks. It keeps the block start in x, computes c = n / x and next = std::min(n / c + 1, n), so the block runs from $m=x$ to $m=$ next - 1. Capping next at n is what stops the last block at $m=N-1$: the line $u=N$ already sits in the fixed part. For each block the program adds F3(next - 1 + c) - F3(x - 1 + c), the first collapsed sum of Section 10; then, if x > c, it subtracts F3(next - c - 2) - F3(x - c - 2), and otherwise it adds F3(c - x) - F3(c - next). Those are precisely the two branches of $H_N$, with next - 1 playing the role of $e$. The indices x, c and next are i64, an alias for long long; at the target bound their values stay near $N$, so only the polynomial values and the running total need wide arithmetic.

run_checkpoints verifies $S(5)=344$ and $S(100)=26709528$ before the real computation and leaves with exit code $2$ if either fails; --skip-checkpoints switches the check off. parse_arguments also accepts --n=, read by parse_i64_after_prefix with std::stoll inside a try block; it rejects unknown options, failed conversions and a bound below $1$. main reduces the exact sum modulo $10^8$ and adds $10^8$ back should the result be negative; that guard is defensive only, since $S(N)$ is positive and each block adds a larger prefix difference before subtracting a smaller one.

The Python file defines only F1, F2, F3 and solve: the last sets $n=10^{12}$, adds the fixed terms, runs the same blocks in a while x < n loop with // for the exact divisions, reduces and returns a string for the entry point to print. It needs no special type, because Python integers have unlimited precision. ^[8] The Java file keeps the polynomial values and the total in BigInteger ^[8] but the block indices in long, with Math.min for the block end; solveExact(long n) takes the bound as a parameter, and solve supplies $n=10^{12}$, reduces and tests the sign with compareTo. Neither port carries the checkpoint routine or the option parser, and neither can be given another bound without an edit.

One design decision deserves a note: all three evaluate the polynomials exactly and reduce only at the very end. Reducing a numerator modulo $10^8$ before its division would not preserve the quotient, and $120$ has no inverse modulo $10^8$ either, since $\gcd(120,10^8)=40$. What does work is to reduce a numerator $P$ modulo $120\cdot10^8$ before dividing: the remainder is still a multiple of the denominator, which divides $120$, and the quotient changes only by a multiple of $10^8$. We checked that this variant reproduces the same residue; the published programs simply keep everything exact.

Complexity and Verification

The loop runs once per distinct value of $\lfloor N/m\rfloor$, so it performs $O(\sqrt N)$ iterations; for $N=10^{12}$ we measured $1\,999\,998$ of them. Each iteration evaluates $F_3$ four times. The arguments never exceed $N$, so the values stay below $F_3(N+1)\approx8.3\cdot10^{57}$, which is $193$ bits, so a block costs a fixed number of multiplications of numbers of $O(\log N)$ bits. The whole computation is therefore $O(\sqrt N)$ big-integer operations, or $O(\sqrt N\log^2N)$ bit operations with schoolbook multiplication. Nothing is stored beyond a few temporaries, so the memory is $O(\log N)$ bits, a handful of machine words. The exact value of $S(10^{12})$ has $48$ decimal digits, or $158$ bits, so the accumulator cannot be a machine integer; the final reduction modulo $10^8$ is there because the problem asks for that residue, not because the exact sum had to be discarded. One run on one machine took $0.75$ s for the C++ version, $0.91$ s for the Java version and $1.66$ s for the Python version, all three printing the same residue; wall-clock time depends on the machine, the compiler and the runtime, so this is one measurement, not a benchmark.

Beyond the two checkpoints built into the C++ program, the following were verified independently for this article.

A direct enumeration over all $|a|,|b|\le N$, testing $a^2+4b$ for squareness and counting the lattice points one abscissa at a time, reproduces the block formula for every $2\le N\le400$ and for $N=1000$. It also gives $S(5)=344$, $S(100)=26709528$, $S(10)=3560$, $S(20)=47021$ and $S(1000)=263967605900$.
$L(a,b)=F_1(d)$ for every pair with $|a|,|b|\le60$ and $a^2+4b=d^2$, the points again being counted abscissa by abscissa; in particular $L(1,2)=8$ and $L(2,-1)=1$.
The line-by-line sum of Section 7, which costs $O(N)$ and uses no blocks, agrees with the $O(\sqrt N)$ block sum at $N=10^4$, $10^5$, $10^6$ and $10^7$.
The identities $F_2(n)-F_2(n-1)=F_1(n)$ and $F_3(n)-F_3(n-1)=F_2(n)$ hold for every integer $-10\le n\le300$, the three binomial forms for every $0\le n\le300$, and the numerators are divisible by $6$, $24$ and $120$ for every $-200\le n\le500$.
Recording every argument the block loop hands to $F_3$ shows that the smallest one is exactly $-1$, where $F_3$ vanishes, for every $3\le N\le3000$ and, from an index-only traversal that evaluates no polynomial, for $N=10^{12}$ as well, whose largest argument is $10^{12}$. At $N=2$ the loop body never runs, so there is nothing to record. No block contains indices from both branches of $H_N$.

Three limitations are worth recording as they stand. First, the C++ parser rejects only $n\lt1$, so --n=1 is accepted and the program returns $16$ instead of $S(1)=9$. The five fixed terms give $3+7+0+2+4=16$, while the true lines are $T_1(0)=3$ and $T_1(1)=6$: the line $u=1$ is counted a second time by $F_1(N)+F_1(N+1)=6$, and $F_2(2)=7$ exceeds $T_1(1)$ by the $F_1(0)=1$ of the pair $(u,v)=(1,-1)$, which $|u-v|\le1$ forbids. Second, parse_i64_after_prefix never asks std::stoll where the conversion stopped, so --n=5junk is taken as $5$ and --n=2.9 as $2$. Third, at the top of the accepted range, $n=2^{63}-1$, the last block has $c=1$ and n / c + 1 overflows before std::min can cap it. Every $N\ge2$ we tested, the target included, is correct.

Footnotes and References

The references below supply the background facts used above. The reduction to integer abscissas, the cubic lattice count, the line-by-line contributions and the collapsed block sums are derived explicitly in this article.

Project Euler 403 — Lattice Points Enclosed by Parabola and Line. The official statement defines D(a, b) as the set of points with x² ≤ y ≤ ax + b, defines L(a, b) as the number of lattice points it contains, gives L(1, 2) = 8 and L(2, −1) = 1, defines S(N) as the sum of L(a, b) over the pairs with |a|, |b| ≤ N whose domain has rational area, states S(5) = 344 and S(100) = 26709528, and asks for S(10¹²) modulo 10⁸.
Wikipedia — Parabola. Its section on the area enclosed between a parabola and a chord states that this area is two thirds of the parallelogram having the chord as one side and the tangent parallel to the chord as the opposite side, and gives A = (2/3)bh when the chord is perpendicular to the axis. It records that Archimedes proved the equivalent statement with triangles. Section 1 uses it as a check on the integral.
Wikipedia — Rational root theorem. A rational root p/q in lowest terms of a polynomial with integer coefficients has p dividing the constant term and q dividing the leading coefficient. The article names the monic case the integral root theorem: there q divides 1, so every rational root is an integer. Section 1 applies this to t² − Δ to conclude that √Δ is rational only when Δ is a perfect square.
OEIS A000125 — Cake numbers. The entry gives a(n) = C(n + 1, 3) + n + 1 = (n + 1)(n² − n + 6)/6 = (n³ + 5n + 6)/6, with the values 1, 2, 4, 8, 15, 26, 42, 64, 93, 130, and reads a(n) as the largest number of pieces that n planar cuts can produce in a cake. This is the function F₁ of Section 3, and the F1 of the three programs.
Wikipedia — Hockey-stick identity. The sum of C(i, r) for i from r to n equals C(n + 1, r + 1). It is proved there by induction with Pascal's rule, by telescoping, and by counting the (r + 1)-element subsets of {1, …, n + 1} by their largest element. Section 6 applies it twice, to turn F₁ into F₂ and F₂ into F₃.
Wikipedia — Divisor summatory function. The article counts the lattice points under the hyperbola jk = x and gives D(x) = 2·Σ⌊x/k⌋ over k ≤ u, minus u², with u = ⌊√x⌋, noting that this allows the sum to be computed in O(√x) time. It is the background for the hyperbolic boundary of Section 5 and for the square-root scale of the work in Section 10; the equal-quotient blocks and their collapsed sums are derived in Section 10 itself.
Boost.Multiprecision — cpp_int. The documentation describes cpp_int as an arbitrary-precision signed integer whose MinBits parameter fixes how many bits live inside the object before dynamic allocation is needed. The C++ solution holds every polynomial value and the running total in this type.
Python — numeric types documents that Python integers have unlimited precision and that // is the floored quotient, and Java — BigInteger documents immutable arbitrary-precision integers together with the add, subtract, multiply, divide, remainder and compareTo methods the Java port uses. Together they explain why neither port needs an external big-integer library.
C++ — Euler403.cpp. The linked immutable C++ revision contains F1, F2, F3, S_exact, run_checkpoints and the argument parser described above. It is the implementation record for this article; the cross-checks listed in the previous section are additional validation, not tests packaged in that file.

Problem 403 source code

C++

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

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using i64 = long long;
using boost::multiprecision::cpp_int;

struct Options {
    i64 n = 1000000000000LL;
    bool run_checkpoints = true;
};

bool parse_i64_after_prefix(const std::string& arg, const std::string& prefix, i64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoll(tail);
    } catch (...) {
        return false;
    }
    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_i64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

cpp_int F1(const cpp_int& n) {
    return (n * n - n + 6) * (n + 1) / 6;
}

cpp_int F2(const cpp_int& n) {
    return (n * n - n + 12) * (n + 1) * (n + 2) / 24;
}

cpp_int F3(const cpp_int& n) {
    return (n * n - n + 20) * (n + 1) * (n + 2) * (n + 3) / 120;
}

cpp_int S_exact(i64 n) {
    const cpp_int N = n;
    cpp_int ans = 0;

    // Closed-form decomposition by discriminant class and floor-division blocks.
    ans += F2(N) + F2(N + 1) + F2(N - 2) + F1(N) + F1(N + 1);

    for (i64 x = 2; x < n;) {
        const i64 c = n / x;
        const i64 next = std::min(n / c + 1, n);

        ans += F3(cpp_int(next - 1 + c)) - F3(cpp_int(x - 1 + c));
        if (x > c) {
            ans -= F3(cpp_int(next - c - 2)) - F3(cpp_int(x - c - 2));
        } else {
            ans += F3(cpp_int(c - x)) - F3(cpp_int(c - next));
        }

        x = next;
    }

    return ans;
}

bool run_checkpoints() {
    if (S_exact(5) != 344) {
        std::cerr << "Checkpoint failed: S(5)\n";
        return false;
    }
    if (S_exact(100) != 26709528) {
        std::cerr << "Checkpoint failed: S(100)\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;
    }

    cpp_int ans = S_exact(options.n) % cpp_int(100000000);
    if (ans < 0) {
        ans += cpp_int(100000000);
    }
    std::cout << ans << '\n';
    return 0;
}

Python

def F1(n):
    return (n * n - n + 6) * (n + 1) // 6

def F2(n):
    return (n * n - n + 12) * (n + 1) * (n + 2) // 24

def F3(n):
    return (n * n - n + 20) * (n + 1) * (n + 2) * (n + 3) // 120

def solve():
    n = 1000000000000
    ans = F2(n) + F2(n + 1) + F2(n - 2) + F1(n) + F1(n + 1)
    
    x = 2
    while x < n:
        c = n // x
        next_x = min(n // c + 1, n)
        ans += F3(next_x - 1 + c) - F3(x - 1 + c)
        if x > c:
            ans -= F3(next_x - c - 2) - F3(x - c - 2)
        else:
            ans += F3(c - x) - F3(c - next_x)
        x = next_x
        
    ans %= 100000000
    if ans < 0:
        ans += 100000000
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;

public class Euler403 {
    private static BigInteger F1(BigInteger n) {
        BigInteger term1 = n.multiply(n).subtract(n).add(BigInteger.valueOf(6));
        BigInteger term2 = n.add(BigInteger.ONE);
        return term1.multiply(term2).divide(BigInteger.valueOf(6));
    }

    private static BigInteger F2(BigInteger n) {
        BigInteger term1 = n.multiply(n).subtract(n).add(BigInteger.valueOf(12));
        BigInteger term2 = n.add(BigInteger.ONE);
        BigInteger term3 = n.add(BigInteger.valueOf(2));
        return term1.multiply(term2).multiply(term3).divide(BigInteger.valueOf(24));
    }

    private static BigInteger F3(BigInteger n) {
        BigInteger term1 = n.multiply(n).subtract(n).add(BigInteger.valueOf(20));
        BigInteger term2 = n.add(BigInteger.ONE);
        BigInteger term3 = n.add(BigInteger.valueOf(2));
        BigInteger term4 = n.add(BigInteger.valueOf(3));
        return term1.multiply(term2).multiply(term3).multiply(term4).divide(BigInteger.valueOf(120));
    }

    private static BigInteger solveExact(long n) {
        BigInteger N = BigInteger.valueOf(n);
        BigInteger ans = BigInteger.ZERO;

        ans = ans.add(F2(N)).add(F2(N.add(BigInteger.ONE))).add(F2(N.subtract(BigInteger.valueOf(2))))
                .add(F1(N)).add(F1(N.add(BigInteger.ONE)));

        long x = 2;
        while (x < n) {
            long c = n / x;
            long next = Math.min(n / c + 1, n);

            ans = ans.add(F3(BigInteger.valueOf(next - 1 + c))).subtract(F3(BigInteger.valueOf(x - 1 + c)));
            if (x > c) {
                ans = ans.subtract(F3(BigInteger.valueOf(next - c - 2))).add(F3(BigInteger.valueOf(x - c - 2)));
            } else {
                ans = ans.add(F3(BigInteger.valueOf(c - x))).subtract(F3(BigInteger.valueOf(c - next)));
            }

            x = next;
        }

        return ans;
    }

    public static String solve() {
        long n = 1000000000000L;
        BigInteger ans = solveExact(n).remainder(BigInteger.valueOf(100000000));
        if (ans.compareTo(BigInteger.ZERO) < 0) {
            ans = ans.add(BigInteger.valueOf(100000000));
        }
        return ans.toString();
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}