Now that we are on the other side of Problem 100, it starts to become exciting. Problem 101 is the first problem we should tackle in Project Euler. A problem which is about polynomials and therefore have links to the currently ongoing problems in UVa Online jugde. The problem reads

If we are presented with the firstkterms of a sequence it is impossible to say with certainty the value of the next term, as there are infinitely many polynomial functions that can model the sequence.

As an example, let us consider the sequence of cube numbers. This is defined by the generating function,

u_{n}=n^{3}: 1, 8, 27, 64, 125, 216, …

Suppose we were only given the first two terms of this sequence. Working on the principle that “simple is best” we should assume a linear relationship and predict the next term to be 15 (common difference 7). Even if we were presented with the first three terms, by the same principle of simplicity, a quadratic relationship should be assumed.

We shall define OP(k,n) to be then^{th}term of the optimum polynomial generating function for the firstkterms of a sequence. It should be clear that OP(k,n) will accurately generate the terms of the sequence forn≤k, and potentially thefirst incorrect term(FIT) will be OP(k,k+1); in which case we shall call it abad OP(BOP).

As a basis, if we were only given the first term of sequence, it would be most sensible to assume constancy; that is, forn≥ 2, OP(1,n) =u_{1}.

Hence we obtain the following OPs for the cubic sequence:

OP(1,n) = 11, 1, 1, 1, …OP(2,n) = 7n-61, 8, 15, …OP(3,n) = 6n^{2}-11n+61, 8, 27, 58, …OP(4,n) =n^{3}1, 8, 27, 64, 125, …

Clearly no BOPs exist fork≥ 4.

By considering the sum of FITs generated by the BOPs (indicated in red above), we obtain 1 + 15 + 58 = 74.

Consider the following tenth degree polynomial generating function:

u_{n}= 1 –n+n^{2}–n^{3}+n^{4}–n^{5}+n^{6}–n^{7}+n^{8}–n^{9}+n^{10}

Find the sum of FITs for the BOPs.

This problem is more or less a standard problem with in control engineering which is my field as well as many many other fields when you want to approximate a set of data with some function.

So if we have the set of n points to which we want to fit a polynomial of degree k. Meaning that we want to find the the set of k+1 coefficients such that we obtain the polynomial . This can be expressed through the following set of linear equations

If we have more datapoints (n > k) then we have an overdetermined system and we don’t have a solution that fits the polynomial exactly. However, we can obtain the best fit through some least squares method.

However, if we try to fit a polynomial of k-1 = n degree the matrix is square, and if all the x values are distinct the matrix will be invertible and thus the set of linear equations will have a solution. So to solve this problem we could use pen and pencil to solve the 10 sets of linear equations starting with

giving us the equation polynomial for 1 point. Followed by the 2 point set of linear equations

and so on. Something which is a rather easy task once you know linear algebra. However, I am certain that I couldn’t do it in 1 minute. So either we should look into linear algebra packages which could most certainly do the trick. Solving a set of equations like this has a Big-O complexity of Or we could find another method get the answer.

## Lagrange interpolating polynomial

Even though there are several other methods which can solve the problem I chose Lagrange interpolation. Due to two things, first of all my trusty old algorithm book Introduction to algorithms described it pretty well and because the information about it on wikipedia and Wolfram Mathworld was pretty good.

The method is called Lagrange interpolating polynomial or just Lagrange polynomial. However, it was first discovered by Waring in 1779, then rediscovered by Euler in 1783 and finally published by Lagrange in 1785. So it should probably have a name like Waring-Euler-Lagrange interpolating polynomial. But I guess fame doesn’t always go to the discoverer.

Given a polynomial P(x) of degree k, this can be represented by a point-value representation using n= k+1 points . Given this representation of the polynomial we can evaluate it for any point with the formula

The formula is not nearly as bad as it looks. It can all be done in a few for loops. According to Introduction to algorithms it runs in time which is a lot better than the first approach. However, for such small problems as we have here, I am pretty certain that it wont mean anything.

## Implementation in C#

For our specific problem we know that , and so on, so that makes the formula a little easier to deal with. Furthermore we know that we are only dealing with integers, so with a few tweaks here and there we can actually get the algorithm to run only on integers.

For a specific evaluation in the point we can evaluate the polynomial and thus get the FIT with the following C# implementation.

long result = 0; for (long i = 1; i <= n; i++) { long temp1 = 1; long temp2 = 1; for (long j = 1; j <= n; j++) { if (i == j) { continue; } else { temp1 *= n+1 - j; temp2 *= i - j; } } result += temp1 * y[i-1]/temp2; }

