模板

update：其实拉格朗日插值法的问题在于求分母的复杂度是 (O(n^2)) 的，要是还要求逆元则再多一个 (logp) 变成 (O(n^2logp))，但是当一个多项式要重复使用的时候，也不必求出他的各个系数，只要预处理出各项分母的逆元之后，(O(nlogp)) 处理分子（求出前缀积和后缀积），然后再插值，渐进复杂度和求出系数值的一样。当分母是等差数列之类的有规律的东西，他的分母可能是可以更快求出的。

---

参考资料：https://attack.blog.luogu.org/solution-p4781

因为高斯消元法是n立方的，有些鬼畜问题需要n平方的拉格朗日插值法。

使用逆元的n平方logn预处理，nlogn单次询问的重心拉格朗日插值法：
（假如使用double的话，理论上预处理和询问都少个log）

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;


namespace Lagrange_Interpolation_Polynomial {
    const int MOD=998244353;
    const int MAXN=2000;
    int n,x[MAXN+5],y[MAXN+5];
    int invd[MAXN+5][MAXN+5];//x_i-x_j的逆元
    int kd[MAXN+5];//k-x_i的值
    int pkd,pinvd[MAXN+5];//k-x的积，invd[i]的积
    //不需要计算k-x_i的逆元，因为复杂度其实一样

    int qpow(int x,int n)
    {
        int ans=1;
        while(n)
        {
            if(n&1)
                ans=((ll)ans*x)%MOD;
            x=((ll)x*x)%MOD;
            n>>=1;
        }
        return ans;
    }

    void init_xi_xj()
    {
        for(int i=1; i<=n; i++)
        {
            for(int j=1; j<=n; j++)
            {
                if(i==j)
                    continue;
                else
                {
                    int d=(x[i]-x[j]);
                    if(d<0)
                        d+=MOD;
                    invd[i][j]=qpow(d,MOD-2);
                }
            }
        }
        for(int i=1;i<=n;i++){
            pinvd[i]=1;
            for(int j=1;j<=n;j++){
                if(i==j)
                    continue;
                pinvd[i]=(ll)pinvd[i]*invd[i][j]%MOD;
            }
        }
    }

    void init_k(int k)
    {
        pkd=1;
        for(int i=1; i<=n; i++)
        {
            kd[i]=(k-x[i]);
            if(kd[i]<0)
                kd[i]+=MOD;
            pkd=(ll)pkd*kd[i]%MOD;
        }
    }

    inline int prod(int i,int k)
    {
        ll ans=(ll)pkd*qpow(kd[i],MOD-2)%MOD;
        ans=ans*(pinvd[i])%MOD;

        return ans;
    }

    int lagrange(int k)
    {
        ll ans=0;
        init_k(k);

        for(int i=1; i<=n; i++)
        {
            ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
        }

        return ans;
    }

}

using namespace Lagrange_Interpolation_Polynomial;



