project euler 101

Problem 101

---

Optimum polynomial

If we are presented with the first k terms 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,
un = n3: 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 the nth term of the optimum polynomial generating function for the first k terms of a sequence. It should be clear that OP(k, n) will accurately generate the terms of the sequence for n ≤ k, and potentially the first incorrect term (FIT) will be OP(k, k+1); in which case we shall call it a bad 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, for n ≥ 2, OP(1, n) = u1.

Hence we obtain the following OPs for the cubic sequence:

OP(1, n) = 1	1, 1, 1, 1, …
OP(2, n) = 7n−6	1, 8, 15, …
OP(3, n) = 6n2−11n+6	1, 8, 27, 58, …
OP(4, n) = n3	1, 8, 27, 64, 125, …

Clearly no BOPs exist for k ≥ 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.

---

最优多项式

如果我们知道了一个数列的前k项，我们仍无法确定地给出下一项的值，因为有无穷个多项式生成函数都有可能是这个数列的模型。

例如，让我们考虑立方数的序列，它可以用如下函数生成，
un = n3: 1, 8, 27, 64, 125, 216, …

如果我们只知道数列的前两项，秉承“简单至上”的原则，我们应当假定这个数列遵循线性关系，并且预测下一项为15（公差为7）。即使我们知道了数列的前三项，根据同样的原则，我们也应当首先假定数列遵循二次函数关系。

给定数列的前k项，定义OP(k, n)是由最优多项式生成函数给出的第n项的值。显然OP(k, n)可以精确地给出n ≤ k的那些项，而可能的第一个不正确项（First Incorrect Term，简记为FIT）将会是OP(k, k+1)；如果事实的确如此，我们称这个多项式为坏最优多项式（Bad OP，简记为BOP）。

在最基本的情况下，如果我们只得到了数列的第一项，我们应当假定数列为常数，也就是说，对于n ≥ 2，OP(1, n) = u1。

由此，我们得到了立方数列的最优多项式如下：

OP(1, n) = 1	1, 1, 1, 1, …
OP(2, n) = 7n−6	1, 8, 15, …
OP(3, n) = 6n2−11n+6	1, 8, 27, 58, …
OP(4, n) = n3	1, 8, 27, 64, 125, …

显然，当k ≥ 4时不存在坏最优多项式。

所有坏最优多项式的第一个不正确项（用红色标示的数）之和为1 + 15 + 58 = 74。

考虑下面这个十阶多项式生成函数：

u n = 1 − n + n 2 − n 3 + n 4 − n 5 + n 6 − n 7 + n 8 − n 9 + n 10

求其所有坏最优多项式的第一个不正确项之和。

package projecteuler;

import java.util.Arrays;

import junit.framework.TestCase;

public class Prj101 extends TestCase {

	
	/**
	 * newton interpolation
	 */
	public void testOptimumPolynomial() {

		long[] fx = new long[12];
		long[] x = new long[12];
		for (int i = 0; i <= 11; i++) {
			fx[i] = getVal(i + 1);
			// fx[i] = (long) Math.pow(i + 1, 3);
			x[i] = i + 1;
		}
		long[] diff = calculate(fx, x);

		
		long sum = 0;
		int j = 1;
		for (int i = 1; i < 11; i++) {
			long[] arr = calN(i, x, diff);
			sum += arr[j ++];
			printArr(arr);
		}
		System.out.println("sum=" + sum);

	}

	long[] calN(int k, long[] x, long[] diff) {

		assert (x.length == diff.length);

		int n = x.length;
		long[] ret = new long[n];

		for (int i = 0; i < n; i++) {
			long mul = 1;
			long sum = diff[0];
			for (int j = 1; j < k; j++) {
				mul *= (x[i] - j);
				sum += mul * diff[j];
			}
			ret[i] = sum;
		}
		return ret;

	}

	long getVal(int n) {
		int q = -1 * n;
		return (long) (1L * (Math.pow(q, 11) - 1) / (q - 1));
	}

	public long[] calculate(long[] fx, long[] x) {

		long[] diff = Arrays.copyOf(fx, fx.length);
		int d = 1;
		for (int i = 0; i < x.length; i++) {

			long[] tp = Arrays.copyOf(diff, diff.length);
			for (int j = i + 1; j < x.length; j++) {
				tp[j] = (diff[j] - diff[j - 1]) / (x[j] - x[j - d]);
			}
			d++;

			diff = Arrays.copyOf(tp, diff.length);
		}

		return diff;
	}

	void printArr(long[] diff) {
		for (int i = 0; i < diff.length; i++) {
			System.out.print(diff[i] + ",");
		}
		System.out.println();
	}
}