The outer for-loop corresponds to the sum, the inner for-loop corresponds to the two products in the formula. As you can see I have dealt with the products on the fraction separately and then done the division in the end. This allowed me to avoid floating point operations.

All we need now is to generate the points which Bjarki made a nice little polynomial class for in order to solve the UVa online judge problems. We also need to loop through all the all the points 1 to 10. However this is more or less trivial and you can see it in the source code.

## Wrapping up

I have shown you one method for finding all the FITs of the problem. It runs in

The sum of FITs is 37076114526 Solution took 0,2098 ms

which is perfectly fine for me. I know there are several other methods for evaluating the polynomial as well such as finite difference which we used in Problem 28.

You can find the source code for this problem right here. And last but not least, since I know there are several ways of solving this problem, how did you do it?

I solved it with a system of linear equations, the same way as you proposed.

My Mathematica code:

Yes, in a language that supports that kind of operations I would definitively go for that kind of solution. However, C# is not well suited for it unless you start implementing your own matrix operations

When you’ve got 10 buddies who have blogs, then you will need to visit those blogs to view any updated content.

You probably have a feed reader, it could actually

examine those 10 blogs each hour (or everytime you need), let

you understand when any have been up to date, and offer you the brand new content material.

I first did it the linear systems way, but after reading your Lagrange-based method (which is much more straightforward), I came up with some optimizations.

Let [latex]$P_{j,n}(x)$[/latex] be the jth term of the n-degree Lagrange polynomial. From Mathworld, we have

[latex]$P_{j,n}(x) = y_j \prod_{k\neq j}^{n} \frac{x-x_k}{x_j-x_k}$[/latex]

We know that [latex]$x_k = k$[/latex]. We also know we only need to evaluate the polynomial at [latex]$x_{n+1}$[/latex]. Expanding this out, we can derive the recurrence relation:

[latex]$P_{j,n}(n+1) = -P_{j,n-1}(n)\frac{n}{n+1-j}, j<n$[/latex]

and

[latex]$P_{n,n}(n+1) = ny_n$[/latex]

With this, we only need two loops. The first goes from n=1..N, for the degree of the Lagrange polynomial. The second goes from j=1..(n-1) to update and sum the terms of the polynomial evaluated at x=n+1. Given a function poly(a,x) which evaluates the generating polynomial, we have

vector<long long> p(N,0); // terms of Lagrange

long long e = 0; // sum of FITs

long long u = poly(a, 1); // u_1

for (int n=1; n<N; ++n) {

long long u_next = poly(a,n+1); // next u

// update evaluated polynomial and sum terms

p[n-1] = n*u;

long long s = p[n-1];

for (int j=1; j<n; ++j) {

p[j-1] = -p[j-1]*n/(n+1-j);

s += p[j-1];

}

cout << "u_" << (n+1) << "=" << u_next << ", Lagrange=" << s << std::endl;

e += s; // add incorrect term

u = u_next;

}

cout << e << std::endl;

I first did it the linear systems way, but after reading your Lagrange-based method (which is much more straightforward), I came up with some optimizations.

Let [latex]P_{j,n}(x)[/latex] be the jth term of the n-degree Lagrange polynomial. From Mathworld, we have

[latex]P_{j,n}(x) = y_j \prod_{k\neq j}^{n} \frac{x-x_k}{x_j-x_k}[/latex]

We know that [latex]x_k = k[/latex]. We also know we only need to evaluate the polynomial at [latex]x_{n+1}[/latex]. Expanding this out, we can derive the recurrence relation:

[latex]P_{j,n}(n+1) = -P_{j,n-1}(n)\frac{n}{n+1-j}, j<n[/latex]

and

[latex]P_{n,n}(n+1) = ny_n[/latex]

With this, we only need two loops. The first goes from n=1..N, for the degree of the Lagrange polynomial. The second goes from j=1..(n-1) to update and sum the terms of the polynomial evaluated at x=n+1. Given a function poly(a,x) which evaluates the generating polynomial, we have

vector<long long> p(N,0); // terms of Lagrange

long long e = 0; // sum of FITs

long long u = poly(a, 1); // u_1

for (int n=1; n<N; ++n) {

long long u_next = poly(a,n+1); // next u

// update evaluated polynomial and sum terms

p[n-1] = n*u;

long long s = p[n-1];

for (int j=1; j<n; ++j) {

p[j-1] = -p[j-1]*n/(n+1-j);

s += p[j-1];

}

cout << "u_" << (n+1) << "=" << u_next << ", Lagrange=" << s << std::endl;

e += s; // add incorrect term

u = u_next;

}

cout << e << std::endl;