数列のk個の項を与えられたときに, 次の項を確実に求めることは不可能である. その数列に合うような多項式が無限個存在するからである.
例として, 立方数の数列を考えよう. これは生成関数 un = n^3 で定義され, 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) = u1.
従って, 立方数の数列について以下のOPを得る.
OP(1, n) = 1 1, 1, 1, 1, ... OP(2, n) = 7n−6 1, 8, 15, ... OP(3, n) = 6n^2−11n+6 1, 8, 27, 58, ... OP(4, n) = n^3 1, 8, 27, 64, 125, ... 明らかに, k ≧ 4 のときにはBOPは存在しない.
BOPのFIT (上の例では赤で示されている) の和は, 1 + 15 + 58 = 74 である.
以下の10次多項式からなる生成関数を考える:
un = 1 − n + n^2 − n^3 + n^4 − n^5 + n^6 − n^7 + n^8 − n^9 + n^10
BOPのFITの総和を求めよ.
65ms．しかしこんな解き方でいいのやら…いいよな．Rだし．
## n個のデータを持つx, y=f(x)にn-1次多項式のモデルを当てはめて係数を返す nth.poly.coefs <- function(x, y){ n <- length(x) model <- "y ~ 1 + x" if(n >= 3){ for(i in 2:(n-1)){ model <- paste(model, "+I(x^", i , ")") } } data <- data.frame(x, y) round(lm(model, data=data)$coefficients) } ## n個の係数をn-1次多項式f(x)に復元したときのf(n+1)の値を求める get.FIT <- function(coefs){ n <- length(coefs) FIT <- coefs[1] if(n > 1){ for(i in 2:n){ FIT <- FIT + coefs[i] * (n + 1)^(i - 1) } } return(FIT) } gen <- function(n) 1 - n + n^2 - n^3 + n^4 - n^5 + n^6 - n^7 + n^8 - n^9 + n^10 n <- 1:11 gn <- gen(n) ## FITの総和 sum.FIT <- 1 for(i in 2:10){ x <- n[1:i] y <- gn[1:i] sum.FIT <- sum.FIT + get.FIT(nth.poly.coefs(x, y)) } sum.FIT
