【MOBAN】多项式算法模板和小记录

多项式乘法，多项式求逆，多项式除法，多项式开方，多项式求导和积分，多项式取模，多项式指数函数，多项式对数函数，多项式快速幂，多项式多点求值，多项式多点插值，常系数线性递推，多项式求复合逆,拆系数ＦＦＴ参考文章：https://blog.csdn.net/kscla/article/details/79356786

预备

不用long long 多件套和LJ NTT

const int maxn = 300005,mod=998244353;
const int inv2 = 499122177;
const int g = 3;
int add(int x,int y) {
    x+=y;return x>=mod?x-mod:x;
}
int sub(int x,int y) {
    x-=y;return x<0?x+mod:x;
}
int mul(int x,int y) {
    return 1ll*x*y%mod;
}
int ksm(int a,int b) {
    int ans=1;
    for(;b;b>>=1,a=mul(a,a))
        if(b&1) ans=mul(ans,a);
    return ans;
}

void ntt(int *a,int s,int dft) {
    for(int i=0,j=0;i<s;i++) {
        if(i<j) swap(a[i],a[j]);
        for(int k=(s>>1);(j^=k)<k;k>>=1);
    }
    for(int st=1;st<s;st<<=1) {
        int dwg = dft==1?ksm(g,(mod-1)/(st<<1)):ksm(g,mod-1-(mod-1)/(st<<1));
        for(int i=0;i<s;i+=(st<<1)) {
            int ng = 1;
            for(int j=i;j<i+st;j++) {
                int x=a[j]; int y=mul(a[j+st],ng);
                ng = mul(ng,dwg);
                a[j]=add(x,y); a[j+st]=sub(x,y);
            }
        }
    }
    if(dft==1) return;
    int invs = ksm(s,mod-2);
    for(int i=0;i<s;i++) a[i]=mul(a[i],invs);
}

多项式乘法

void MUL(int *a,int *b,int *c,int le) {
    for(int i=0;i<le;i++) ta[i]=b[i];
    for(int j=0;j<le;j++) tb[j]=c[j];
    ntt(ta,le,1); ntt(tb,le,1);
    for(int i=0;i<le;i++) ta[i]=mul(ta[i],tb[i]);
    ntt(ta,le,-1);
    for(int i=0;i<le;i++) a[i] = ta[i];
}

多项式求逆：

求Ａ(x)*B(x)==1(mod x^n) 设定A(x)*G(x)==1(mod x^(n/2) 那么就有B(x) – G(x) == 0 (mod x^(n/2) ) 就有B(x)^2 + G(x)^2 – 2B(x)G(x) == 0 (mod x^(n)) 两边同乘上A，就最终可以得到 B(x) = G(x)(2-AG(x))了

void ginv(int deg,int *a,int *b) {
    if(deg==1) { b[0]=ksm(a[0],mod-2); return; }
    ginv(deg>>1,a,b);
    for(int i=0;i<deg;i++) ta[i]=a[i];
    for(int i=deg;i<(deg<<1);i++) ta[i]=0;
    ntt(ta,deg<<1,1); ntt(b,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i])));
    ntt(b,(deg<<1),-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}

多项式除法

我们设定Ar（x) = x^n*A(1/x) 比如举个例子A(x) = x^3 + 2x^2 + 4x + 1 Ar(x) = 1+2x+4*x^3+x^3 这恰好等于A(x)的系数反转。考虑原式我们只要将x替换成1/x，会发现一些有趣的事情。 F(1/x)x^n == Q(1/x)G(1/x)x^n + x^nR(x) Fr(x) = Qr(x) + Gr(x) + Rr(x)x^(n-m+1) 那么Ar(x) == Dr(x)Br(x) ( mod x^(n-m+1))

#include<stdio.h>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
const int mod = 998244353;
const int g = 3;
const int maxn = 700005;
int mul(int x,int y) { return 1ll*x*y%mod; }
int add(int x,int y) { x+=y; return x>=mod?x-mod:x;}
int sub(int x,int y) { x-=y; return x<0?x+mod:x;}
int ksm(int a,int b) {
    int ans=1;
    for(;b;b>>=1,a=mul(a,a))
        if(b&1) ans=mul(ans,a);
    return ans;
}
int n,m; int tmp[maxn],ta[maxn],tb[maxn];
void ntt(int *a,int s,int dft) {
    for(int i=0,j=0;i<s;i++) {
        if(i<j) swap(a[i],a[j]);
        for(int k=(s>>1);(j^=k)<k;k>>=1);
    }
    for(int st=1;st<s;st<<=1) {
        int dwg = dft==1?ksm(g,(mod-1)/(st<<1)):ksm(g,mod-1-(mod-1)/(st<<1));
        for(int i=0;i<s;i+=(st<<1)) {
            int ng = 1;
            for(int j=i;j<i+st;j++) {
                int x=a[j]; int y=mul(a[j+st],ng);
                ng = mul(ng,dwg);
                a[j]=add(x,y); a[j+st]=sub(x,y);
            }
        }
    }
    if(dft==1) return;
    int invs = ksm(s,mod-2);
    for(int i=0;i<s;i++) a[i]=mul(a[i],invs);
}
void ginv(int deg,int *a,int *b) {
    if(deg==1) { b[0]=ksm(a[0],mod-2); return; }
    ginv(deg>>1,a,b);
    for(int i=0;i<deg;i++) ta[i]=a[i];
    for(int i=deg;i<(deg<<1);i++) ta[i]=0;
    ntt(ta,deg<<1,1); ntt(b,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i])));
    ntt(b,(deg<<1),-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}
