模板记录代码模板：FFT:NTT:分治NTT:多项式求逆:多项式开根:扩展lucas:多项式大全

71 阅读 0 评论 47 点赞

我是靠谱客的博主敏感龙猫，最近开发中收集的这篇文章主要介绍模板记录代码模板：FFT:NTT:分治NTT:多项式求逆:多项式开根:扩展lucas:多项式大全，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

代码模板：

#include<bits/stdc++.h>
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define sz(a) int(a.size())
#define clr(a) memset(a,0,sizeof(a))
#define all(a) a.begin(),a.end()
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,int> pli;
const int inf=1e9;
const ll Inf=1e18;
const int mod=0;
int gi() {
    int x=0,o=1;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
    if(ch=='-') o=-1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    return x*o;
}
template<typename T> bool chkmax(T &a,T b) { return a<b?a=b,1:0; };
template<typename T> bool chkmin(T &a,T b) { return a>b?a=b,1:0; };
int add(int a,int b) { return a+b>=mod?a+b-mod:a+b; }
int sub(int a,int b) { return a-b<0?a-b+mod:a-b; }
void inc(int &a,int b) { a=(a+b>=mod?a+b-mod:a+b); }
void dec(int &a,int b) { a=(a-b<0?a-b+mod:a-b); }
string to_string(string s) { return '"'+s+'"'; }
string to_string(const char *s) { return to_string((string)s); }
template<typename A,typename B> string to_string(pair<A,B> p) {
    return "("+to_string(p.fi)+","+to_string(p.se)+")";
}
template<typename T> string to_string(T v) {
    int fst=1;string ret="{";
    for(auto x:v) {
        if(!fst) ret+=",";
        fst=0,ret+=to_string(x);
    }
    ret+="}";return ret;
}
void dbg_out() { cerr<<endl; }
template<typename Head,typename... Tail> void dbg_out(Head H,Tail... T) {
    cerr<<" "<<to_string(H);
    dbg_out(T...);
}
#define dbg(...) cerr<<"{"<<#__VA_ARGS__<<"}:",dbg_out(__VA_ARGS__)
int main() {
#ifndef ONLINE_JUDGE
    freopen(".in","r",stdin);
    freopen(".out","w",stdout);
#endif
    
    return 0;
}

FFT:

const double pi=acos(-1);
typedef complex<double> E;
namespace FFT {
    int n,m,l,r[N];
    void FFT(E *a,int f) {
        for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
        for(int i=1;i<n;i<<=1) {
            E wn(cos(pi/i),f*sin(pi/i));
            for(int p=i<<1,j=0;j<n;j+=p) {
                E w(1,0);
                for(int k=0;k<i;k++,w*=wn) {
                    E x=a[j+k],y=w*a[j+k+i];
                    a[j+k]=x+y;a[j+k+i]=x-y;
                }
            }
        }
        if(f==-1) for(int i=0;i<=n;i++) a[i]/=n;
    }
    void calc(E *a,E *b,int _n,int _m) {
        n=_n,m=_m,l=0;
        m+=n;for(n=1;n<=m;n<<=1) l++;
        for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        FFT(a,1);FFT(b,1);
        for(int i=0;i<=n;i++) a[i]*=b[i];
        FFT(a,-1);
    }
}

NTT:

namespace NTT {
    int n,m,l,r[N];
    void NTT(int *a,int f) {
        for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
        for(int i=1;i<n;i<<=1) {
            int wn=q_pow(pr,phi/(i<<1));
            if(f==-1) wn=q_pow(wn,mod-2);
            for(int p=i<<1,j=0;j<n;j+=p) {
                int w=1;
                for(int k=0;k<i;k++,w=1ll*w*wn%mod) {
                    int x=a[j+k],y=1ll*w*a[j+k+i]%mod;
                    a[j+k]=(x+y)%mod;a[j+k+i]=(x-y+mod)%mod;
                }
            }
        }
        if(f==-1) {
            int inv=q_pow(n,mod-2);
            for(int i=0;i<=n;i++) a[i]=1ll*a[i]*inv%mod;
        }
    }
    int calc(int *a,int *b,int _n,int _m,int type) {//type==2是多项式求逆时用的
        n=_n,m=_m,l=0;
        m+=n;for(n=1;n<=m;n<<=1) l++;
        for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        NTT(a,1);NTT(b,1);
        for(int i=0;i<=n;i++) {
            if(type==1) a[i]=1ll*a[i]*b[i]%mod;
            else a[i]=1ll*a[i]*b[i]%mod*b[i]%mod;
        }
        NTT(a,-1);return n;
    }
}

分治NTT:

void cdq(int l,int r) {
    if(l==r) {
        f[l]=((1ll*q_pow(2,1ll*l*(l-1)/2)*q_pow(fac[l-1],mod-2)%mod-f[l])%mod+mod)%mod;//具体视题目推出的式子而定，这里是bzoj3456城市规划的代码
        return;
    }
    int mid=(l+r)>>1;
    cdq(l,mid);
    for(int i=0;i<=mid-l;i++) a[i]=f[i+l];
    for(int i=0;i<=r-l;i++) b[i]=g[i];
    int tmp=NTT::calc(a,b,mid-l,r-l);
    for(int i=mid+1-l;i<=r-l;i++) f[i+l]=(f[i+l]+a[i])%mod;
    for(int i=0;i<=tmp;i++) a[i]=b[i]=0;
    cdq(mid+1,r);
}

