Problem 101: Optimum Polynomial

Detailed mathematical approach

Problem Summary

Infinitely many polynomials fit any finite list of values, so the problem follows the principle that simple is best. For the first \(k\) terms \(u_1,\dots,u_k\) of a sequence, \(\operatorname{OP}(k,n)\) is the \(n\)-th term of the optimum polynomial for these terms, the polynomial of the lowest degree that produces them; for a single term it is the constant \(\operatorname{OP}(1,n)=u_1\). By construction \(\operatorname{OP}(k,n)=u_n\) for \(n\le k\), so the first term that can be wrong is \(\operatorname{OP}(k,k+1)\). If it differs from \(u_{k+1}\), the optimum polynomial is a bad OP (BOP), and \(\operatorname{OP}(k,k+1)\) is its first incorrect term (FIT). For the cubes \(u_n=n^3\) the optimum polynomials are \(1\), \(7n-6\), \(6n^2-11n+6\), and \(n^3\) for every \(k\ge4\); only the first three are bad, with the FITs \(1\), \(15\) and \(58\), which add up to \(74\). ^[1] We must find the sum of the FITs of all BOPs for the sequence generated by the polynomial of degree ten

$$\begin{aligned}u_n={}&1-n+n^2-n^3+n^4-n^5\\&+n^6-n^7+n^8-n^9+n^{10}.\end{aligned}$$

The whole solution in one sentence: write the first eleven terms in a row, below them the differences of neighbouring numbers, below those their differences, and so on; the optimum polynomial for the first \(k\) terms predicts the next term minus the first number of the \(k\)-th row of differences, so it is bad exactly when that number is not \(0\), which for this sequence happens for every \(k\) from \(1\) to \(10\) and for no larger \(k\), and the requested sum is the sum of these ten predictions.

Mathematical Approach

We first see why a few terms cannot determine a sequence and why the optimum polynomial is unique, then compute FITs in three ways and prove that they agree; the last way needs nothing but subtractions.

1. Many rules fit the same terms

Take the first two cubes, \(1\) and \(8\). The line \(7n-6\) passes through both and predicts \(15\) as the third term. But for every number \(c\) the polynomial

$$7n-6+c\,(n-1)(n-2)$$

passes through both as well, because the added part vanishes at \(n=1\) and at \(n=2\). At \(n=3\) it predicts \(15+2c\): the choice \(c=6\) gives the true cube \(27\), the choice \(c=-7\) gives \(1\), and every other value is possible too. So a choice has to be made, and the problem takes the polynomial of the lowest degree. For two terms this is the line, since no constant takes both values \(1\) and \(8\), while every \(c\ne0\) raises the degree to \(2\). Incidentally, \(c=6\) gives \(6n^2-11n+6\), the optimum polynomial for the first three cubes.

2. Exactly one polynomial of degree at most \(k-1\)

For any values \(u_1,\dots,u_k\) there is exactly one polynomial of degree at most \(k-1\) that takes the value \(u_n\) at \(n\) for \(n=1,\dots,k\). ^[2] If \(P\) and \(Q\) were two such polynomials, their difference \(P-Q\) would have degree at most \(k-1\) and the \(k\) zeros \(1,2,\dots,k\). But a polynomial that is not identically zero has at most as many zeros as its degree, because every zero \(a\) splits off a factor \(n-a\); so \(P=Q\). Existence follows from Lagrange's formula in Section 4. Since a polynomial of degree at most \(k-1\) fits, the lowest possible degree is at most \(k-1\), and the optimum polynomial is this unique one. Its degree may be lower, as for the cubes with \(k\ge5\).

3. First approach: solve for the coefficients

The most direct way to find \(\operatorname{OP}(3,n)=a+bn+cn^2\) for the cubes is to write down the three conditions for \(n=1,2,3\):

$$\begin{aligned}a+b+c&=1,\\a+2b+4c&=8,\\a+3b+9c&=27.\end{aligned}$$

Subtracting each equation from the next removes \(a\) and leaves \(b+3c=7\) and \(b+5c=19\); one more subtraction gives \(2c=12\). So \(c=6\), \(b=-11\) and \(a=6\), and the prediction \(\operatorname{OP}(3,4)=96-44+6=58\) differs from \(4^3=64\): the third optimum polynomial is bad, and its FIT is \(58\). For \(k\) terms the same method solves \(k\) linear equations, the one for \(n\) with the coefficients \(1,n,n^2,\dots,n^{k-1}\) (a Vandermonde system), by Gaussian elimination, which carries out such subtractions systematically and costs about \(k^3\) operations with fractions for every \(k\). ^[2] The right-hand sides \(7\) and \(19\) produced by the subtractions are the differences of \(1\), \(8\) and \(27\), and \(12\) is the difference of these differences; Section 5 works with such differences directly.