int Q[maxn],F[maxn],G[maxn],Fr[maxn],Gr[maxn],Gri[maxn];
void MUL(int *a,int *b,int *c,int le) {
    for(int i=0;i<le;i++) ta[i]=b[i];
    for(int j=0;j<le;j++) tb[j]=c[j];
    ntt(ta,le,1); ntt(tb,le,1);
    for(int i=0;i<le;i++) ta[i]=mul(ta[i],tb[i]);
    ntt(ta,le,-1);
    for(int i=0;i<le;i++) a[i] = ta[i];
}
int main() {
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++) {
        scanf("%d",&F[i]); Fr[n-i]=F[i];
    }
    for(int i=0;i<=m;i++) {
        scanf("%d",&G[i]); Gr[m-i]=G[i];
    }
    int le = 1;
    for(;le<=(n-m);le<<=1);
    for(int i=n-m+1;i<=m;i++) Gr[i]=0;
    ginv(le,Gr,Gri);
    le = 1; for(;le<=(n<<1);le<<=1);
    memset(ta,0,sizeof ta); memset(tb,0,sizeof tb);
    MUL(Q,Gri,Fr,le);
    reverse(Q,Q+n-m+1);
    for(int i=n-m+1;i<le;i++) Q[i]=0;
    for(int i=0;i<=n-m;i++) printf("%d ",Q[i]);
    puts("");
    le = 1; for(;le<=(n<<1);le<<=1);
    memset(ta,0,sizeof ta); memset(tb,0,sizeof tb);
    MUL(tmp,Q,G,le);
    for(int i=0;i<m;i++) printf("%d ",sub(F[i],tmp[i]));
}

多项式取模（源自多项式除法）

我们已知除数和商和被除数，那么余就等于被除数－商×除数

int Q[maxn],Fr[maxn],Gr[maxn],Gri[maxn],Yy[maxn];
void gmod(int *O,int *F,int *G,int n,int m) {//F%G==O (F n G m)
    for(int i=0;i<=n;i++) Fr[i] = F[n-i];
    for(int i=0;i<=m;i++) Gr[i] = G[m-i];
    for(int i=n-m+2;i<=m;i++) Gr[i] = 0;
    int le = 1; for(;le<=(n-m);le<<=1);
    for(int i=n-m+1;i<=m;i++) Gr[i]=0;
    ginv(le,Gr,Gri);
    le = 1; for(;le<=(n<<1);le<<=1);
    MUL(Q,Gri,Fr,le);
    reverse(Q,Q+n-m+1);
    for(int i=n-m+1;i<le;i++) Q[i]=0;
    MUL(tmp,Q,G,le);
    for(int i=0;i<m;i++) O[i] = sub(F[i],tmp[i]);
}

多项式开方

考虑和多项式求逆用到的同样的倍增的方法如果我已知B(x)*B(x) == A(x) mod(x^m) 求C(x)使得C(x)C(x) == A(x) mod(x^(2m)) 则(B(x)+C(x))*(B(x)-C(x))==0(mod x^(m)) 那么我们讨论B(x)-C(x)==0情况最后又可以得到C(x)==(A(x)+B(x)^2)/(2*B(x)) （mod x^(2m)) code://代入a，b^2 = a

void ginv(int deg,int *a,int *b) {
    if(deg==1) { b[0]=ksm(a[0],mod-2); return; }
    ginv(deg>>1,a,b);
    for(int i=0;i<deg;i++) ta[i]=a[i];
    for(int i=deg;i<(deg<<1);i++) ta[i]=0;
    ntt(ta,deg<<1,1); ntt(b,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i])));
    ntt(b,(deg<<1),-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}
