[玲珑杯#Round5]Grid Point 拉格朗日插值

答案为
$g(n) = \sum_{i=0}^n(2n+1)^k-(2i+1)^k$

令这个k+1次多项式为$f(x)$
令s=k+1
根据拉格朗日插值公式
$f(x) = \sum_{i=0}^s p_i*y_i$
其中
$p_i=\frac {\prod_{i=0}^s (x-xi)} {\prod_{i=0}^s(xj-xi)}$
令$x_i=i$
$f(x) = \sum_{i=0}^s p_i(i)*y_i$
答案即为$f(n)$
$f(n) = \sum_{i=0}^s \frac {\prod_{j=0}^sn-j} {\prod_{j=0}^si-j}y_i$
$f(n) = \sum_{i=0}^s (-1)^{d-i}\frac {n(n-1)...(n-s)} {(s-i)!i!(n-i)} y_i$

今天上课完全没听懂他在说啥
在下面重新整理了一下这道题
手写的书写水平一般拍摄水平垃圾
凑合着看吧

#include <iostream>
#include <cstdio>
#include <cstring>
#define K 2000050
#define p 1000000007LL

using namespace std;
typedef long long LL;
int n,k,s;
LL invjc[K/2],jc[K/2],xk[K],sum_xk[K];
LL y[K],pr[K],cnt;
bool prime[K];
int inv[K],l[K],r[K];

inline void inc(LL &x,LL y){ x = (x + y + p) % p; }

LL qp(LL a,LL b) {
    LL r = 1;
    while (b) {
        if (b&1) r = r * a % p;
        b >>= 1 ,a = a * a % p;
    } 
    return r;
}

void _inv() {
    inv[1] = 1;
    for (LL i=2;i<=s;i++)
        inv[i] = (p - p / i) % p * inv[p % i] % p , inv[i] = (inv[i] + p) % p;

    invjc[0] = invjc[1] = 1;
    for (LL i=2;i<=s;i++)
        invjc[i] = invjc[i-1] * inv[i] % p;
}

void _pre() {
    jc[0] = jc[1] = 1;
    for (LL i=2;i<=s;i++) jc[i] = jc[i-1] * i % p;

    l[0] = n % p;
    for (LL _=1;_<=s;_++) l[_] = 1LL * l[_-1] * (n-_) % p;
    r[s] = (n-s) % p; r[s+1] = 1LL;
    for (LL _=s-1;_>=0;_--) r[_] = 1LL * r[_+1] * (n-_) % p;

    xk[1] = 1;
    for (int i=1;i<=2*s+1;i++) {
        if (!prime[i]) pr[++cnt] = i , xk[i] = qp(i,k);
        for (int j=1;j<=cnt && pr[j]<=i;j++) {
            xk[ i*pr[j] ] = xk[i] * xk[ pr[j] ] % p;
            prime[ i*pr[j] ] = 1;
            if (i % pr[j] == 0) break;
        }
    }

//  LL pres = 0;
//  for(int i=0;i<=k+1;i++) {
//      LL cur = xk[2 * i + 1];
//      y[i] = (1LL * i * cur % p + p - pres) % p;
//      (pres += cur) %= ;
//  }

    sum_xk[1] = 1;
    for (int i=3;i<=2*s+1;i+=2) sum_xk[i] = ( sum_xk[i-2] + xk[i] ) % p;

    for (LL i=1;i<=s;i++) y[i] = (i+1) * xk[2*i+1] % p - sum_xk[2*i+1] , y[i] = (y[i] + p) % p;
}

void init() {
    cnt = 0;
    memset(xk,0,sizeof(xk));
    memset(sum_xk,0,sizeof(sum_xk));
    memset(prime,0,sizeof(prime));
    memset(l,0,sizeof(l));
    memset(r,0,sizeof(r));
    memset(jc,0,sizeof(jc));
    memset(invjc,0,sizeof(invjc));
    memset(y,0,sizeof(y));
}

int main() {
//  freopen("1.in","r",stdin);
    ios::sync_with_stdio(false); 
    int T=0; cin >> T;
    while (T--) {
        init();
        cin >> k >> n; cnt = 0;
        s = k + 1;
        _pre(); 
        _inv();
        LL ans = 0LL;
        ans = 1LL * r[0] * y[0] % p * invjc[s] % p * invjc[0] % p;
        for (int i=1;i<=s;i++){
            LL fh = ( (s-i)%2 == 0 ) ? 1 : -1;
            inc( ans , fh * l[i-1] * r[i+1] % p * y[i] % p * invjc[s-i] % p * invjc[i] % p );
        }
        cout << ans << "\n";
    }
    return 0;
}