4. Second approach: Lagrange's formula at the next place

Lagrange's formula adds up one piece for each known term, a piece that vanishes at the other places. ^[3] For the first two cubes it reads

$$\begin{aligned}\operatorname{OP}(2,x)&=1\cdot\frac{x-2}{1-2}+8\cdot\frac{x-1}{2-1}\\&=-(x-2)+8(x-1)=7x-6.\end{aligned}$$

The first fraction is \(1\) at \(x=1\) and \(0\) at \(x=2\), the second one the other way round, so the sum takes the values \(1\) and \(8\) there; at \(x=3\) the fractions are \(-1\) and \(2\), which gives \(-1+2\cdot8=15\). For \(k\) terms each piece is a product, and \(\prod\) multiplies its factors as \(\sum\) adds its terms:

$$\operatorname{OP}(k,x)=\sum_{i=1}^{k}u_i\prod_{j\ne i}\frac{x-j}{i-j},$$

where \(j\) runs over \(1,\dots,k\) without \(i\). The product that belongs to \(u_i\) is \(1\) at \(x=i\) and \(0\) at the other places \(1,\dots,k\), where one of its factors vanishes, so the sum takes the value \(u_i\) at \(x=i\); and it has degree at most \(k-1\), because each product has \(k-1\) factors. This proves existence. We only need the value at \(x=k+1\). There the numerators \(k+1-j\) run through \(k,k-1,\dots,1\) except \(k+1-i\), so they multiply to \(k!/(k+1-i)\). The denominators \(i-j\) are \(i-1,\dots,1\) for \(j\lt i\) and \(-1,\dots,-(k-i)\) for \(j\gt i\), so they multiply to \((-1)^{k-i}(i-1)!\,(k-i)!\). In terms of the binomial coefficient \(\binom{k}{m}=k!/(m!\,(k-m)!)\), the number of ways to choose \(m\) of \(k\) things, ^[4] the quotient is \((-1)^{k-i}\binom{k}{i-1}\), because \((k+1-i)\cdot(k-i)!=(k+1-i)!\); so

$$\operatorname{OP}(k,k+1)=\sum_{i=1}^{k}(-1)^{k-i}\binom{k}{i-1}\,u_i.$$

The weights are row \(k\) of Pascal's triangle without its last entry, with alternating signs: \(-1,2\) for \(k=2\), \(1,-3,3\) for \(k=3\) and \(-1,4,-6,4\) for \(k=4\). For the cubes, \(\operatorname{OP}(3,4)\)\({}=1-3\cdot8+3\cdot27\)\({}=58\) is the last FIT of the statement, while \(\operatorname{OP}(4,5)=-1+4\cdot8\)\({}-6\cdot27+4\cdot64=125\) is the true fifth cube, so the fourth optimum polynomial is not bad.

5. Third approach: differences lower the degree

Write \(\Delta u_n=u_{n+1}-u_n\) for the difference of two neighbouring terms, and \(\Delta^2u_n=\Delta(\Delta u_n)\), \(\Delta^3u_n\) and so on for differences of differences. ^[5] For the cubes, each row below holds the differences of the row above it:

$$\begin{matrix}1&&8&&27&&64&&125\\&7&&19&&37&&61&\\&&12&&18&&24&&\\&&&6&&6&&&\\&&&&0&&&&\end{matrix}$$

The third differences are constant and the fourth are zero, and the same happens for every polynomial: after as many steps as its degree the differences are constant, and one step later they vanish. If \(p\) is a polynomial of degree \(r\ge1\) with leading coefficient \(a\), the binomial theorem ^[4] gives \((n+1)^r=n^r+r\,n^{r-1}+\dots\), where the dots stand for lower powers of \(n\). So \(\Delta(a\,n^r)=a\,r\,n^{r-1}+\dots\), and the lower terms of \(p\) only produce lower powers as well. Each difference therefore lowers the degree by exactly one and multiplies the leading coefficient by the old degree, while a constant has the difference \(0\). After \(r\) steps only the constant \(r!\,a\) is left, and the next differences vanish. For the cubes \(r=3\) and \(a=1\), and indeed \(3!=6\).