void gsqrt(int deg,int *a,int *b) {
    if(deg==1) { b[0]=1; return; }
    gsqrt(deg>>1,a,b);
    memset(d,0,sizeof d);
    ginv(deg,b,d);
    for(int i=0;i<deg;i++) c[i]=a[i];
    for(int i=deg;i<(deg<<1);i++) c[i]=0;
    ntt(c,deg<<1,1); ntt(b,deg<<1,1); ntt(d,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(add(b[i],mul(c[i],d[i])),inv2);
    ntt(b,deg<<1,-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}

多项式求导和积分

有关微积分的参考：某博客https://www.cnblogs.com/zwfymqz/p/9127894.html 高中数学书选修2-2 http://jiaofu.yousi.com/compontent/pdf/?url=http://jfpdf.yousi.com/160424030303585634.pdf#page=11] 若[latex] F(x) = \sum_{i=0}^{n}a_{i} * x^i [/latex] 则[latex] F'(x) = \sum_{i=1}^{n} i * a_{i} * x^{i-1} [/latex] 则[latex] \int F(x) = \sum_{i=1}^{n} \frac{a_{i}* x^{i+1} }_{i+1} [/latex] 所以我们如果把多项式看做一个函数就可以求导和积分了。

int direv(int *a,int *b,int len) { //求导
    for(int i=1;i<len;i++) b[i-1]=mul(a[i],i); b[len-1]=0;
}
void inter(int *a,int *b,int len) { //积分
    for(int i=1;i<len;i++) b[i]=mul(a[i-1],ksm(i,mod-2)); b[0]=0; //这里也可以预处理逆元
}

多项式对数函数

给出

次多项式

，求一个mod x^n下的多项式B(x)==ln A(x) 我们知道[latex]\ln x[/latex]的导数为[latex]\frac{1}{x}[/latex],我们设[latex]F(x) = \ln x [/latex] [latex]B(x) = F(A(x))[/latex] [latex]B(x)\prime = F(A(x))\prime * A(x)\prime[/latex] [latex]B(x) = \int \frac {A\prime}{A} [/latex] 直接代入搞就可以了

void direv(int *a,int *b,int len) { //求导
    for(int i=1;i<len;i++) b[i-1]=mul(a[i],i); b[len-1]=0;
}
void inter(int *a,int *b,int len) { //积分
    for(int i=1;i<len;i++) b[i]=mul(a[i-1],ksm(i,mod-2)); b[0]=0; //这里也可以预处理逆元
}
int F[maxn],G[maxn],A[maxn],B[maxn];
void getln(int deg,int *a,int *b) {
    direv(a,A,deg); ginv(deg,a,B);
    ntt(A,deg<<1,1); ntt(B,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) A[i]=mul(A[i],B[i]);
    ntt(A,deg<<1,-1); inter(A,b,deg);
        for(int i=0;i<(deg<<1);i++) A[i]=B[i]=0;
}

多项式指数指数

已知函数

，求

G==e^F (mod)

由于牛顿迭代公式 x_i= x_i-1 - f(x_i-1) / f ' (x_i-1) 那么初始x0 = e^a₀ 显然只有a0等于0时有意义。因为lnG - F = 0 套用牛顿迭代公式得到G_i = G_i-1'(1-lnG_i-1'+F) 就可以搞了。

#include<stdio.h>
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>

using namespace std;
const int mod = 998244353;
const int g = 3;
const int maxn = 400005;
int mul(int x,int y) { return 1ll*x*y%mod; }
int add(int x,int y) { x+=y; return x>=mod?x-mod:x;}
int sub(int x,int y) { x-=y; return x<0?x+mod:x; }
int n; int a[maxn],b[maxn],c[maxn];
int ksm(int a,int b){
    int ans=1; 
    for(;b;b>>=1,a=mul(a,a))
        if(b&1) ans=mul(ans,a);
    return ans;
}
void ntt(int *a,int s,int dft) {
    for(int i=0,j=0;i<s;i++) {
        if(i<j) swap(a[i],a[j]);
        for(int k=(s>>1);(j^=k)<k;k>>=1);
    }
    for(int st=1;st<s;st<<=1) {
        int dwg = dft==1?ksm(g,(mod-1)/(st<<1)):ksm(g,mod-1-(mod-1)/(st<<1));
        for(int i=0;i<s;i+=(st<<1)) {
            int ng = 1;
            for(int j=i;j<i+st;j++) {
                int x=a[j]; int y=mul(a[j+st],ng);
                ng = mul(ng,dwg);
                a[j]=add(x,y); a[j+st]=sub(x,y);
            }
        }
    }
    if(dft==1) return;
    int invs = ksm(s,mod-2);
    for(int i=0;i<s;i++) a[i]=mul(a[i],invs);
}
int ta[maxn];
void ginv(int deg,int *a,int *b) {
    if(deg==1) { b[0]=ksm(a[0],mod-2); return; }
    ginv(deg>>1,a,b);
    for(int i=0;i<deg;i++) ta[i]=a[i];
    for(int i=deg;i<(deg<<1);i++) ta[i]=0;
    ntt(ta,deg<<1,1); ntt(b,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i])));
    ntt(b,(deg<<1),-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}
void direv(int *a,int *b,int len) { //求导
    for(int i=1;i<len;i++) b[i-1]=mul(a[i],i); b[len-1]=0;
}
void inter(int *a,int *b,int len) { //积分
    for(int i=1;i<len;i++) b[i]=mul(a[i-1],ksm(i,mod-2)); b[0]=0; //这里也可以预处理逆元
}
int F[maxn],G[maxn],A[maxn],B[maxn];
void getln(int deg,int *a,int *b) {
    direv(a,A,deg); ginv(deg,a,B);
    ntt(A,deg<<1,1); ntt(B,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) A[i]=mul(A[i],B[i]);
    ntt(A,deg<<1,-1); inter(A,b,deg);
    for(int i=0;i<(deg<<1);i++) A[i]=B[i]=0;
}
int f[maxn];
void EXP(int *a,int *b,int deg) {
    if(deg==1) {b[0]=1; return;}
    EXP(a,b,deg>>1); getln(deg,b,f);
    f[0]=sub(a[0]+1,f[0]);
    for(int i=1;i<deg;i++) f[i]=sub(a[i],f[i]);
    ntt(f,deg<<1,1); ntt(b,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i] = mul(b[i],f[i]);
    ntt(b,deg<<1,-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=f[i]=0;
}
int main() {
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%d",&F[i]);
    int len = 1; for(;len<n;len<<=1);
    EXP(F,G,len);
    for(int i=0;i<n;i++) printf("%d ",G[i]);
}

多项式快速幂

啥，直接二分快速幂? 不不不，这是O（nlognlogk）的，他不够优秀我们知道F(x)^k == exp(ln(F(x)^k)) 并且ln(F(x)^k)==k*ln(F(x)) 我们对多项式求ln，然后系数乘以k之后exp还原就好注意细节！我们发现求ln的时候用到的求导和反积分会损失掉常数项！所以在我们回代回exp之前记得单独处理常数项。

    int len = 1; for(;len<=(2*max(m,n));len<<=1);
    getln(len,F,OF);//求出F的k次快速幂
    memset(F,0,sizeof F);
    for(int i=0;i<len;i++) F[i] = mul(k%mod,OF[i]);
    memset(OF,0,sizeof OF);
    OF[0] = ksm(orz,k%(mod-1));//由于求导和积分会损失掉常数项，为保证牛顿迭代正确，必须要常数项带入
    EXP(F,OF,len);

常系数线性递推

求一个满足

阶齐次线性递推数列

a_i

的第

项。即：[latex]a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i} [/latex] 首先矩阵卡速米的做法很显然，把这个递推矩阵写出来，然后直接爆推就可以了。然而k这么大，一次矩阵乘法就GG了。我们如果能够构造一个数列ci使得[latex] A^{n} = \sum_{i=0}^{k} c_{i} * A_{i} [/latex]。但其实我们完全不用真的算出Ａ^n，因为实际上我们想要得到一个递推矩阵，然后最后乘上一个列向量ＳＴ（给定的初始那ｋ个数)。我们把这个向量乘进去[latex] ST * A^{n} = \sum_{i=0}^{k} c_{i} * ST * A_{i} [/latex]。我们只需要他最低的那项。我们仔细观察发现这里的ST*Ai的最低项不就是STi吗？所以我们只需要构造出这样一个数列ci,我们就可以得到，答案就是[latex] ANS = \sum_{i=0}^{k-1} c_{i} * ST_{i} [/latex]。现在我们考虑怎么构造这个数列。如果我们把一个这个递推矩阵的特征多项式手算出来：

