myh的超级多项式 - 分治FFT

题目大意:
给你 n , k , { v k } , { f k } n,k,\{v_k\},\{f_k\} n,k,{vk},{fk},并且已知 f n = ∑ i = 1 k a i v i n f_n=\sum_{i=1}^k a_iv_i^{n} fn=i=1kaivin,求 f n , n − k ≤ 1000 , k ≤ 1 0 5 , O ( k ( n − k ) ) f_n,n-k\le1000,k\le10^5,O(k(n-k)) fn,nk1000,k105,O(k(nk))可过。
题解:考虑构造 f n f_n fn关于 f n − k … f n − 1 f_{n-k}\ldots f_{n-1} fnkfn1的递推式,也就是要求系数序列 { c k } \{c_k\} {ck},使得: f n = ∑ i = 1 k c i f n − i f_n=\sum_{i=1}^kc_if_{n-i} fn=i=1kcifni成立。
c 0 = − 1 c_0=-1 c0=1,等价于:
∑ i = 0 k c i f n − i = 0 ∑ i = 0 k c i ∑ j = 1 k a j v j n − i = 0 \sum_{i=0}^kc_if_{n-i}=0\\ \sum_{i=0}^kc_i\sum_{j=1}^ka_jv_j^{n-i}=0 i=0kcifni=0i=0kcij=1kajvjni=0
其一个充分条件是:
∀ j ∈ [ 1 , k ] , ∑ i = 0 k c i a j v j n − i = 0 \forall j\in[1,k],\sum_{i=0}^kc_ia_jv_j^{n-i}=0 j[1,k],i=0kciajvjni=0
亦即:
∀ j ∈ [ 1 , k ] , ∑ i = 0 k c i v j k − i = 0 \forall j\in[1,k],\sum_{i=0}^kc_iv_j^{k-i}=0 j[1,k],i=0kcivjki=0
F ( x ) = ∑ i = 0 k c k − i x i F(x)=\sum_{i=0}^kc_{k-i}x^i F(x)=i=0kckixi,则上式等价于 { v k } \{v_k\} {vk} F ( x ) F(x) F(x)的零点,再由 [ x k ] F ( x ) = c 0 = − 1 [x^k]F(x)=c_0=-1 [xk]F(x)=c0=1可知: F ( x ) = − ∏ i = 1 k ( x − v i ) F(x)=-\prod_{i=1}^k(x-v_i) F(x)=i=1k(xvi)
使用分治FFT求出 F ( x ) F(x) F(x)进而得到 { c k } \{c_k\} {ck}进而求出 f n f_n fn即可。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<vector>
#define p 1004535809
#define N 100010
#define lint long long
#define gc getchar()
#define clr(a,n) memset(a,0,sizeof(int)*(n))
using namespace std;
inline int inn()
{
    int x,ch;while((ch=gc)<'0'||ch>'9');
    x=ch^'0';while((ch=gc)>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^'0');return x;
}
int a[N<<2],b[N<<2],r[N<<2],c[N],x[N],v[N];vector<int> f[N];
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%p) (k&1)?ans=(lint)ans*x%p:0;return ans; }
inline int NTT(int *a,int n,int s)
{
    for(int i=1;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=2;i<=n;i<<=1)
    {
        int wn=fast_pow(3,s>0?(p-1)/i:p-1-(p-1)/i);
        for(int j=0,t=i>>1;j<n;j+=i) for(int k=0,w=1,x,y;k<t;k++,w=(lint)w*wn%p)
            x=a[j+k],y=(lint)w*a[j+k+t]%p,((a[j+k]=x+y)>=p?a[j+k]-=p:0),((a[j+k+t]=x-y)<0?a[j+k+t]+=p:0);
    }
    if(s<0) { int ninv=fast_pow(n,p-2);for(int i=0;i<n;i++) a[i]=(lint)a[i]*ninv%p; }
    return 0;
}
inline int tms(const vector<int> &A,const vector<int> &B,vector<int> &C)
{
    int m1=(int)A.size(),m2=(int)B.size();
    int n=1,L=0;for(;n<m1+m2-1;n<<=1,L++);
    for(int i=1;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
    clr(a,n);for(int i=0;i<m1;i++) a[i]=A[i];NTT(a,n,1);
    clr(b,n);for(int i=0;i<m2;i++) b[i]=B[i];NTT(b,n,1);
    for(int i=0;i<n;i++) a[i]=(lint)a[i]*b[i]%p;NTT(a,n,-1);
    C.resize(m1+m2-1);for(int i=0;i<m1+m2-1;i++) C[i]=a[i];
    return 0;
}
int solve(int l,int r)
{
    if(l==r) return f[l].resize(2),f[l][0]=(v[l]?p-v[l]:0),f[l][1]=1;
    int mid=(l+r)>>1;solve(l,mid),solve(mid+1,r),tms(f[l],f[mid+1],f[l]);
    return 0;
}
int main()
{
    int n=inn(),k=inn();
    for(int i=1;i<=k;i++) v[i]=inn();
    for(int i=1;i<=k;i++) x[i]=inn();
    solve(1,k);
    for(int i=0;i<=k;i++) c[i]=(f[1][i]?p-f[1][i]:0);
    for(int i=0;i<=k/2;i++) swap(c[i],c[k-i]);
    for(int i=k+1;i<=n;i++)
    {
        lint w=0;
        for(int j=1;j<=k;j++) w+=(lint)c[j]*x[i-j]%p;
        x[i]=(int)(w%p);
    }
    return !printf("%d\n",x[n]);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值