6. The \(k\)-th difference as a signed sum

Unfolding the definition gives \(\Delta^2u_1=u_3-2u_2+u_1\) and \(\Delta^3u_1=u_4-3u_3+3u_2-u_1\). In general,

$$\Delta^ku_1=\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}\,u_{1+i}.$$

The proof is by induction on \(k\), for all sequences at once; for \(k=1\) it is the definition \(\Delta u_1=u_2-u_1\). For the step we write \(\Delta^{k+1}u_1\)\({}=\Delta^ku_2-\Delta^ku_1\) and use the formula for both terms, for the first one applied to the sequence \(u_2,u_3,\dots\). With the convention \(\binom{k}{r}=0\) for \(r\lt0\) and for \(r\gt k\), the term \(u_{1+i}\), for \(0\le i\le k+1\), then receives the coefficients \(\binom{k}{i-1}\) and \(\binom{k}{i}\), both with the sign \((-1)^{k+1-i}\), and Pascal's rule \(\binom{k}{i-1}+\binom{k}{i}\)\({}=\binom{k+1}{i}\), ^[6] which with this convention also holds for \(i=0\) and \(i=k+1\), combines them into the coefficient that the formula for \(k+1\) requires. For the cubes, \(\Delta^3u_1=64-3\cdot27+3\cdot8-1\)\({}=6\), the first third difference in the table above.

7. The key identity: the prediction is \(u_{k+1}-\Delta^ku_1\)

Let \(P\) be the optimum polynomial for the first \(k\) terms. It has degree at most \(k-1\), so by Section 5 its \(k\)-th differences are zero, and by Section 6

$$\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}\,P(1+i)=0.$$

The same sum with \(u_{1+i}\) in place of \(P(1+i)\) is \(\Delta^ku_1\). Since \(P(1+i)=u_{1+i}\) for \(i\lt k\), the two sums differ only in their last terms, \(P(k+1)\) and \(u_{k+1}\), which both have the coefficient \(1\). Subtracting the two sums, and writing \(\operatorname{OP}(k,k+1)\) for \(P(k+1)\), gives

$$\operatorname{OP}(k,k+1)=u_{k+1}-\Delta^ku_1.$$

So the optimum polynomial for \(k\) terms is bad exactly when \(\Delta^ku_1\ne0\), and its FIT is then \(u_{k+1}-\Delta^ku_1\); solving the zero sum for \(P(k+1)\) instead gives back the weights of Section 4. For the cubes the first numbers of the rows are \(1,7,12,6,0\), so the FITs are \(8-7=1\), \(27-12=15\) and \(64-6=58\), and \(\Delta^4u_1=0\) shows that the fourth optimum polynomial is not bad. For a sequence given by a polynomial of degree \(d\) with leading coefficient \(a\), \(\Delta^du_1=d!\,a\ne0\) and all later differences vanish: the optimum polynomial for \(d\) terms is always bad, and those for more than \(d\) terms, which are the generating polynomial itself, never are.

8. Newton's form: one new term for every new value

The binomial polynomials \(\binom{x}{j}=x(x-1)\cdots(x-j+1)/j!\), with \(\binom{x}{0}=1\), have degree \(j\). Pascal's rule, which holds for every number in place of \(x\), ^[6] gives \(\Delta\binom{x}{j}=\binom{x}{j-1}\) for \(j\ge1\), while \(\Delta\binom{x}{0}=0\), and \(\binom{0}{j}\) is \(1\) for \(j=0\) and \(0\) for \(j\ge1\). Every polynomial \(p\) of degree at most \(d\) is a combination \(\sum_jb_j\binom{x}{j}\) with \(j\le d\): subtracting \(d!\,a\binom{x}{d}\), where \(a\) is the coefficient of \(x^d\) in \(p\), leaves degree at most \(d-1\), and we repeat. Taking \(j\) differences turns \(\binom{x}{m}\) into \(\binom{x}{m-j}\) for \(m\ge j\) and into \(0\) for \(m\lt j\), so setting \(x=0\) afterwards leaves exactly \(b_j\). Hence \(b_j=\Delta^jp(0)\), which is Newton's forward difference formula: ^[5]