引入特征多项式

我们已知对于一个矩阵A,若有非０行（列）向量y，和常数x，使得Ay = xy，那么就称y为A的特征向量，x为Ａ的特征值。我们移一个项，使得|A-Ex|y=0 (E为单位矩阵)，由于ｙ行列式不为０，那么一定|A-Ex| = 0 ，同时我们称f(x) = |A-Ex|为矩阵Ａ的特征多项式。由于某个神奇的定理（Cayley-Hamilton定理），f(A) = 0。（啥？你说为什么Ａ是一个矩阵而不是一个常数？我们其实把这个多项式f的ｘ替换成矩阵变量就可以了，感性理解的话|A-A*E| = 0) 所以手算出来后发现，这个多项式恰好是[latex] f(A) = (-1)^{k} * (A^k - \sum_{i=1}^{k} a_{i}*A^{k-i}) [/latex] 由于f(A)==0 , 那么　直接消掉(-1)^{k}，得到[latex]f(A) = (A^k - \sum_{i=1}^{k} a_{i}*A^{k-i}) [/latex] 也就是说A^n == A^n % (f(A)) , 之后我们就将这个以矩阵为变量的多项式愉快地降低到了次数k-1。也就是说，我们只要构出了这个多项式，也就是我们要的ci数列对应构造出的多项式。多项式取模即可 luogu 线性递推

