https://blog.youkuaiyun.com/Lee_w_j__/article/details/81135539
https://blog.youkuaiyun.com/w571523631/article/details/81133769
#include<bits/stdc++.h>
#define LL long long
#define INF 0x3f3f3f3f
#define ll long long
#define ULL unsigned long long
using namespace std;
const int maxn = 1000+4;
const LL mod = 1e9+7;
LL a[maxn];
int n;
LL b[maxn];
/*
测试用例 (函数)x^3-2*x+7
a[]={7,6,11,28,63}
第0-4项和 115
最高次项3,代入项数4
*/
namespace polysum {
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
const int D=2010;
ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
ll powmod(ll a,ll b){ll res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}
//函数用途:给出数列的(d+1)项,其中d为最高次方项
//求出数列的第n项,数组下标从0开始
ll calcn(int d,ll *a,ll n) { // a[0].. a[d] a[n]
if (n<=d) return a[n];
p1[0]=p2[0]=1;
rep(i,0,d+1) {
ll t=(n-i+mod)%mod;
p1[i+1]=p1[i]*t%mod;
}
rep(i,0,d+1) {
ll t=(n-d+i+mod)%mod;
p2[i+1]=p2[i]*t%mod;
}
ll ans=0;
rep(i,0,d+1) {
ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
if ((d-i)&1) ans=(ans-t+mod)%mod;
else ans=(ans+t)%mod;
}
return ans;
}
void init(int M) {
f[0]=f[1]=g[0]=g[1]=1;
rep(i,2,M+5) f[i]=f[i-1]*i%mod;
g[M+4]=powmod(f[M+4],mod-2);
per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;
}
//函数用途:给出数列的(m+1)项,其中m为最高次方
//求出数列的前(n-1)项的和
ll polysum(ll m,ll *a,ll n) { // a[0].. a[m] \sum_{i=0}^{n-1} a[i]
ll b[D];
for(int i=0;i<=m;i++) b[i]=a[i];
b[m+1]=calcn(m,b,m+1);
rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
return calcn(m+1,b,n-1);
}
ll qpolysum(ll R,ll n,ll *a,ll m) { // a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
if (R==1) return polysum(n,a,m);
a[m+1]=calcn(m,a,m+1);
ll r=powmod(R,mod-2),p3=0,p4=0,c,ans;
h[0][0]=0;h[0][1]=1;
rep(i,1,m+2) {
h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
h[i][1]=h[i-1][1]*r%mod;
}
rep(i,0,m+2) {
ll t=g[i]*g[m+1-i]%mod;
if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
}
c=powmod(p4,mod-2)*(mod-p3)%mod;
rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
rep(i,0,m+2) C[i]=h[i][0];
ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
if (ans<0) ans+=mod;
return ans;
}
} // polysum::init();
LL pow2(LL a,LL b)
{
LL res = 1;
while(b){
if(b&1) res = res%mod*a%mod%mod;
a = a%mod*a%mod%mod;
b>>=1;
}
return res%mod;
}
int main()
{
LL ans;
LL now;
polysum::init(1000+32);
cin>>n;
for(int i=0;i<=n;i++)
{
cin>>a[i];
}
cout<<polysum::polysum(n,a,n+2)<<endl;
// while(~scanf("%d",&n))
// {
// for(int i = 1;i<=n;i++){
// scanf("%lld",&a[i]);
// }
// now = 1;
// ans = 0;
// sort(a+1,a+1+n);
// for(int i = 1;i<=n;i++){
// if(a[i]==a[i-1]) continue;
// b[0] = 0;
// for(int j = 1;j<=n-i+1;j++){
// b[j] = j%mod*(pow2(j,n-i+1)%mod-pow2(j-1,n-i+1)%mod+mod)%mod;
// }
// ans = (ans%mod+now*(polysum::polysum(n-i+1,b,a[i]+1)%mod-polysum::polysum(n-i+1,b,a[i-1]+1)%mod+mod)%mod)%mod;
// now = (now%mod*a[i]%mod)%mod;
// }
// printf("%lld\n",ans);
// }
}
例外一种n^2的更为灵活的多项式拉格朗日插值法
http://www.ebola.pro/special/Lagrange
#include<bits/stdc++.h>
#define ha 998244353
using namespace std;
int Pow(int a,int b)
{
int ans=1;
for(;b;b>>=1,a=1ll*a*a%ha)
if(b&1) ans=1ll*ans*a%ha;
return ans;
}
int lagrange(int *arrx,int *arry,int n,int x)
{
int res=0;
for(int i=0;i<n;i++)
{
int s1=1,s2=1;
for(int j=0;j<n;j++)
{
if(i==j) continue;
s1=1ll*s1*(x-arrx[j])%ha;
s2=1ll*s2*(arrx[i]-arrx[j])%ha;
}
res=(res+1ll*s1*Pow(s2,ha-2)%ha*arry[i])%ha;
}
return res;
}
int a[2010],b[2010];
int main()
{
int n,x;
scanf("%d",&n);
for(int i=0;i<=n;i++)
scanf("%d%d",a+i,b+i);
scanf("%d",&x);
printf("%d\n",lagrange(a,b,n+1,x));
return 0;
}