luogu P5667 拉格朗日插值2 拉格朗日插值 多项式多点求值 NTT
LINK：P5667 拉格朗日插值2
给出了n个连续的取值的自变量的点值 求 f(m+1),f(m+2),...f(m+n).
如果我们直接把f这个函数给插值出来就变成了了多项式多点求值 这个难度好像有点大。
不妨 直接考虑拉格朗日插值。
设此时要求f(k) 那么则有 \(\sum_{i=0}^nf(i)\frac{\Pi_{i\neq j}(k-j)}{\Pi_{i\neq j} (i-j)}\)
可以化简一下 \(f(k)=\sum_{i=0}^nf(i)\frac{ \Pi_{i\neq j}(k-j) \cdot (-1)^{k-i} }{fac_i\cdot fac_{n-i}}\)
其实还是可以化简的 \(f(k)=\sum_{i=0}^nf(i)\frac{ \Pi (k-j) \cdot (-1)^{k-i} }{fac_i\cdot fac_{n-i}}\cdot \frac{1}{k-i}\)
\(f(k)=\sum_{i=0}^nf(i)\frac{k!\cdot (-1)^{k-i} }{fac_i\cdot fac_{n-i}}\cdot \frac{1}{(k-i)\cdot (k-n-1)!}\)
再提出来一些项 \(f(k)=k!\cdot \frac{1}{(k-n-1)!}\sum_{i=0}^nf(i)\frac{(-1)^{k-i} }{fac_i\cdot fac_{n-i}}\cdot \frac{1}{(k-i)}\)
容易发现这类似于卷积的形式。
不过 如果直接做卷积会出现问题 因为卷积的时候i<=j 而并非i<=n.
需要解决这个问题 比较容易想到 把j放到靠后的位置就能得到贡献了。
那么其实就是 所以维护一个长度为2n的序列进行卷积即可。
这个时候 最好的方法就是 把两个多项式写下来 看一下卷积的过程。
分析一下每一个位置都应该是什么数字即可。
注意 卡常的话中间有一个求2n长度逆元的东西 可以采用前缀积后缀积的方法来优化成O(n).
但是 由于多项式长度2^20 还是跑的很慢.. 勉强卡过。
码力还行 出错的地方是 数组开小了 开了1e6 少开了一点 wa了3,4发 我tcl.
const int MAXN=1100010,G=3;
int lim,n,m,N;
int rev[MAXN];
ll f[MAXN],b[MAXN],a[MAXN];
ll fac[MAXN],pre[MAXN],suf[MAXN],inv[MAXN];
inline int ksm(ll b,int p)
{
ll cnt=1;
while(p)
{
if(p&1)cnt=cnt*b%mod;
p=p>>1;b=b*b%mod;
}
return cnt;
}
inline void NTT(ll *a,int op)
{
rep(0,lim-1,i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int len=2;len<=lim;len=len<<1)
{
int mid=len>>1;
int wn=ksm(G,op==1?(mod-1)/len:mod-1-(mod-1)/len);
for(int j=0;j<lim;j+=len)
{
ll d=1;
for(int i=0;i<mid;++i)
{
ll x=a[i+j],y=a[i+j+mid]*d%mod;
a[i+j]=(x+y)%mod;a[i+j+mid]=(x-y+mod)%mod;
d=d*wn%mod;
}
}
}
if(op==-1)
{
int IN=ksm(lim,mod-2);
rep(0,lim-1,i)a[i]=a[i]*IN%mod;
}
}
int main()
{
//freopen("1.in","r",stdin);
get(n);get(m);fac[0]=1;
rep(0,n,i)get(f[i]);
rep(1,n,i)fac[i]=fac[i-1]*i%mod;
inv[n]=ksm(fac[n],mod-2);
fep(n-1,0,i)inv[i]=inv[i+1]*(i+1)%mod;
rep(0,n,i)a[i]=(f[i]*(((n-i)&1)?-1:1)%mod*inv[i]%mod*inv[n-i]%mod+mod)%mod;
pre[0]=1;N=n<<1|1;suf[N+1]=1;
rep(1,N,i)pre[i]=pre[i-1]*(m-n+i-1)%mod;
fep(N,1,i)suf[i]=suf[i+1]*(m+i-n-1)%mod;
ll IN=ksm(pre[N],mod-2);
rep(1,N,i)b[i]=pre[i-1]*suf[i+1]%mod*IN%mod;
lim=1;while(lim<N+N)lim=lim<<1;
rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
NTT(a,1);NTT(b,1);
rep(0,lim-1,i)a[i]=a[i]*b[i]%mod;
ll w=1;rep(m-n,m,i)w=w*i%mod;
b[n+1]=w;rep(n+2,N,i)b[i]=b[i-1]*(m+i-n-1)%mod*IN%mod*pre[i-n-2]%mod*suf[i-n]%mod;
NTT(a,-1);rep(n+1,N,i)printf("%lld ",a[i]*b[i]%mod);return 0;
}
中间的地方看起来确实比较ex... 这也是没有办法的事情。