分治(NTT)可以在(O(nlog^2n))时间内求出形如(f[n]=(sum_{i=0}^{n-1}f[i]*a[n-i])+C)的式子的值，其中(C)是一个只与(n)有关的常数。具体做法是(CDQ)分治，每次用([l,mid])的(f[])值去更新([mid+1,r])的(f[])值，当递归到(l==r)时，(f[l])即被更新完毕，直接算出它最终的值即可。

多项式求逆:

void Inv(int *f,int *g,int len) {
    if(len==1) {
        g[0]=q_pow(f[0],mod-2);
        return;
    }
    Inv(f,g,len>>1);
    for(int i=0;i<len;i++) a[i]=f[i],b[i]=g[i];
    int tmp=NTT::calc(a,b,len-1,len-1,2);
    for(int i=0;i<len;i++) g[i]=((-a[i]+2*g[i])%mod+mod)%mod;
    for(int i=0;i<=tmp;i++) a[i]=b[i]=0;
}

多项式开根:

不会，留坑

扩展lucas:

int qpow(int a,int b,int mod) {
  int ret=1;
  while(b) {
    if(b&1) ret=1ll*ret*a%mod; a=1ll*a*a%mod; b>>=1;
  }
  return ret;
}
void exgcd(int a,int b,int &x,int &y) {
  if(!b) { x=1,y=0;return; }
  exgcd(b,a%b,x,y);
  int tmp=x;x=y;y=tmp-a/b*y;
}
int inv(int a,int p) {
  int x,y;exgcd(a,p,x,y);
  return (x+p)%p;
}
int fac(int n,int p,int pk) {
  if(n==0) return 1;
  int ret=1;
  for(int i=2;i<=pk;i++) if(i%p) ret=1ll*ret*i%pk;
  ret=qpow(ret,n/pk,pk);
  for(int i=2;i<=n%pk;i++) if(i%p) ret=1ll*ret*i%pk;
  return 1ll*ret*fac(n/p,p,pk)%pk;
}
int C(int n,int m,int p,int pk) {
  if(n<m) return 0;
  int a=fac(n,p,pk),b=fac(m,p,pk),c=fac(n-m,p,pk),k=0;
  for(int i=n;i;i/=p) k+=i/p;
  for(int i=m;i;i/=p) k-=i/p;
  for(int i=n-m;i;i/=p) k-=i/p;
  int ret=1ll*a*inv(b,pk)%pk*inv(c,pk)%pk*qpow(p,k,pk)%pk;
  return 1ll*ret*(mod/pk)%mod*inv(mod/pk,pk)%mod;
}
int lucas(int n,int m) {
  int ret=0,p=mod;
  for(int i=2;i<=p;i++) {
    if(p%i==0) {
      int pk=1;
      while(p%i==0) p/=i,pk*=i;
      ret=(ret+C(n,m,i,pk))%mod;
    }
  }
  return ret;
}

多项式大全