$$p(x)=\sum_{j=0}^{d}\binom{x}{j}\,\Delta^jp(0).$$

Now let \(P\) be the optimum polynomial for the first \(k\) terms of any sequence. The formula with \(d=k-1\) for \(p(x)=P(1+x)\), whose differences at \(0\) are the numbers \(\Delta^jP(1)\), gives \(P(n)\) at \(x=n-1\) as a combination of these numbers for \(j\lt k\). By Section 6, \(\Delta^jP(1)\) uses only \(P(1),\dots,P(j+1)\), which for \(j\lt k\) are the terms \(u_1,\dots,u_{j+1}\); so \(\Delta^jP(1)=\Delta^ju_1\), and for every sequence

$$\operatorname{OP}(k,n)=\sum_{j=0}^{k-1}\binom{n-1}{j}\,\Delta^ju_1.$$

For the cubes, \(\operatorname{OP}(3,n)=1+7(n-1)\)\({}+12\binom{n-1}{2}\), which multiplies out to \(6n^2-11n+6\). Going from \(k\) to \(k+1\) terms adds the single term \(\binom{n-1}{k}\Delta^ku_1\), which equals \(\Delta^ku_1\) at \(n=k+1\): the identity of Section 7 once more. At \(n=k+1\) the formula becomes

$$\operatorname{OP}(k,k+1)=\sum_{j=0}^{k-1}\binom{k}{j}\,\Delta^ju_1,$$

for example \(1+3\cdot7+3\cdot12=58\) for the cubes and \(k=3\). The C++ program evaluates this form.

9. All FITs at once

When the optimum polynomials for \(1,2,\dots,d\) terms are all bad, the FITs of Section 7 add up to a closed form, because \(\Delta^0u_1=u_1\) can be added to both sums:

$$\begin{aligned}&\sum_{k=1}^{d}\bigl(u_{k+1}-\Delta^ku_1\bigr)\\&\qquad=\sum_{n=1}^{d+1}u_n-\sum_{k=0}^{d}\Delta^ku_1.\end{aligned}$$

In words: the sum of the FITs is the sum of the first \(d+1\) terms minus the sum of the first numbers of all rows of differences. For the cubes, \((1+8+27+64)\)\({}-(1+7+12+6)\)\({}=100-26=74\).

10. A complete example of degree four

The same family with degree \(4\) can be done entirely by hand. The polynomial \(u_n=1-n+n^2-n^3+n^4\) gives the terms \(1,11,61,205,521\) for \(n=1,\dots,5\), and the table of differences is

$$\begin{matrix}1&&11&&61&&205&&521\\&10&&50&&144&&316&\\&&40&&94&&172&&\\&&&54&&78&&&\\&&&&24&&&&\end{matrix}$$

The last number is \(4!=24\), as Section 5 predicts. The first numbers of the rows are \(1,10,40,54,24\); since \(\Delta^ku_1\ne0\) for \(k=1,\dots,4\), all four optimum polynomials are bad, with the FITs \(11-10=1\), \(61-40=21\), \(205-54=151\) and \(521-24=497\). Their sum is \(670\), and the closed form of Section 9 agrees: \(799-129=670\). The same computation for the degrees \(2\) and \(3\) gives the sums \(6\) and \(-55\).

11. The polynomial of the problem

Each \(u_n\) is a geometric sum with the ratio \(-n\); its eleven summands add up to ^[7]

$$\begin{aligned}u_n&=\sum_{r=0}^{10}(-n)^r=\frac{1-(-n)^{11}}{1+n}\\&=\frac{n^{11}+1}{n+1}.\end{aligned}$$

This makes the first terms easy to check by hand: \(u_1=2/2=1\), \(u_2=2049/3=683\), \(u_3=177148/4=44287\) and \(u_4=4194305/5=838861\). The polynomial has degree \(10\) and leading coefficient \(1\), so the tenth differences are all \(10!=3628800\) and the eleventh are \(0\). Only \(k=1,\dots,10\) can contribute, and the eleven terms \(u_1,\dots,u_{11}\) contain all the information needed; the largest, \(u_{11}=(11^{11}+1)/12\), has eleven digits. The first numbers of the rows of differences begin with \(1,682,42922,708048\), so the first three FITs are \(\operatorname{OP}(1,2)=683-682=1\), \(\operatorname{OP}(2,3)=44287-42922\)\({}=1365\) and \(\operatorname{OP}(3,4)=838861-708048\)\({}=130813\).

