zoukankan      html  css  js  c++  java
  • 【知识点】多项式乘法

    FFT:

    没啥好说的吧。。

    证明应该都会,写的时候记住两个点就行:

    1.怎么定义复数?千万别写成

    complex<double> w=(1,0);

    可以自己试一下这样输出什么东西……

    2.枚举len,遍历前一半,用原来的$a_{i},a_{i+len/2}$值计算新的$a_{i},a_{i+len/2}$值。

    我觉得我还是解释一下这么做的原因好一点。。

    我们希望得到多项式$A(x)$在$w_{n}^{1},w_{n}^{2},cdots,w_{n}^{n}$处的点值表示。

    而我们有$w_{n}^{k}=w_{frac{n}{2}}^{frac{k}{2}},w_{n}^{k+frac{n}{2}}=-w_{frac{n}{2}}^{frac{k}{2}}$。

    那么把偶数下标的系数和奇数下标的系数分别拆出来组成两个多项式$A_{0}(x)$和$A_{1}(x)$,则有

    $A(w_{n}^{k})=A_{0}(w_{n}^{2k})+w_{n}^{k}A_{1}(w_{n}^{2k})=A_{0}(w_{frac{n}{2}}^{k})+w_{n}^{k}A_{1}(w_{frac{n}{2}}^{k})$

    $A(w_{n}^{k+frac{n}{2}})=A_{0}(w_{n}^{2k})-w_{n}^{k}A_{1}(w_{n}^{2k})=A_{0}(w_{frac{n}{2}}^{k})-w_{n}^{k}A_{1}(w_{frac{n}{2}}^{k})$

    那么我们可以先递归处理$A_{0}(x)$和$A_{1}(x)$,然后爆枚$k$,得到$A(x)$的$n$个点值。

    这东西很像线段树,倒数第$i$层的$n$是$2^{i}$,有$frac{n}{2^{i}}$个多项式需要爆枚$2^{i}$次。

    每层爆枚的复杂度是$O(n)$,一共$log{n}$层,总复杂度$O(nlog{n})$。

    发现位置$i$最终会被换到$i$的二进制位反过来的位置,于是可以不用递归,直接模拟线段树上爆枚的过程即可。

    强烈推荐伟大的石神给出的里程碑式的证明

    代码:

    #include<bits/stdc++.h>
    #define maxn 3000005
    #define maxm 500005
    #define inf 0x7fffffff
    #define ll long long
    #define cp complex<double>
    #define pi acos(-1)
    
    using namespace std;
    cp a[maxn],b[maxn],c[maxn];
    int ind[maxn];
    
    inline int read(){
        int x=0,f=1; char c=getchar();
        for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
        for(;isdigit(c);c=getchar()) x=x*10+c-'0';
        return x*f;
    }
    
    inline void fft(cp *t,int n,int op){
        for(int i=0;i<n;i++) ind[i]=(i&1)?((ind[i>>1]>>1)|(n>>1)):(ind[i>>1]>>1);
        for(int i=0;i<n;i++) if(ind[i]<i) swap(t[i],t[ind[i]]);
        for(int len=2;len<=n;len<<=1)
            for(int j=0;j<n;j+=len){
                cp w(1,0),p(cos(op*2*pi/len),sin(op*2*pi/len));
                for(int k=0;k<(len>>1);k++,w*=p){
                    cp t1=t[j+k+(len>>1)],t2=t[j+k];
                    t[j+k+(len>>1)]=t2-t1*w,t[j+k]=t2+t1*w;
                }
            }
        return;
    }
    
    int main(){
        int n=read(),m=read();
        for(int i=0;i<=n;i++) a[i]=read();
        for(int i=0;i<=m;i++) b[i]=read();
        int len=1; while(len<=n+m) len<<=1;
        fft(a,len,1),fft(b,len,1);
        for(int i=0;i<len;i++) c[i]=a[i]*b[i];
        fft(c,len,-1);
        for(int i=0;i<=n+m;i++) cout<<(int)(c[i].real()/len+0.5)<<" ";
        cout<<endl;
        return 0;
    }
    FFT

    NTT:

    原根的性质与单位根一模一样。只是$IDFT$的时候得求个逆元。

    题里有模数就用$NTT$,没有模数也尽量用$NTT$,避免出现时间精度双爆炸的惨剧。

    注意$w_{n}^{k}=g^{frac{mod-1}{n}cdot k}$,证明时将k分离出来考虑即可。

    我还是证一下吧。。虽然这也不算是证明了。。

    根据原根的性质,有$g_{i}mod{p}$的值互不相同$(0leq i<p)$。

    那么实际上$w_{n}^{k}=g^{frac{k}{n}}$,都满足首尾值相等且中间值互不相同。

    但是等一下,这个$frac{k}{n}$不是整数啊?

    而且指数没有同余定律,只有费马小定理$a^{p-1}equiv apmod p$。

    那我们只能取$n|(mod-1)$的$n$来做了,此时$g^{frac{k}{n}}=g^{frac{k(mod-1)}{n}}$。

    现在我们就理解为什么$NTT$只能用形如$2^{k}+1$的模数了。

    伟大的石神在刚才那篇里也更了$NTT$!

    代码:

    #include<bits/stdc++.h>
    #define maxn 3000005
    #define maxm 500005
    #define inf 0x7fffffff
    #define ll long long
    #define mod 998244353
    #define g 3
    
    using namespace std;
    ll a[maxn],b[maxn],res[maxn],ind[maxn];
    
    inline ll read(){
        ll x=0,f=1; char c=getchar();
        for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
        for(;isdigit(c);c=getchar()) x=x*10+c-'0';
        return x*f;
    }
    
    inline ll power(ll a,ll b){ll ans=1;while(b)ans=(b&1)?ans*a%mod:ans,a=a*a%mod,b>>=1;return ans;}
    inline void ntt(ll *t,ll n,ll op){
        for(ll i=0;i<n;i++) ind[i]=(i&1)?((ind[i>>1]>>1)|(n>>1)):(ind[i>>1]>>1);
        for(ll i=0;i<n;i++) if(ind[i]<i) swap(t[i],t[ind[i]]); 
        for(ll len=1;len<=n;len<<=1){
            ll p=power(g,(mod-1)/len);
            if(op==-1) p=power(p,mod-2);
            for(ll i=0;i<n;i+=len)
                for(ll j=i,w=1,tp;j<i+(len>>1);j++,w=w*p%mod)
                    tp=t[j+(len>>1)]*w%mod,t[j+(len>>1)]=(t[j]-tp+mod)%mod,t[j]=(t[j]+tp)%mod;
        }
        return;
    }
    
    int main(){
        ll n=read(),m=read();
        for(ll i=0;i<=n;i++) a[i]=read();
        for(ll i=0;i<=m;i++) b[i]=read();
        ll len=1; while(len<=n+m) len<<=1;
        ntt(a,len,1),ntt(b,len,1);
        for(ll i=0;i<len;i++) res[i]=a[i]*b[i];
        ntt(res,len,-1);
        ll mo=power(len,mod-2);
        for(ll i=0;i<=n+m;i++) cout<<res[i]*mo%mod<<" ";
        cout<<endl;
        return 0;
    }
    NTT

    (upd:数组一定要开两倍,点数少了就没法表示一个多项式)

  • 相关阅读:
    CSS之旅——第二站 如何更深入的理解各种选择器
    CSS之旅——第一站 为什么要用CSS
    记录一些在用wcf的过程中走过的泥巴路 【第一篇】
    asp.net mvc 之旅—— 第二站 窥探Controller下的各种Result
    asp.net mvc 之旅—— 第一站 从简单的razor入手
    Sql Server之旅——终点站 nolock引发的三级事件的一些思考
    Sql Server之旅——第十四站 深入的探讨锁机制
    Sql Server之旅——第十三站 对锁的初步认识
    Sql Server之旅——第十二站 sqltext的参数化处理
    Sql Server之旅——第十一站 简单说说sqlserver的执行计划
  • 原文地址:https://www.cnblogs.com/YSFAC/p/11738259.html
Copyright © 2011-2022 走看看