// luogu-judger-enable-o2
#include<stdio.h>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<cmath>
using namespace std;
const int maxn = 32005*6;
const int mod = 998244353;
const int g = 3;
int add(int x,int y) { x+=y; return x>=mod?x-mod:x; }
int sub(int x,int y) { x-=y; return x<0?x+mod:x; }
int mul(int x,int y) { return 1ll*x*y%mod; }
int ksm(int a,int b) {
    int ans=1;
    for(;b;b>>=1,a=mul(a,a))
        if(b&1) ans=mul(ans,a);
    return ans;
}

void ntt(int *a,int s,int dft) {
    for(int i=0,j=0;i<s;i++) {
        if(i<j) swap(a[i],a[j]);
        for(int k=(s>>1);(j^=k)<k;k>>=1);
    }
    for(int st=1;st<s;st<<=1) {
        int dwg; if(dft==1) dwg=ksm(g,(mod-1)/(st<<1) ); else dwg=ksm(g,mod-1-(mod-1)/(st<<1));
        for(int i=0;i<s;i+=(st<<1)) {
            int ng = 1;
            for(int j=i;j<i+st;j++) {
                int x = a[j]; int y = mul(a[j+st],ng);
                ng = mul(ng,dwg);
                a[j] = add(x,y); a[j+st] = sub(x,y);
            }
        }
    }
    if(dft==1) return;
    int invs = ksm(s,mod-2);
    for(int i=0;i<s;i++) a[i]=mul(a[i],invs);
}
int ta[maxn],tb[maxn];
void ginv(int *a,int *b,int deg) {// in>>a out<<b
    if(deg==1) { b[0] = ksm(a[0],mod-2); return;}
    ginv(a,b,deg>>1);
    for(int i=0;i<deg;i++) ta[i] = a[i];
    for(int i=deg;i<(deg<<1);i++) ta[i]=0;
    ntt(b,deg<<1,1); ntt(ta,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i])));
    ntt(b,deg<<1,-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}
void MUL(int *a,int *b,int *c,int deg){ //in>>a>>b out<<c
    for(int i=0;i<deg;i++) ta[i]=a[i];
    for(int i=0;i<deg;i++) tb[i]=b[i];
    ntt(ta,deg,1); ntt(tb,deg,1);
    for(int i=0;i<deg;i++) c[i]=mul(ta[i],tb[i]);
    ntt(c,deg,-1);
}
int k;
void clr(int *a) {
    for(int i=0;i<=2*k;i++) a[i]=0;
}
int ar[maxn],br[maxn],bri[maxn],Q[maxn],tmp[maxn];
void gmod(int *a,int *b,int *c,int n,int m) {//in>>a%b out<<c a n b m
    clr(ar); clr(br); clr(bri); clr(Q); clr(tmp);
    for(int i=0;i<=n;i++) ar[i]=a[n-i];
    for(int i=0;i<=m;i++) br[i]=b[m-i];
    for(int i=n-m+2;i<=m;i++) br[i]=0;
    int le=1; for(;le<=(n-m);le<<=1);
    ginv(br,bri,le);
    le = 1; for(;le<=2*n;le<<=1);
    MUL(bri,ar,Q,le);
    reverse(Q,Q+n-m+1);
    for(int i=n-m+1;i<le;i++) Q[i]=0;
    MUL(Q,b,tmp,le);
    for(int i=0;i<m;i++) c[i]=sub(a[i],tmp[i]);
    for(int i=m;i<=2*k;i++) c[i]=0;
}
int n;
int xs[maxn],st[maxn],a[maxn];
int sg[maxn],ans[maxn];
int main() {
  //  freopen("orz.in","r",stdin);
    scanf("%d%d",&n,&k);
    for(int i=1;i<=k;i++) scanf("%d",&xs[i]),xs[i]%=mod,xs[i]+=(xs[i]<0)?mod:0;
    for(int i=0;i<k;i++) {
        scanf("%d",&st[i]),st[i]%=mod,st[i]+=(st[i]<0)?mod:0;
    }
    for(int i=1;i<=k;i++) sg[k-i] = mod-xs[i]; sg[k]=1;
    int le=1; for(;le<=2*k;le<<=1);
    a[1]=1; ans[0]=1;
    while(n) {
        if(n&1) {
            MUL(ans,a,ans,le<<1);
            int zz=(k<<1); while(ans[--zz]==0);
            if(zz>=k) gmod(ans,sg,ans,zz,k);
        }
        MUL(a,a,a,le);
        int zz = (k<<1); while(a[--zz]==0);
        if(zz>=k)gmod(a,sg,a,zz,k);
        n>>=1;
    }
    int AS=0;
    for(int i=0;i<k;i++) AS=add(AS,mul(st[i],ans[i]));
    printf("%d\n",AS);
   // cout<<ans[0]<<endl;
}