12. Why all ten predictions fail

It remains to show that \(\Delta^ku_1\ne0\) for every \(k\) from \(1\) to \(10\), and there is a short proof. Group the terms of \(u_n\) in pairs,

$$u_n=1+(n^2-n)+(n^4-n^3)+\dots+(n^{10}-n^9),$$

and put \(n=1+x\). Each pair \(n^{2q}-n^{2q-1}=n^{2q-1}(n-1)\) becomes \(x\,(1+x)^{2q-1}\), so

$$u_{1+x}=1+x\sum_{q=1}^{5}(1+x)^{2q-1}.$$

This polynomial in \(x\) has the coefficients \(c_0,\dots,c_{10}\) equal to \(1,5,25,70,130,166,\)\(148,91,37,9,1\), all positive. Next, a function from a set of \(j\) elements into a set of \(x\) elements is a list of \(j\) values, each one of \(x\) choices, so there are \(x^j\) of them. The values that occur form its image, and a function is onto a set if every element of the set occurs. Let \(\sigma(j,k)\) be the number of functions from a set of \(j\) elements onto a set of \(k\) elements, \(k!\) times a Stirling number of the second kind; ^[8] so \(\sigma(j,k)\ge0\) and \(\sigma(k,k)=k!\). Sorting the \(x^j\) functions by their image, one of the \(\binom{x}{k}\) subsets of \(k\) elements, gives

$$x^j=\sum_{k=0}^{j}\sigma(j,k)\binom{x}{k}.$$

For \(j=3\) and \(x=2\), of the \(8\) lists of three values from \(\{1,2\}\), \(2\) are constant, \(\binom{2}{1}=2\) choices of the value times \(\sigma(3,1)=1\), and \(6\) use both values, \(\binom{2}{2}=1\) times \(\sigma(3,2)=2^3-2=6\); so \(8=2\cdot1+1\cdot6\). With \(\sigma(3,3)=3!=6\) the case \(j=3\) reads \(x^3=\binom{x}{1}+6\binom{x}{2}\)\({}+6\binom{x}{3}\). Both sides of the general formula are polynomials in \(x\) that agree at every whole number \(x\ge0\), so by Section 2 they are equal. Hence \(u_{1+x}=\sum_kB_k\binom{x}{k}\) with \(B_k=\sum_jc_j\,\sigma(j,k)\), and Section 8 shows that \(B_k=\Delta^ku_1\). No term of \(B_k\) is negative, and the term with \(j=k\) is \(c_k\,k!\), so

$$\Delta^ku_1\ge c_k\,k!\gt0\quad(1\le k\le10).$$

In the example of degree four, \(u_{1+x}=1+2x+4x^2+3x^3+x^4\), and \(\sigma(2,2)=2\), \(\sigma(3,2)=6\), \(\sigma(4,2)=14\) give \(\Delta^2u_1=4\cdot2+3\cdot6+1\cdot14\)\({}=40\), as in Section 10. So all ten optimum polynomials of the problem are bad, and by Section 9 the requested sum is the sum of \(u_1,\dots,u_{11}\) minus the sum of \(\Delta^0u_1,\dots,\Delta^{10}u_1\).

13. A pitfall: a correct next term

A correct next term does not make an optimum polynomial correct forever. For \(u_n=n^2-3n+3\), with the terms \(1,1,3,7\), the constant \(\operatorname{OP}(1,n)=1\) predicts \(u_2=1\) correctly, as \(\Delta u_1=0\) says, and fails first at \(n=3\). The statement calls an optimum polynomial bad when its first incorrect term is the next one, \(\operatorname{OP}(k,k+1)\), so this one is not a BOP, although it is wrong later. Only the optimum polynomial for two terms, again the constant \(1\), is bad, with the FIT \(1\). Adding \(\operatorname{OP}(k,k+1)\) for \(k=1,2\) without comparing gives \(2\), and so does the closed form of Section 9, which assumes that every \(k\) is bad. All three programs compare before they add, and only at \(k+1\). For the polynomial of the problem the case cannot occur: by Section 12 each optimum polynomial for \(k\le10\) terms is already wrong at \(k+1\), and those for more terms are the polynomial itself, so a reading that also counted later failures would find the same ten BOPs and FITs.

How the Code Works