#include<bits/stdc++.h>
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define sz(a) int(a.size())
#define clr(a) memset(a,0,sizeof(a))
#define all(a) a.begin(),a.end()
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,int> pli;
const int inf=1e9;
const ll Inf=1e18;
const int N=4e5+10;
const int mod=998244353;
int gi() {
    int x=0,o=1;char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
    if(ch=='-') o=-1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    return x*o;
}
template<typename T> bool chkmax(T &a,T b) { return a<b?a=b,1:0; };
template<typename T> bool chkmin(T &a,T b) { return a>b?a=b,1:0; };
int add(int a,int b) { return a+b>=mod?a+b-mod:a+b; }
int sub(int a,int b) { return a-b<0?a-b+mod:a-b; }
void inc(int &a,int b) { a=(a+b>=mod?a+b-mod:a+b); }
void dec(int &a,int b) { a=(a-b<0?a-b+mod:a-b); }
int qpow(int a,int b) {
    int ret=1;
    while(b) {
        if(b&1) ret=1ll*ret*a%mod;
        a=1ll*a*a%mod,b>>=1;
    }
    return ret;
}
struct poly {
    vector<int> v;
    int& operator[](int x) { return v[x]; }
    void set(int n) { v.resize(n+1); }
    int get() { return sz(v)-1; }
    void rev() { reverse(all(v)); }
    void push_back(int x) { v.pb(x); }
    void clear() { v.clear(); }
};
void NTT(poly &a,int n,int f,vector<int> &r) {
    for(int i=0;i<n;i++)
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<n;i<<=1) {
        int wn=qpow(3,(mod-1)/(i<<1));
        if(f==-1) wn=qpow(wn,mod-2);
        for(int p=i<<1,j=0;j<n;j+=p) {
            int w=1;
            for(int k=0;k<i;k++,w=1ll*w*wn%mod) {
                int x=a[j+k],y=1ll*w*a[j+k+i]%mod;
                a[j+k]=add(x,y),a[j+k+i]=sub(x,y);
            }
        }
    }
    if(f==-1) {
        int inv=qpow(n,mod-2);
        for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
    }
}
poly operator+(poly a,poly b) {
    int n=max(a.get(),b.get());
    a.set(n),b.set(n);
    for(int i=0;i<=n;i++) inc(a[i],b[i]);
    return a;
}
poly operator-(poly a,poly b) {
    int n=max(a.get(),b.get());
    a.set(n),b.set(n);
    for(int i=0;i<=n;i++) dec(a[i],b[i]);
    return a;
}
poly operator-(int b,poly a) {
    poly c;c.set(0);c[0]=b;
    return c-a;
}
poly operator*(int b,poly a) {
    int n=a.get();
    for(int i=0;i<=n;i++) a[i]=1ll*a[i]*b%mod;
    return a;
}
poly operator*(poly a,poly b) {
    int n=a.get(),m=b.get(),l=0;
    m+=n;for(n=1;n<=m;n<<=1) ++l;
    vector<int> r(n);
    for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    a.set(n-1),b.set(n-1);
    NTT(a,n,1,r),NTT(b,n,1,r);
    for(int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%mod;
    NTT(a,n,-1,r);
    a.set(m);return a;
}
poly getinv(poly a,int n) {
    poly ret;
    if(n==1) { ret.set(0);ret[0]=qpow(a[0],mod-2);return ret; }
    ret=getinv(a,n>>1);
    a.set(n-1);ret=2*ret-a*ret*ret;ret.set(n-1);
    return ret;
}
int get2(int n) {
    int ret=1;
    while(ret<=n) ret<<=1;
    return ret;
}
pair<poly,poly> operator/(poly a,poly b) {
    poly d,r;
    int n=a.get(),m=b.get();
    a.rev(),b.rev();
    poly invb=getinv(b,get2(n-m));
    d=a*invb;d.set(n-m);
    a.rev(),b.rev(),d.rev();
    r=a-b*d;r.set(m-1);
    return mp(d,r);
}
poly operator%(poly a,poly b) {
    return (a/b).se;
}
#define ls (o<<1)
#define rs (o<<1|1)
void solve1(vector<int> &v,poly *A,int o,int l,int r) {
    if(l==r) {
        A[o].set(1);
        A[o][0]=sub(0,v[l]),A[o][1]=1;
        return;
    }
    int mid=(l+r)>>1;
    solve1(v,A,ls,l,mid);
    solve1(v,A,rs,mid+1,r);
    A[o]=A[ls]*A[rs];
}
void solve2(vector<int> &v,poly *A,poly B,int o,int l,int r) {
    if(l==r) { v.pb(B[0]);return; }
    int mid=(l+r)>>1;
    solve2(v,A,B%A[ls],ls,l,mid);
    solve2(v,A,B%A[rs],rs,mid+1,r);
}
vector<int> getval(poly a,vector<int> v) {
    static poly A[N<<2];
    int n=a.get(),m=v.size();n=max(n,m);
    a.set(n),v.resize(n);
    solve1(v,A,1,0,n-1);
    v.clear();
    solve2(v,A,a,1,0,n-1);
    v.resize(m);return v;
}
poly getdao(poly a) {
    int n=a.get();
    for(int i=1;i<=n;i++) a[i-1]=1ll*i*a[i]%mod;
    a.set(n-1);return a;
}
poly getjifen(poly a) {
    int n=a.get();a.set(n+1);
    vector<int> inv(n+2);
    inv[1]=1;for(int i=2;i<=n+1;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    for(int i=n;~i;i--) a[i+1]=1ll*inv[i+1]*a[i]%mod;
    a[0]=0;return a;
}
poly getln(poly a,int n) {
    a=getjifen(getdao(a)*getinv(a,get2(n)));
    a.set(n);return a;
}
poly getexp(poly a,int n) {
    poly ret;
    if(n==1) { ret.set(0);ret[0]=1;return ret; }
    ret=getexp(a,n>>1);
    a.set(n-1);
    ret=ret*(1-getln(ret,n-1)+a);
    ret.set(n-1);return ret;
}
poly getpoly(int o,int l,int r,poly *A,vector<int> &v) {
    if(l==r) {
        poly ret;ret.set(0);
        ret[0]=v[l];return ret;
    }
    int mid=(l+r)>>1;
    poly L=getpoly(ls,l,mid,A,v),R=getpoly(rs,mid+1,r,A,v);
    return L*A[rs]+R*A[ls];
}
poly getpoly(vector<int> &x,vector<int> &y) {
    static poly A[N<<2];
    int n=x.size();
    solve1(x,A,1,0,n-1);
    poly daoA=getdao(A[1]);
    vector<int> v=getval(daoA,x);
    for(int i=0;i<n;i++) v[i]=1ll*y[i]*qpow(v[i],mod-2)%mod;
    return getpoly(1,0,n-1,A,v);
}
int main() {
    
    return 0;
}

转载于:https://www.cnblogs.com/gczdajuruo/p/8455771.html