プロジェクトオイラー 問１０１

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１０１「最適多項式」
数列のk個の項を与えられたときに, 次の項を確実に求めることは不可能である. 
その数列に合うような多項式が無限個存在するからである.

例として, 立方数の数列を考えよう. これは生成関数 u_n = n³ で定義され, 1, 8, 27, 64, 125, 216, ...となる.

この数列の最初の2項のみが与えられているとしよう. 
"Simple is best"の法則にのっとり, 線形の関係があると仮定し, 3つ目の項が15であると予想する (差分が7). もし最初の3項のみが与えられていたとしても, 同じ原則により, 二次の関係があると仮定して次の項を予測する.

数列の最初のk項を生成できる最適な多項式のn項を OP(k, n) で表すことにする. 
明らかに, n ≦ k について OP(k, n) は正しい. 
最初の異なる項 (First Incorrect Term, FIT) は OP(k, k+1) であろう. 
これを bad OP (BOP) と呼ぶことにする.

原則より, 最初の項しか与えられていない場合には, 定数項とするのが理に適っているだろう;即ち, n ≧ 2, OP(1, n) = u₁.

従って, 立方数の数列について以下のOPを得る.

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 のときにはBOPは存在しない.
BOPのFIT (上の例では赤で示されている) の和は, 1 + 15 + 58 = 74 である.

以下の10次多項式からなる生成関数を考える:
u_n = 1 - n + n² - n³ + n⁴ - n⁵ + n⁶ - n⁷ + n⁸ - n⁹ + n¹⁰
BOPのFITの総和を求めよ.






注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。






私の回答例は以下の通りです。

def g(Y):
	x, y = len(Y), 0
	Xi = [i for i in xrange(x)]
	for i in xrange(x):
		a, b = Y[i], 1
		for j in xrange(x):
			if i==j: continue
			a *= (x - Xi[j])
			b *= (Xi[i] - Xi[j])
		y += a//b
	return y

def f(s):
	L, FIT, BOP, t, n = [], 0, 0, -1, 0
	while FIT!=t:
		BOP += FIT
		n += 1
		FIT, t = g(L), eval(s)
		L.append(t)
	return BOP

u = "1 - n + n**2 - n**3 + n**4 - n**5 + n**6 - n**7 + n**8 - n**9 + n**10"
BOP = f(u)
print "BOP:", BOP

この問題は解法が思いつかず、ネット検索したところ、「ラグランジュの補間多項式」の問題と判明。
ラグランジュ補間(wiki) を参考にロジックを組みました。

１．関数g(Y)
・ラグランジュの補間多項式を組み込んだ関数です。
　数値リストYから要素数より１つ少ない次元の補間多項式による、次の要素を返します。
・数値xはリストYの要素数で、求める要素の要素番号です。
　リストXiは、リストYの要素番号のリストです。
・ループ変数iは補間式内の項の番号です。
　ループ変数jは補間式の各項内の掛け算要素の番号です。
・変数a, bは、その掛け算要素の分子と分母です。
　項１つ分の分子と分母の計算が終わったら、その商を変数yに累積していきます。

２．関数f(s)
・文字型で数式を受取り、問題文にあるBOPのFITの総和を返します。
　数式の変数は１種類、nだけとします。
・リストLは、受け取り数式sによる実値をためていきます。
・whileループ内は、現在までの実値リストに基づく推定値FITと
　数式による実値tとが不一致の間、処理を実行します。
　また、初回n=1のときは一致することがあるので、実値リストLがカラのときは処理を行うようにします。
・受け取り数式sの変数nに1から順に値を代入していき、関数g()を使って推定値FITを求めます。
　文字列である数式をeval関数で処理する数式として認識させ実行します。



解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
BOP: 37076114526