bzoj 4161 O(k^2logn)

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
using namespace std;
const int maxn = 60006 , mod = 998244353;
int k; long long n;
int a[maxn],h[maxn];
int mul(int x,int y) { return 1ll*x*y%mod; }
int add(int x,int y) { x+=y; return x>=mod?x-mod:x; }
int b[maxn],c[maxn],tmp[maxn];
void qmul(int *x,int *y) {
    for(int i=0;i<k;i++) {
        for(int j=0;j<k;j++) tmp[i+j] = add(tmp[i+j],mul(x[i],y[j]));
    }
    for(int i=k*2-2;i>=k;i--) {
        for(int j=0;j<=k;j++) 
            tmp[i-k+j] = add(tmp[i-k+j],mul(tmp[i],a[k-j]));
        tmp[i] = 0;
    }
    for(int i=0;i<k;i++) x[i] = tmp[i],tmp[i]=0;
}
int main() {
    scanf("%d%lld",&k,&n); n--;
    for(int i=0;i<k;i++) scanf("%d",&h[i]);
    for(int i=1;i<=k;i++) scanf("%d",&a[i]);
    b[0] = 1; c[1] = 1;
    if(n<k) { printf("%d",h[n]); return 0; }
    for(;n;n>>=1,qmul(c,c))
        if(n&1) qmul(b,c);
    int ans = 0;
    for(int i=0;i<k;i++) ans = add(ans,mul(b[i],h[i]));
    printf("%d",ans);
}

多项式多点求值

给定一个

次多项式

f (x)

，现在请你对于

i \in [1, m]

，求出

f (a_{i})

。我们考虑分治处理，对于

i \in [1,m/2] , 我们构造一个g(x) = \prod（x-ai) (1<=i<=m/2),我们发现对于[1,m/2]中所有的值代入这个函数值都为0，那么我们直接对于[1,m/2]的原函数的f都%上g(x)应该是不会影响（即消除掉若干个g(x)不会影响)。这样就可以同时达到消除大于等于x^(m/2)次项的目的，我们发现到最后剩下的那个常数项的值恰好就是我们要求的f(x)（因为除了常数项都被我们消除掉了）, 时间复杂度为O（nlog^2n)

我们写成一个类似线段树的结构把插值的多项式存下来。

