首页 > 其他 > 详细

多项式入门

时间:2019-11-18 21:05:31      阅读:86      评论:0      收藏:0      [点我收藏+]

FFT/NTT

这里令 \(A(x)\) 表示多项式 \(A\)\(x\) 处的取值,\(A[x]\) 表示多项式 \(A\) 的第 \(x\) 项.

卷积形如:

\(C[r]=\sum_{p,q}[(p+q)\mod n=r]A[p]\times B[q]\)

\(w^{i}\)\(1\)\(i\) 次单位根.

则有 \(\frac{1}{n}\sum_{k=0}^{n-1}w^{vk}=[v\mod n=0]\)

而上面 \(C[r]=\sum_{p,q}[(p+q)\mod n=r]A[p]\times B[q]\) 中有一个特判:\([(p+q)\mod n=r]\)

我们考虑将这个特判换掉:

\(\Rightarrow\) \([(p+q-r)\mod n=0]\)

\(\Rightarrow \frac{1}{n}\sum_{k=0}^{n-1}w^{(p+q-r)k}\)

展开,得 \(\frac{1}{n}\sum_{k=0}^{n-1}w^{-rk}w^{pk}w^{qk}\)

那么,\(C[r]=\sum_{p,q}\frac{1}{n}\sum_{k=0}^{n-1}w^{-rk}w^{pk}w^{qk} A[p]\times B[q]\)

整理得好看一点,可得 \(C[r]=\frac{1}{n}\sum_{k=0}^{n-1}(w^{-r})^k\sum_{p=0}^{n-1}(w^k)^pA[p]\sum_{q=0}^{n-1}(w^k)^qB[q]\)

\(C[r]=\frac{1}{n}\sum_{k=0}^{n-1}(w^{-r})^k A(w^k)B(w^k)\)

构造多项式 \(G\), 满足 \(G[k]=A(w^k)B(w^k)\) ,然后求出 \(\frac{G(w^{-r})}{n}\) 就是 \(C[r]\) 了.

注意:由于平时我们写多项式的时候初始的 \(n\) 未必是 2 的整数次幂,所以不能完成循环卷积.

第二类斯特林数

\(S(n,m)\) 表示将 \(n\) 个互不相同的元素装进 \(m\) 个相同盒子,且盒子不能为空的方案数.

朴素递推:\(S(n,m)=S(n-1,m-1)+S(n-1,m)\times m\)

因为盒子不为空,所以第一种情况只能放入最新的盒子,而第二种情况的话 \(m\) 个盒子可以随便放.

这个太慢了,考虑二项式反演:

假设 \(n\) 固定,我们要求 \(S(n,m)\) ,则令 \(f[k]\) 表示钦定 \(k\) 个盒子为空,其余随便选的方案,\(g[k]\) 表示恰好.

则有 \(f[k]=\binom{m}{k}(m-k)^n=\sum_{i=k}^{m}\binom{i}{k}g[i]\)

反演,得 \(g[0]=\sum_{i=0}^{m}(-1)^if[i]\).

注意:这里实际上默认将不同得盒子看作不同,所以还要除一个阶乘.

那么,就有 \(S(n,m)=\frac{1}{m!}\sum_{i=0}^{m}(-1)^i\binom{m}{i}(m-i)^n\)

我们把组合数展开,然后约掉最外面的 \(\frac{1}{m!}\) 就会得到:\(S(n,m)=\sum_{i=0}^{m}\frac{(-1)^i}{i!}\frac{(m-i)^n}{(m-i)!}\)

这个用 NTT 加速多项式乘法来算即可.

#include <bits/stdc++.h> 
#define ll long long
#define setIO(s) freopen(s".in","r",stdin) 
using namespace std;         
int n;          
const ll mod=167772161,G=3,N=400006;       
ll f[N<<1],g[N<<1],fac[N],inv[N];     
ll qpow(ll x,ll y)
{     
    ll tmp=1ll;
    while(y)
    {
        if(y&1) tmp=tmp*x%mod;        
        y>>=1,x=x*x%mod;
    } 
    return tmp;   
}
void NTT(ll *a,int len,int flag)
{
    int i,j,k,mid;
    for(i=k=0;i<len;++i)
    {
        if(i>k) swap(a[i],a[k]);
        for(j=len>>1;(k^=j)<j;j>>=1);       
    }
    for(mid=1;mid<len;mid<<=1)         
    {
        ll wn=qpow(G,(mod-1)/(mid<<1));   
        if(flag==-1) wn=qpow(wn,mod-2);    
        for(i=0;i<len;i+=(mid<<1))        
        { 
            ll w=1ll;  
            for(j=0;j<mid;++j)
            {
                ll x=a[i+j],y=w*a[i+mid+j]%mod;       
                a[i+j]=(x+y)%mod,a[i+j+mid]=(x-y+mod)%mod;  
                w=w*wn%mod; 
            }
        }
    }
    if(flag==-1)
    {
        ll re=qpow(len,mod-2);        
        for(i=0;i<len;++i) a[i]=a[i]*re%mod; 
    }
}              
int main()
{
    // setIO("input");
    scanf("%d",&n);                    
    fac[0]=1ll;
    inv[0]=1ll;        
    int i,j,limit=1;
    for(i=1;i<=n;++i)   fac[i]=fac[i-1]*1ll*i%mod,  inv[i]=qpow(fac[i],mod-2);   
    for(i=0;i<=n;++i)  
    {
        g[i]=qpow(i,n)*inv[i]%mod;               
        if(i&1)  f[i]=mod-inv[i];           
        else f[i]=inv[i];     
    }  
    for(;limit<=2*(n+1);limit<<=1);                                        
    NTT(f,limit,1),NTT(g,limit,1);                                                                                     
    for(i=0;i<limit;++i)  f[i]=f[i]*g[i]%mod;          
    NTT(f,limit,-1);                                              
    for(i=0;i<=n;++i)      printf("%lld ",f[i]);   
    return 0;  
}   

多项式入门

原文:https://www.cnblogs.com/guangheli/p/11885107.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!