The C++ program ^[9] follows Sections 5 and 8 with exact integers of type __int128_t, abbreviated as i128; this is another name for the 128-bit type __int128 that GCC offers as an extension. ^[10] main calls parse_arguments, which accepts --skip-checkpoints and --degree=D for \(D\ge1\) written in decimal digits, by default \(10\), and rejects every other argument. It does not check whether \(D\) fits into the int that collects its digits, and the overflow of a signed integer is undefined behaviour in C++, ^[11] so a far too large \(D\) is not cleanly rejected. Unless the checkpoints are skipped, run_checkpoints applies the method to the cubes, given by the coefficients \(0,0,0,1\), and stops the program with an error if the result is not \(74\). Then sum_of_fits_alternating builds the coefficients \(1,-1,1,\dots\) of \(\sum_r(-n)^r\) and passes them to sum_of_fits_from_coeffs.

That function computes \(u_1,\dots,u_{d+1}\) with term_from_coeffs, which adds each coefficient times the current power of \(n\) and then multiplies the power by \(n\). It builds the table of differences row by row and keeps the first number of each row in first_diff_entries, which then holds \(\Delta^0u_1,\dots,\Delta^du_1\). For each \(k\) from \(1\) to \(d\) it evaluates Newton's form at \(n=k+1\) as the sum of binom_i128(n - 1, j) times \(\Delta^ju_1\) over \(j\lt k\), and adds the prediction only if it differs from \(u_{k+1}\), exactly as the definition of a BOP requires. binom_i128 uses the multiplicative rule of the binomial coefficients: ^[4] step \(i\) multiplies \(\binom{n}{i-1}\) by \(n-i+1\), which gives \(i\binom{n}{i}\), and divides by \(i\) without remainder. Finally i128_to_string writes the result in decimal digits, because the standard output streams cannot print a 128-bit integer.

The Python and Java programs use Section 4 instead. Python's solve computes the eleven terms with its inner function u, using unbounded integers and computing every power n**k on its own, and keeps them in the list values; Java's main fills a long[] values with the method u, which multiplies a running power by \(n\) and finally forms the unused power \(n^{11}\), and 64-bit long values are plenty for the largest such number, \(11^{11}\approx2.9\cdot10^{11}\). For each \(k\) from \(1\) to \(10\), two nested loops over \(i\) and \(j\ne i\) evaluate Lagrange's formula at \(k+1\) factor by factor, without the simplification to binomial coefficients; Python reads the places from the list xs and the values from ys, Java writes the places as i + 1.0 and j + 1.0. The divisions produce floating-point numbers, through Python's operator / and Java's double variables, so each prediction is rounded with round or Math.round and added to total only if it differs from values[k], which is \(u_{k+1}\). Both programs have the degree \(10\) built in and contain no self-check.

Every exact prediction is an integer, by Section 4 a sum of integers \(u_i\) times integer weights, so rounding recovers it as long as the floating-point error stays below \(1/2\). A double has a significand of \(53\) bits and represents every integer of absolute value up to \(2^{53}\approx9.0\cdot10^{15}\) exactly, ^[12] but the divisions by \(i-j\) still produce rounded fractions. Compared with exact fractions, no intermediate value in these loops, not even a product before its division, exceeds about \(3.1\cdot10^{11}\), and no prediction is more than about \(3\cdot10^{-5}\) away from the exact one, far below \(1/2\). The margin shrinks fast: on the same family of degree \(12\) the error before rounding grows to about \(0.016\), and for degree \(13\) a rounded prediction is wrong by \(5\). The C++ program needs no rounding, and its limit is the 128-bit range: with --degree=D every value stays inside it up to \(D=25\). At \(D=26\) the last multiplication in term_from_coeffs forms the unused power \(27^{27}\), which exceeds \(2^{127}-1\), and from \(D=27\) on values that the predictions use overflow as well. The program does not check for overflow, and since signed overflow is undefined behaviour, nothing guarantees its output from \(D=26\) on.

Complexity and Verification