#include<stdio.h>
#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
const int maxn = 64005;
const int mod = 998244353;
const int g = 3;
int mul(int x,int y) { return 1ll*x*y%mod; }
int add(int x,int y) { x+=y; return x>=mod?x-mod:x; }
int sub(int x,int y) { x-=y; return x<0?x+mod:x; }
int ksm(int a,int b) {
    int ans = 1;
    for(;b;b>>=1,a=mul(a,a))
        if(b&1) ans = mul(ans,a);
    return ans;
}
void ntt(int *a,int deg,int dft) {
    for(int i=0,j=0;i<deg;i++) {
        if(i<j) swap(a[i],a[j]);
        for(int k=(deg>>1);(j^=k)<k;k>>=1);
    }
    for(int st=1;st<deg;st<<=1) {
        int dwg = (dft==1) ? ksm(g,(mod-1)/(st<<1)) : ksm(g,mod-1-(mod-1)/(st<<1));
        for(int i=0;i<deg;i+=(st<<1)) {
            int ng = 1;
            for(int j=i;j<i+st;j++) {
                int x = a[j]; int y = mul(ng,a[j+st]);
                ng = mul(ng,dwg);
                a[j] = add(x,y); a[j+st] = add(x,mod-y);
            }
        }
    }
    if(dft==1) return;
    int invs = ksm(deg,mod-2);
    for(int i=0;i<deg;i++) a[i] = mul(invs,a[i]);
}
int F[18][maxn*4];
int n,m;
int ta[maxn*4],tb[maxn*4];
void MULL(int *a,int *b,int *c,int le) {
    for(int i=0;i<le;i++) ta[i]=b[i];
    for(int j=0;j<le;j++) tb[j]=c[j];
    ntt(ta,le,1); ntt(tb,le,1);
    for(int i=0;i<le;i++) ta[i]=mul(ta[i],tb[i]);
    ntt(ta,le,-1);
    for(int i=0;i<le;i++) a[i] = ta[i];
}
void ginv(int deg,int *a,int *b) {
    if(deg==1) { b[0]=ksm(a[0],mod-2); return; }
    ginv(deg>>1,a,b);
    for(int i=0;i<deg;i++) ta[i]=a[i];
    for(int i=deg;i<(deg<<1);i++) ta[i]=0;
    ntt(ta,deg<<1,1); ntt(b,deg<<1,1);
    for(int i=0;i<(deg<<1);i++) b[i]=mul(b[i],sub(2,mul(ta[i],b[i])));
    ntt(b,(deg<<1),-1);
    for(int i=deg;i<(deg<<1);i++) b[i]=0;
}
int Q[maxn*4],Fr[maxn*4],Gr[maxn*4],Gri[maxn*4],tmp[maxn*4];
void gmod(int *O,int *F,int *G,int n,int m) {//F%G==O (F n G m)
    for(int i=0;i<=n;i++) Fr[i] = F[n-i];
    for(int i=0;i<=m;i++) Gr[i] = G[m-i];
    for(int i=n-m+2;i<=m;i++) Gr[i] = 0;
    int le = 1; for(;le<=(n-m);le<<=1);
    for(int i=n-m+1;i<=m;i++) Gr[i]=0;
    ginv(le,Gr,Gri);
    le = 1; for(;le<=(n<<1);le<<=1);
    MULL(Q,Gri,Fr,le);
    reverse(Q,Q+n-m+1);
    for(int i=n-m+1;i<le;i++) Q[i]=0;
    MULL(tmp,Q,G,le);
    for(int i=0;i<m;i++) O[i] = sub(F[i],tmp[i]);
    for(int i=0;i<le;i++) tmp[i] = Fr[i] = Gr[i] = Gri[i] = Q[i] = 0;
}
vector<int>ve[maxn*2];
int tot,rt,ls[maxn*2],rs[maxn*2],X[maxn];
int tpa[maxn*4],tpb[maxn*4],tpc[maxn*4];
void maketree(int &p,int l,int r) {
    p = ++tot;
    int len = 1; for(;len<r-l+2;len<<=1);
    ve[p].resize(len); // len ci
    if(l==r) {
        ve[p][0] = (mod-(X[l]%mod+mod)%mod); ve[p][1] = 1; return; 
    }
    int mid = (l+r)>>1;
    maketree(ls[p],l,mid); maketree(rs[p],mid+1,r);
    
    int ss = ve[ls[p]].size();
    for(int i=0;i<ss;i++) tpa[i] = ve[ls[p]][i]; for(int i=ss;i<len;i++) tpa[i]=0;
    
        ss = ve[rs[p]].size();
    for(int i=0;i<ss;i++) tpb[i] = ve[rs[p]][i]; for(int i=ss;i<len;i++) tpb[i]=0;

    MULL(tpc,tpa,tpb,len);
    for(int i=0;i<len;i++) ve[p][i] = tpc[i];
}
int G[maxn*4];
int ANS[maxn];
void DC(int dep,int p,int l,int r,int cs) {
    int m = r-l+1;
    for(int i=0;i<=m;i++) G[i] = ve[p][i];
    if(cs>=m) gmod(F[dep],F[dep-1],G,cs,m);
    else {
        for(int i=0;i<=m-1;i++) F[dep][i] = F[dep-1][i];
    }
    for(int i=0;i<=m;i++) G[i] = 0;
    if(l==r) {
        ANS[l] = F[dep][0];
        return;
    }
    int mid = (l+r)>>1;
    DC(dep+1,ls[p],l,mid,m-1);
    DC(dep+1,rs[p],mid+1,r,m-1);
}
int main() {
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++) {
        int x; scanf("%d",&x);
        F[0][i] = x;
    }
    for(int i=1;i<=m;i++) scanf("%d",&X[i]);
    maketree(rt,1,m);
    for(int i=0;i<=m;i++) G[i] = ve[rt][i];
    DC(1,rt,1,m,n);
    for(int i=1;i<=m;i++) printf("%d\n",ANS[i]);
}

多项式多点插值

给定n个点(x1,y1),(x2,y2),(x3,y3)......(xn,yn)求一个n-1次多项式F满足F(x_i) = y_i 如果我们直接考虑拉格朗日插值就会得到一个n^2的做法 [latex] {F(x)=\sum_{i=1}^n{\frac{\prod_{{j}\neq{i}}{({x}-{x_j})}}{\prod_{{j}\neq{i}}{({x_i}-{x_j})}}{y_i}}} [/latex]

多项式求复合逆

预备:拉格朗日反演

若两个多项式，都没有常数项，有A(B(x))=x,已知A，称B为A的复合逆。显然若A(B(x))=x则B(A(x))，且一定满足一次项系数互为逆元。那么如果已知Ａ要求Ｂ的ｎ次项系数，有 [latex][x^n]B(x)=\frac{1}{n} [x^{n-1}] (\frac{x}{A(x)})^{n}[/latex] 同时拉格朗日反演有扩展有 [latex][x^n]A(B(x))=\frac{1}{n}[x^{n-1}]A'(x)(\frac{x}{B(x)})^n[/latex] 就这样我们就可以在nlogn的时间复杂度内完成求解多项式复合逆的某一项系数了

任意模数ＮＴＴ

由于某些毒瘤出题人出的模数并不棒，想毒瘤一下大家，没法ＮＴＴ，ＦＦＴ又精度上有点出事。

多模数ＮＴＴ

一般对于长度为１e5，数范围在1e9的卷积，可以想到其卷积之后范围应该是大约最大1e23（1e51e91e9)的。因此我们可以考虑三个１ｅ９左右的模数来ＮＴＴ之后利用中国剩余定理合并。把这三个模数背下来吧，因为他们很优秀——原根都是３．　998244353,1004535809,469762049

(1<<26)*7,(1<<23)*119,(1<<21)*479