int main() {
    int k;
    scanf("%d%d",&n,&k);
    for(int i=1; i<=n; i++)
    {
        scanf("%d%d",&x[i],&y[i]);
    }
    init_xi_xj();
    printf("%d
",lagrange(k));
}

通过n个不同的点直接构造出多项式。

$f(x)=sumlimits_{i=1}^{n} y_i prodlimits_{i e j}frac{x-x_j}{x_i-x_j} $

正确性：代入 (x_i) 只有 (y_i) 右边是1，其他分子都有0，消掉了。

那么可以n平方构造。

需要小心溢出以及卡常数。假如进行同一个多项式的多次求 (f(k)) 的值 ，(x_i-x_j) 的逆元需要多次使用，可以像下面的类似做法预处理他们的差，当然kd和pkd要重新处理。不需要用map存，这样和直接求qpow一个鬼样。
开个invd数组，保存每两个项之间的差的逆元。

注意到prod里面有一项是 (k-x_j) 的积，这里也可以选择直接记录这个积。毕竟模运算比较费时，卡出来的常数可能差几倍

重心拉格朗日插值法

因为是要求逆元的，所以预处理 (O(n^2logn)) 。
然后每次求新k时，预处理重心拉格朗日插值法的g， (O(nlogn)) 。

多次求同一多项式的值快了不少的，比原版的拉格朗日插值。就是要 (O(n^2)) 的额外空间。

先卡掉重复求差的逆元的，已经可以通过了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int MOD=998244353;

int n,x[2005],y[2005];

int qpow(int x,int n){
    int ans=1;
    while(n){
        if(n&1)
            ans=((ll)ans*x)%MOD;
        x=((ll)x*x)%MOD;
        n>>=1;
    }
    return ans;
}

int invd[2005][2005];
int kd[2005];

void init_prod(int k){
    for(int i=1;i<=n;i++){
        kd[i]=(k-x[i]);
        if(kd[i]<MOD)
            kd[i]+=MOD;
        for(int j=1;j<=n;j++){
            if(i==j)
                continue;
            else{
                int d=x[i]-x[j];
                if(d<MOD)
                    d+=MOD;
                invd[i][j]=qpow(d,MOD-2);;
            }
        }
    }
}

int prod(int i,int k) {
    ll ans=1;
    for(int j=1;j<=n;j++){
        if(j==i)
            continue;
        else{
            ans=ans*(kd[j])%MOD*(invd[i][j])%MOD;
        }
    }

    return ans;
}

int lagrange(int k) {
    ll ans=0;
    init_prod(k);

    for(int i=1;i<=n;i++){
        ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
    }

    return ans;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    int k;
    scanf("%d%d",&n,&k);
    for(int i=1;i<=n;i++){
        scanf("%d%d",&x[i],&y[i]);
    }
    printf("%d
",lagrange(k));
}

再卡掉n倍logn，效果不是很明显，毕竟是n平方的算法：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int MOD=998244353;

int n,x[2005],y[2005];

int qpow(int x,int n){
    int ans=1;
    while(n){
        if(n&1)
            ans=((ll)ans*x)%MOD;
        x=((ll)x*x)%MOD;
        n>>=1;
    }
    return ans;
}

int invd[2005][2005];
int kd[2005],pkd;

void init_prod(int k){
    pkd=1;
    for(int i=1;i<=n;i++){
        kd[i]=(k-x[i]);
        if(kd[i]<MOD)
            kd[i]+=MOD;
        pkd=(ll)pkd*kd[i]%MOD;
        for(int j=1;j<=n;j++){
            if(i==j)
                continue;
            else{
                int d=x[i]-x[j];
                if(d<MOD)
                    d+=MOD;
                invd[i][j]=qpow(d,MOD-2);;
            }
        }
    }
}

int prod(int i,int k) {
    ll ans=(ll)pkd*qpow(kd[i],MOD-2)%MOD;
    for(int j=1;j<=n;j++){
        if(j==i)
            continue;
        else{
            ans=ans*(invd[i][j])%MOD;
        }
    }

    return ans;
}

int lagrange(int k) {
    ll ans=0;
    init_prod(k);

    for(int i=1;i<=n;i++){
        ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
    }

    return ans;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    int k;
    scanf("%d%d",&n,&k);
    for(int i=1;i<=n;i++){
        scanf("%d%d",&x[i],&y[i]);
    }
    printf("%d
",lagrange(k));
}

也可以通过之前用高斯消元法做的https://codeforces.com/contest/1155/problem/E，但是高斯消元法在处理多次询问的时候貌似复杂度好一些？（其实是我写挂了，询问变成n平方log的了，已订正）

实测比高斯消元快一些，估计是高斯消元写得不好。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;


const int MOD=1000003;
int inv[MOD+5];

namespace Lagrange_Interpolation_Polynomial {
    const int MAXN=10;
    int n,x[MAXN+5],y[MAXN+5];
    int invd[MAXN+5][MAXN+5];//x_i-x_j的逆元
    int kd[MAXN+5];//k-x_i的值
    int pkd,pinvd[MAXN+5];//k-x的积，invd[i]的积
    //不需要计算k-x_i的逆元，因为复杂度其实一样

    int qpow(int x,int n)
    {
        int ans=1;
        while(n)
        {
            if(n&1)
                ans=((ll)ans*x)%MOD;
            x=((ll)x*x)%MOD;
            n>>=1;
        }
        return ans;
    }

    void init_xi_xj()
    {
        for(int i=1; i<=n; i++)
        {
            for(int j=1; j<=n; j++)
            {
                if(i==j)
                    continue;
                else
                {
                    int d=(x[i]-x[j]);
                    if(d<0)
                        d+=MOD;
                    invd[i][j]=inv[d];
                }
            }
        }
        for(int i=1;i<=n;i++){
            pinvd[i]=1;
            for(int j=1;j<=n;j++){
                if(i==j)
                    continue;
                pinvd[i]=(ll)pinvd[i]*invd[i][j]%MOD;
            }
        }
    }

    void init_k(int k)
    {
        pkd=1;
        for(int i=1; i<=n; i++)
        {
            kd[i]=(k-x[i]);
            if(kd[i]<0)
                kd[i]+=MOD;
            pkd=(ll)pkd*kd[i]%MOD;
        }
    }

    inline int prod(int i,int k)
    {
        ll ans=(ll)pkd*inv[kd[i]]%MOD;
        ans=ans*(pinvd[i])%MOD;

        return ans;
    }

    int lagrange(int k)
    {
        ll ans=0;
        init_k(k);

        for(int i=1; i<=n; i++)
        {
            ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
        }

        return ans;
    }

}

using namespace Lagrange_Interpolation_Polynomial;


ll jury(ll x) {
    ll a[11]= {10002,0,1};
    ll ret=0;
    ll tx=1;
    for(int i=0; i<=10; i++) {
        ret=(ret+tx*a[i]%MOD)%MOD;
        tx=tx*x%MOD;
    }
    //cout<<"jury: "<<ret<<endl;
    return ret;
}


int query(int x) {
    printf("? %d
",x);
    fflush(stdout);
    scanf("%d",&x);
    //x=jury(x);
    return x;
}

void answer(int x) {
    printf("! %d
",x);
}


int main() {
    n=11;

    inv[1]=1;
    for(int i=2;i<MOD;i++){
        inv[i]=(ll)inv[MOD%i]*(MOD-MOD/i)%MOD;
    }

    for(int i=0; i<=10; i++) {
        x[i+1]=i;
        int res=query(i);
        if(res==0) {
            answer(i);
            exit(0);
        }
        y[i+1]=res;
    }

    init_xi_xj();
    for(int i=11;i<MOD;i++){
        int res=lagrange(i);
        //cout<<"res="<<res<<endl;
        //cout<<"jury="<<jury(i)<<endl;
        if(res==0){
            answer(i);
            exit(0);
        }
    }
    answer(-1);
}

一般我们代入的点都是连续自然数，当 (x_i) 的取值是连续自然数时，插值法因为分母变得容易处理而变成n复杂度。

以从0开始取举例。

假设当前n为4，从1开始取。其实分子无论怎么取值都是可以用前缀积和后缀积预处理出来的，然后两段乘在一起。关键是分母，当分母是连续自然数的时候不必两两枚举。

(y_1frac{1*(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)}+y_2frac{(x-1)*(x-3)(x-4)}{(2-1)(2-3)(2-4)}+y_3frac{(x-1)(x-2)*(x-4)}{(3-1)(3-2)(3-4)}+y_4frac{(x-1)(x-2)(x-3)*1}{(4-1)(4-2)(4-3)})

分母是：(1*(-1)(-2)(-3),(1)*(-1)*(-2),(2)(1)*(-1),(3)(2)(1))

也就是前面阶乘，后面是负数的阶乘，负数的阶乘也很好想，比如这里n=4，那么第一项有3个负数，阶乘为负数，第二项为正数，第三项为负数……也就是n-i为奇数的项是负数。

可不断插点的重心拉格朗日插值法。

原式
$f(x)=sumlimits_{i=1}^{n} y_i prodlimits_{i e j}frac{x-x_j}{x_i-x_j} $

把分子补上缺的项，提出来，记 (g=prodlimits_{i=1}^{n} x-x_i)
$f(x)=gsumlimits_{i=1}^{n} frac{1}{x-x_i} prodlimits_{i e j}frac{y_i}{x_i-x_j} $

再记 (h_i= prodlimits_{i e j}frac{y_i}{x_i-x_j})
$f(x)=gsumlimits_{i=1}^{n} frac{h_i}{x-x_i} $

那么每次多加一个点的时候，更新 (g) 以及每一个 (h_i) 就可以了。

---

经典例题：求 $sumlimits_{i=1}^{n} i^k $ n很大但k只有几千万（需要连续前缀优化）
https://www.luogu.org/problemnew/show/P4593
推出计算式之后，就是计算上面那个东西。注意步骤，首先要确定多项式的系数，一个k次多项式的和一般都是k+1次的，这里是m+1次多项式的和，也就是m+2次多项式，需要m+3个点。

注意用拉格朗日插值时，判断k与x_i不重复再插，不然返回y_i就可以了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int MOD=1e9+7;

namespace Lagrange_Interpolation_Polynomial {
    const int MAXN=80;
    int n,x[MAXN+5],y[MAXN+5];
    int invd[MAXN+5][MAXN+5];//x_i-x_j的逆元
    int kd[MAXN+5];//k-x_i的值
    int pkd,pinvd[MAXN+5];//k-x的积，invd[i]的积
    //不需要计算k-x_i的逆元，因为复杂度其实一样

    int qpow(int x,int n)
    {
        int ans=1;
        while(n)
        {
            if(n&1)
                ans=((ll)ans*x)%MOD;
            x=((ll)x*x)%MOD;
            n>>=1;
        }
        return ans;
    }

    void init_xi_xj()
    {
        for(int i=1; i<=n; i++)
        {
            //printf("x[%d]=%d y[%d]=%d
",i,x[i],i,y[i]);
            for(int j=1; j<=n; j++)
            {
                if(i==j)
                    continue;
                else
                {
                    int d=(x[i]-x[j]);
                    if(d<0)
                        d+=MOD;
                    invd[i][j]=qpow(d,MOD-2);
                }
            }
        }
        for(int i=1;i<=n;i++){
            pinvd[i]=1;
            for(int j=1;j<=n;j++){
                if(i==j)
                    continue;
                pinvd[i]=(ll)pinvd[i]*invd[i][j]%MOD;
            }
            //cout<<"pinvd["<<i<<"]="<<pinvd[i]<<endl;
        }
    }

    void init_k(int k)
    {
        pkd=1;
        for(int i=1; i<=n; i++)
        {
            kd[i]=(k-x[i]);
            if(kd[i]<0)
                kd[i]+=MOD;
            pkd=(ll)pkd*kd[i]%MOD;
        }
        //cout<<"pkd="<<pkd<<endl;
    }

    inline int prod(int i,int k)
    {
        ll ans=(ll)pkd*qpow(kd[i],MOD-2)%MOD;
        ans=ans*(pinvd[i])%MOD;

        //cout<<"in prod ans="<<ans<<endl;
        return ans;
    }

    int lagrange(int k)
    {
        for(int i=1;i<=n;i++){
            if(x[i]==k)
                return y[i];
        }
        ll ans=0;
        //cout<<"k="<<n<<endl;
        init_k(k);

        for(int i=1; i<=n; i++)
        {
            ans=(ans+(ll)y[i]*(prod(i,k)))%MOD;
            //cout<<"in lag ans="<<ans<<endl;
        }
        //cout<<"(1~n)i^k="<<ans<<endl;
        return ans;
    }

}

using namespace Lagrange_Interpolation_Polynomial;



ll init(int pn,int k){
    //拉格朗日插值的n需要的是多项式的点的数目
    n=pn;

    for(int i=1;i<=n;i++){
        x[i]=i;
        y[i]=0;
        for(int j=1;j<=i;j++){
            y[i]=((ll)y[i]+qpow(j,k))%MOD;
        }
    }
    //cerr<<"n="<<n<<endl;;
    init_xi_xj();
}

ll calc(ll r){
    //cerr<<"lag"<<endl;
    ll t=lagrange(r);
    //cout<<"calc"<<" r="<<t<<endl;
    return t;
}

vector<int> vm;


int main() {
    int t;
    scanf("%d",&t);
    while(t--){
        ll n;
        int m;
        scanf("%lld%d",&n,&m);

        int k=m+1;
        //m+1次方的和，m+2次多项式，需要m+3个点
        vm.clear();
        for(int i=0;i<m;i++){
            int tt;
            scanf("%d",&tt);
            vm.push_back(tt);
        }
        sort(vm.begin(),vm.end());
        n%=MOD;
        //m+1次方的和

        //m+2次多项式，需要m+3个点
        init(k+2,k);
        ll ans=0;
        while(n){
            ans=(ans+calc(n))%MOD;
            for(auto &mi:vm){
                ans=(ans-qpow(mi,k))%MOD;
                if(ans<0)
                    ans+=MOD;
            }
            //cout<<"ans=!!!   "<<ans<<endl;
            if(vm.size()==0)
                break;
            vector<int> tv;
            for(int i=1;i<vm.size();i++){
                tv.push_back(vm[i]-vm[0]);
            }

            n-=vm[0];
            vm=tv;
            //vm.erase(ve.begin());
        }
        cout<<ans%MOD<<endl;
    }
}