Write \(d\) for the degree, here \(10\), and count every arithmetic operation as one step, which is fair here because all numbers involved stay below \(10^{12}\). The \(d+1\) terms cost about \(d^2\) operations in C++ and Java and a little more in Python, which computes every power on its own, and the table of differences costs \(d(d+1)/2\) subtractions. The C++ program spends about \(k^2/2\) multiplications and as many divisions on the prediction for \(k\), since every binomial coefficient is built afresh, about \(d^3/6\) in all, exactly \(220\) multiplications and \(165\) divisions for \(d=10\); the Lagrange loops of Python and Java spend \(k(k-1)\) multiplications and as many divisions, about \(d^3/3\) in all, \(330\) of each for \(d=10\). Solving the systems of Section 3 would cost of the order of \(d^4\). With the identity of Section 7 each prediction is one subtraction, and the whole computation is quadratic in \(d\). The memory is linear in \(d\). For \(d=10\) these are a few hundred operations, and each program finishes within a few hundredths of a second, almost all of it start-up time.

Besides the self-check of the C++ program with the cubes, further checks were run separately for this article. All three programs print the same number, equal to a computation with exact fractions. For \(2000\) random polynomials of degree at most \(9\) with integer coefficients between \(-20\) and \(20\), and every \(k\) up to two more than the degree, five methods give the same prediction: the Vandermonde system and Lagrange's formula with fractions, the binomial weights, Newton's form and \(u_{k+1}-\Delta^ku_1\). For the polynomial of the problem, all ten numbers \(\Delta^ku_1\) are positive and at least \(c_k\,k!\), the eleventh differences vanish, and the closed form of Section 9 equals the sum of the FITs.

Footnotes and References

The references below supply the background facts; the weights of Lagrange's formula at the next place, the identity \(\operatorname{OP}(k,k+1)\)\({}=u_{k+1}-\Delta^ku_1\), the closed form of Section 9 and the proof of Section 12 are derived in this article.

1. Project Euler 101 — Optimum Polynomial. The official statement defines OP(k, n), the first incorrect term (FIT) and the bad OP (BOP), works out the cubes with the FITs 1, 15 and 58 and their sum 74, and asks for the sum of the FITs of the BOPs for u_n = 1 − n + n² − n³ + … + n¹⁰.
2. Wikipedia — Polynomial interpolation. Through n + 1 points with distinct nodes there is a unique polynomial of degree at most n; one proof uses that a non-zero polynomial of degree at most n has at most n zeros. The coefficients satisfy a linear system with a Vandermonde matrix, and Gaussian elimination solves it in O(n³) operations.
3. Wikipedia — Lagrange polynomial. The basis polynomial of a node is the product of (x − x_m)/(x_j − x_m) over the other nodes, and the values times their basis polynomials add up to the unique interpolating polynomial of lowest degree.
4. Wikipedia — Binomial coefficient. C(n, k) = n!/(k!(n − k)!) is the number of ways to choose k of n elements; the page also gives Pascal's triangle, the binomial formula for (x + y)ⁿ and the multiplicative rule C(n, k) = C(n, k − 1)(n − k + 1)/k.
5. Wikipedia — Finite difference. The n-th forward difference is the sum of (−1)^(n−j) C(n, j) f(x + j); for a polynomial of degree n with leading coefficient a it is the constant a·n!, and Newton's series f(x) = Σ C(x − x₀, k) Δ^k f(x₀) holds for every polynomial.
6. Wikipedia — Pascal's rule. C(n − 1, k) + C(n − 1, k − 1) = C(n, k), and the rule holds more generally when n is replaced by any complex number.
7. Wikipedia — Geometric series. For a ratio r ≠ 1 the sum a + ar + … + arⁿ, from the 0th to the nth term, equals a(1 − r^(n+1))/(1 − r).
8. Wikipedia — Stirling numbers of the second kind. There are k! S(n, k) surjections from an n-element set onto a k-element set, and x^n is the sum of S(n, k) times the falling factorials x(x − 1)⋯(x − k + 1); since each of these is k! C(x, k), this gives the expansion of Section 12.
9. C++ — Euler101.cpp. The linked immutable revision contains the functions described above and the options --skip-checkpoints and --degree=D; it is the implementation record for this article.
10. GCC — 128-bit Integers. As an extension, GCC supports __int128 and unsigned __int128 on targets that have an integer mode wide enough to hold 128 bits.
11. C++ working draft — [expr.pre], paragraph 4. The behaviour of evaluating an arithmetic expression is undefined if its mathematical result is not in the range of values that its type can represent (infinities and NaNs of floating-point types aside); an overflowing operation on signed integers is such a case.
12. Wikipedia — Double-precision floating-point format. The significand has 53 bits, and every integer from −2⁵³ to 2⁵³ is represented exactly.