拆系数ＦＦＴ

多模数ＮＴＴ太慢了，那么，快乐的拆系数ＦＦＴ来啦！我们把每个数拆成aM+b的形式（Ｍ是常数），这样的话，就可以把一个Ａ×Ｂ　转化为　(Ａa + Ab) × (Ba + Bb)。然后分为三个区域AaBa , AbBa + AaBb，Aa*Bb分别对应乘M^2，乘M到M^2，不乘任何数。一般为了方便，选择取Ｍ为(1<<15)=32768.缺点是空间消耗比较爆炸，在一些情况下精度也并不友好。手写复数套装code:

struct cd{
    db x,y;
}fa[maxm],fb[maxm];
cd operator+(cd aa,cd bb) {
    return (cd){aa.x+bb.x,aa.y+bb.y};
}
cd operator-(cd aa,cd bb) {
    return (cd){aa.x-bb.x,aa.y-bb.y};
}
cd operator*(cd aa,cd bb) {
    return (cd){aa.x*bb.x-aa.y*bb.y,aa.x*bb.y+aa.y*bb.x};
}
cd operator/(cd aa,db bb) {
    return (cd){aa.x/bb,aa.y/bb};
}
void fft(cd *a,int n,int dft)
{
    for(int j=0,i=0;i<n;i++) {
        if(i<j) swap(a[i],a[j]);
        for(int k=(n>>1);((j^=k)<k);k>>=1);
    }
    for(int st=1;st<n;st<<=1)
    {
        cd dwfg=(cd){cos(dft*pi/st),sin(dft*pi/st)};
        for(int i=0;i<n;i+=(st<<1))
        {
            cd nfg=(cd){1,0};
            for(int j=i;j<i+st;j++)
            {
                cd x=a[j],y=a[j+st]*nfg;
                nfg=nfg*dwfg;
                a[j]=x+y; a[j+st]=x-y;
            }
        }
    }
    if(dft==-1) for(int i=0;i<n;i++) a[i]=a[i]/n;
}
int fta[maxm],ftb[maxm],ftc[maxm],ftd[maxm];
int oz1[maxm],oz2[maxm],oz3[maxm],oz4[maxm];
void MUL(int *c,int *a,int *b,int deg) {
    for(int i=0;i<deg;i++) fa[i]=(cd){1.0*a[i],0},fb[i]=(cd){1.0*b[i],0};
    fft(fa,deg,1); fft(fb,deg,1);
    for(int i=0;i<deg;i++) fa[i]=fa[i]*fb[i];
    fft(fa,deg,-1);
    for(int i=0;i<deg;i++) c[i]=( (long long)(fa[i].x + 0.5) )%mod;
}
void anymul(int *a,int *b,int *c,int deg) { // c = a*b
    for(int i=0;i<deg;i++) {
        fta[i] = a[i]>>10; ftb[i] = a[i]&((1<<10)-1);
        ftc[i] = b[i]>>10; ftd[i] = b[i]&((1<<10)-1);
    }
    MUL(oz1,ftb,ftd,deg);MUL(oz2,ftb,ftc,deg);
    MUL(oz3,ftd,fta,deg);MUL(oz4,fta,ftc,deg);
    for(int i=0;i<deg;i++) c[i]=add(add(oz1[i],mul(add(oz2[i],oz3[i]),1<<10)),mul(oz4[i],1<<20));
}

posted @ 2018-12-04 15:08 Newuser233 阅读(10) 评论(0) 收藏举报

刷新页面返回顶部

Newuser233

来自高中南开，大学华科的蒟蒻。qq:947534376 cf账号:newusers 个人主页 www.newuser.top 从wordpress刚转到cnblog，文章都还没有分类ovo

【MOBAN】多项式算法模板和小记录

预备

多项式乘法

多项式求逆：

多项式除法

多项式取模（源自多项式除法）

多项式开方

多项式求导和积分

多项式对数函数

多项式指数指数

多项式快速幂

常系数线性递推

引入特征多项式

多项式多点求值

多项式多点插值

多项式求复合逆

预备:拉格朗日反演

任意模数ＮＴＴ

多模数ＮＴＴ

拆系数ＦＦＴ

公告

Newuser233

来自高中南开，大学华科的蒟蒻。qq:947534376 cf账号:newusers 个人主页 www.newuser.top 从wordpress刚转到cnblog，文章都还没有分类ovo

【MOBAN】多项式算法模板和小记录

预备

多项式乘法

多项式求逆：

多项式除法

多项式取模（源自多项式除法）

多项式开方

多项式求导和积分

多项式对数函数

多项式指数指数

多项式快速幂

常系数线性递推

引入 特征多项式

多项式多点求值

多项式多点插值

多项式求复合逆

预备:拉格朗日反演

任意模数ＮＴＴ

多模数ＮＴＴ

拆系数ＦＦＴ

公告

引入特征多项式