码迷,mamicode.com
首页 > 其他好文 > 详细

多项式FFT相关模板

时间:2016-07-14 23:59:24      阅读:530      评论:0      收藏:0      [点我收藏+]

标签:

自己码了一个模板...有点辛苦...常数十分大,小心使用

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <algorithm>
#include <vector>
using namespace std;
#define ll long long
#define pb push_back
ll MOD=998244353;
#define SZ 666666
ll w[2][SZ];
ll qp(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%MOD;
        a=a*a%MOD; b>>=1;
    }
    return ans;
}
int K;
void fftinit(int n)
{
    for(K=1;K<n;K<<=1);
    w[0][0]=w[0][K]=1;
    ll g=qp(3,(MOD-1)/K); //3是原根
    for(int i=1;i<K;i++) w[0][i]=w[0][i-1]*g%MOD;
    for(int i=0;i<=K;i++) w[1][i]=w[0][K-i];
}
void fft(int* x,int v)
{
    for(int i=0,j=0;i<K;i++)
    {
        if(i>j) {x[i]^=x[j]; x[j]^=x[i]; x[i]^=x[j];}
        for(int l=K>>1;(j^=l)<l;l>>=1);
    }
    for(int i=2;i<=K;i<<=1)
    {
        for(int j=0;j<K;j+=i)
        {
            for(int l=0;l<i>>1;l++)
            {
                ll t=(ll)x[j+l+(i>>1)]*w[v][K/i*l]%MOD;
                x[j+l+(i>>1)]=(x[j+l]-t+MOD)%MOD;
                x[j+l]=(x[j+l]+t)%MOD;
            }
        }
    }
    if(!v) return;
    ll rv=qp(K,MOD-2);
    for(int i=0;i<K;i++) x[i]=x[i]*rv%MOD;
}
struct poly
{
    vector<int> ps;
    int cs() {return ps.size()-1;}
    int& operator [] (int x) {return ps[x];} //ps.at(x)
    void sc(int x) {ps.resize(x+1);}
    void dbg()
    {
        bool fi=0;
        for(int i=cs();i>=0;i--)
        {
            if(!ps[i]) continue;
            if(fi)
            {
                if(i==0) printf("+%d",ps[i]);
                else if(ps[i]==1) printf("+");
                else if(ps[i]==-1) printf("-");
                else printf("+%d",ps[i]);
            }
            else
            {
                if(i==0) printf("%d",ps[i]);
                else if(ps[i]==1);
                else if(ps[i]==-1) printf("-");
                else printf("%d",ps[i]);
            }
            if(i>1) printf("x^%d",i);
            else if(i==1) printf("x");
            fi=1;
        }
        if(!fi) printf("0");
        putchar(10);
    }
    void clr()
    {
        int p=cs()+1;
        while(p&&!ps[p-1]) --p;
        sc(p-1);
    }
};
ll gm(ll x)
{
    x=x%MOD;
    if(x<0) x+=MOD;
    return x;
}
namespace PolyMul{int ta[SZ],tb[SZ],tc[SZ];}
poly operator * (poly a,poly b)
{
    using namespace PolyMul;
    if(a.cs()<200||b.cs()<200)
    {
        poly g;
        g.sc(a.cs()+b.cs());
        for(int i=0;i<=a.cs();i++)
        {
            for(int j=0;j<=b.cs();j++) g[i+j]=gm(g[i+j]+a[i]*(ll)b[j]%MOD);
        }
        return g;
    }
    poly c;
    int t=a.cs()+b.cs();
    c.sc(t); fftinit(t+1);
    memset(ta,0,sizeof(int)*K);
    memset(tb,0,sizeof(int)*K);
    memset(tc,0,sizeof(int)*K);
    for(int i=a.cs();i>=0;i--) ta[i]=a[i];
    for(int i=b.cs();i>=0;i--) tb[i]=b[i];
    fft(ta,0); fft(tb,0);
    for(int i=0;i<K;i++) tc[i]=(ll)ta[i]*tb[i]%MOD;
    fft(tc,1);
    for(int i=t;i>=0;i--) c[i]=tc[i];
    c.clr();
    return c;
}
namespace PolyInv{int ay[SZ],a0[SZ],tmp[SZ];}
void ginv(int t)
{
    using namespace PolyInv;
    if(t==1) {a0[0]=qp(ay[0],MOD-2); return;}
    ginv((t+1)>>1); fftinit(t+t+3);
    memset(tmp,0,sizeof(int)*K);
    for(int i=t;i<K;i++) tmp[i]=a0[i]=0; 
    for(int i=0;i<t;i++) tmp[i]=ay[i];
    fft(tmp,0); fft(a0,0);
    for(int i=0;i<K;i++) a0[i]=gm((2-(ll)tmp[i]*a0[i])%MOD*a0[i]);
    fft(a0,1);
    for(int i=t;i<K;i++) a0[i]=0;
}
poly inv(poly x)
{
    using namespace PolyInv;
    poly y; y.sc(x.cs());
    for(int i=x.cs();i>=0;i--) ay[i]=x[i];
    ginv(x.cs()+1);
    for(int i=x.cs();i>=0;i--) y[i]=a0[i];
    y.clr();
    return y;
}
poly operator + (poly a,poly b)
{
    poly w; w.sc(max(a.cs(),b.cs()));
    for(int i=a.cs();i>=0;i--) w[i]=a[i];
    for(int i=b.cs();i>=0;i--) w[i]+=b[i], w[i]=gm(w[i]);
    return w;
}
poly operator - (poly a,poly b)
{
    poly w; w.sc(max(a.cs(),b.cs()));
    for(int i=a.cs();i>=0;i--) w[i]=a[i];
    for(int i=b.cs();i>=0;i--) w[i]-=b[i], w[i]=gm(w[i]);
    w.clr();
    return w;
}
void div(poly a,poly b,poly& d,poly& r)
{
    int n=a.cs(),m=b.cs();
    if(n<m) {d.sc(0); d[0]=0; r=a; return;}
    fftinit(2*n);
    poly aa=a; reverse(aa.ps.begin(),aa.ps.end());
    poly bb=b; reverse(bb.ps.begin(),bb.ps.end());
    bb.sc(n-m); bb=inv(bb); d=aa*bb; d.sc(n-m);
    reverse(d.ps.begin(),d.ps.end()); r=a-b*d;
    r.clr();
}
poly operator / (poly a,poly b)
{
    poly d,r; div(a,b,d,r); return d;
}
poly operator % (poly a,poly b)
{
    poly d,r; div(a,b,d,r); return r;
}
poly dev(poly x)
{
    for(int i=1;i<=x.cs();i++) x[i-1]=(ll)x[i]*i%MOD;
    x.sc(x.cs()-1); return x;
}
poly inte(poly x) //C=0
{
    x.sc(x.cs()+1);
    for(int i=x.cs();i>=1;i--) x[i]=x[i-1]; x[0]=0;
    for(int i=x.cs();i>=1;i--) x[i]=(ll)x[i]*qp(i,MOD-2)%MOD;
    return x;
}
ll qz_(poly& a,ll x)
{
    ll ans=0;
    for(int i=a.cs();i>=0;i--) ans=(ans*x%MOD+a[i])%MOD;
    return gm(ans);
}
namespace PolyGetv{int xs[SZ],anss[SZ];};
void gv(poly f,int m,int* x,int* ans)
{
    //f.clr();
    if(f.cs()<=7)
    {
        for(int i=0;i<=m;i++) ans[i]=qz_(f,x[i]);
        return;
    }
    poly m0,m1,tmp;
    m0.sc(0); m1.sc(0); tmp.sc(1);
    m0[0]=m1[0]=1; tmp[1]=1;
    int hf=m/2;
    for(int i=0;i<=hf;i++) tmp[0]=gm(-x[i]), m0=m0*tmp;
    for(int i=hf+1;i<=m;i++) tmp[0]=gm(-x[i]), m1=m1*tmp;
    gv(f%m0,hf,x,ans);
    gv(f%m1,m-hf,x+hf+1,ans+hf+1);
}
vector<int> getv(poly a,vector<int> x)
{
    using namespace PolyGetv;
    a.clr();
    if(!x.size()) return vector<int>();
    int m=x.size()-1;
    for(int i=0;i<=m;i++) xs[i]=x[i];
    gv(a,m,xs,anss);
    vector<int> ans; ans.resize(m+1);
    for(int i=0;i<=m;i++) ans[i]=anss[i];
    return ans;
}
int main()
{
}

加减乘逆元除取模求导积分多点求值...感觉够用了。

大部分运算没有用题目测试过...都是小数据/目测啥的...有问题求评论告知。

相关介绍请见picks博客及上一篇FFT入门。

http://picks.logdown.com/posts/177631-fast-fourier-transform

http://picks.logdown.com/posts/189620-inverse-element-of-polynomial

http://picks.logdown.com/posts/197262-polynomial-division

http://www.cnblogs.com/zzqsblog/p/5665654.html

多项式FFT相关模板

标签:

原文地址:http://www.cnblogs.com/zzqsblog/p/5672002.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!