[学习笔记] 多项式全家桶系列

最新推荐文章于 2021-01-03 11:20:58 发布
最新推荐文章于 2021-01-03 11:20:58 发布
阅读量332 收藏 2
点赞数
原文链接：http://www.cnblogs.com/YoungNeal/p/10209424.html
版权

题外话：
本来这段时间都在洛谷写博客的...
但是因为洛谷猝不及防的改版了导致洛谷博客的搜索功能一度变的很垃圾..
所以又回来了...
等什么时候变好了再回去吧...
哦对了我去洛谷写博客的唯一理由是自己高兴。嗯，真的。

多项式全家桶系列

本来觉得这些没啥用，但是想学新东西了...那就随便学学叭

预处理单位根

upd：学了一下预处理单位根（实际上就是又学了一遍单位根到底是什么）

贴上\(\text{NTT}\)预处理单位根的部分代码：

首先是预处理部分：

w[0][0]=w[1][0]=1;
w[0][1]=ksm(3,(mod-1)/maxn);
w[1][1]=ksm((mod+1)/3,(mod-1)/maxn);
for(int i=2;i<maxn;i++) 
    w[0][i]=1ll*w[0][i-1]*w[0][1]%mod,w[1][i]=1ll*w[1][i-1]*w[1][1]%mod;

这个很显然

然后就是\(\text{NTT}\)里面了，直接把整个函数弄上来吧

void ntt(int *f,int g){
    for(int i=1;i<lim;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        for(int R=mid<<1,j=0;j<lim;j+=R){
            for(int k=0;k<mid;k++){
                int x=f[j+k],y=1ll*w[g][maxn/R*k]*f[j+k+mid]%mod;
                f[j+k]=x+y>=mod?x+y-mod:x+y,f[j+k+mid]=x-y<0?mod+x-y:x-y;
            }
        }
    } if(g)
        for(int in=ksm(lim),i=0;i<lim;i++) f[i]=1ll*f[i]*in%mod;
}

一定一定要记得 \(lim\) 先除再乘！！！！

这个是\(\text{FFT}\)的：

for(int i=0;i<2;i++)
    for(int j=0;j<maxn;j++)
        w[i][j]=cp(cos(2*pi/maxn*j),(i?1:-1)*sin(2*pi/maxn*j));
void fft(cp *f,int opt){
    for(int i=1;i<lim;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        for(int R=mid<<1,j=0;j<lim;j+=R){
            for(int k=0;k<mid;k++){
                cp x=f[j+k],y=f[j+k+mid]*w[opt][maxn/R*k];
                f[j+k]=x+y;f[j+k+mid]=x-y;
            }
        }
    }
}

多项式求逆

Description

给定多项式 \(F(x)\)，请求出一个多项式 \(G(x)\)，满足 \(F(x)G(x)\equiv 1\; \operatorname{mod} x^n\)。对 \(998244353\) 取模。

Sol

推导过程挺简单也挺显然的...

结论就是 \[B(x)\equiv 2B'(x)-A(x)B'^2(x)\; \operatorname{mod} x^n\]

其中 \[A(x)B'(x)\equiv 1\; \operatorname{mod} x^{\frac n2}\]

一个多项式有没有逆取决于常数项有没有逆

Code

void solveinv(int *a,int *b,int len){
    if(len==1) return b[0]=ksm(a[0]),void();
    solveinv(a,b,len>>1); lim=len+len;
    for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
    for(int i=0;i<len;i++) tmpa[i]=a[i];
    for(int i=len;i<lim;i++) tmpa[i]=0;
    ntt(tmpa,3),ntt(b,3);
    for(int i=0;i<lim;i++) b[i]=1ll*b[i]*(2ll-1ll*tmpa[i]*b[i]%mod+mod)%mod;
    ntt(b,g); for(int i=len;i<lim;i++) b[i]=0;
}

多项式除法

Description

给出 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次多项式 \(G(x)\)，请求出多项式 \(Q(x),R(x)\)，满足:

  • \(deg(Q)=n-m,deg(R)<m\)
  • \(F(x)=G(x)Q(x)+R(x)\)

对 \(998244353\) 取模。

Sol

首先，对于任意多项式

\[P(x)=\sum_{i=0}^na_ix^i\]

显然有

\[P^{rev}(x)=x^nP(\frac 1x)=\sum_{i=0}^na_{n-i}x^i\]

也就是说通过这个可以把多项式的系数翻转。

对 \(F(x)=G(x)Q(x)+R(x)\) 进行这个操作：

\[F(x)=G(x)Q(x)+R(x)\]

\[F(\frac 1x)=G(\frac 1x)Q(\frac 1x)+R(\frac 1x)\]

\[x^nF(\frac 1x)=x^nG(\frac 1x)Q(\frac 1x)+x^nR(\frac 1x)\]

\[x^nF(\frac 1x)=x^mG(\frac 1x)x^{n-m}Q(\frac 1x)+x^{n-m+1}x^{m-1}R(\frac 1x)\]

\[F^{rev}(x)=G^{rev}(x)Q^{rev}(x)+x^{n-m+1}R^{rev}(x)\]

注意到 \(x^{n-m+1}R^{rev}(x)\) 的最低项次数至少为 \(n-m+1\)，所以我们对式子取模：

\[F^{rev}(x)=G^{rev}(x)Q^{rev}(x)\; \operatorname{mod}x^{n-m+1}\]

因为 \(deg(Q)=n-m\)，所以这样是不会出问题的。

所以就可以算出 \(Q^{rev}(x)=F^{rev}(x)(G^{rev}(x))^{-1}\)

多项式求逆即可。

整个做法的关键就是消去余数。

Code

void solve(){
    n=getint(),m=getint();
    for(int i=0;i<=n;i++) f[i]=getint();
    for(int i=0;i<=m;i++) g[i]=getint();
    std::reverse(f,f+n+1),std::reverse(g,g+m+1); //翻转系数
    for(int i=0;i<=n-m;i++) ing[i]=g[i];  
    //因为对 x^{n-m+1} 取模，所以后边的系数都应该为0
    getni(ing,q,n-m+1);  //求出g的逆q
    for(int i=0;i<=n-m;i++) c[i]=f[i]; //后边的系数都应该为0
    mul(q,c,n-m+1);
    for(int i=0;i<lim;i++) c[i]=0;
    std::reverse(q,q+n-m+1);
    for(int i=0;i<=n-m;i++) printf("%d ",q[i]);puts("");
    for(int i=n-m+1;i<lim;i++) q[i]=0;
    std::reverse(f,f+n+1),std::reverse(g,g+m+1);
    mul(q,g,n);
    for(int i=0;i<m;i++) printf("%d ",(f[i]-q[i]+mod)%mod);
}

多项式开根

Description

给定多项式 \(A(x)\)，若存在一个多项式 \(B(x)\),满足：

  • \(deg(B)\leq deg(A)\)
  • \(B^2(x)\equiv A(x)\;\operatorname{mod} x^n\)

则称 \(B(x)\) 为 \(A(x)\) 在模 \(x^n\) 意义下的开根。

Sol

懒得写了...

感谢AlphaINF

F49oJH.png

补充一下，上边那个式子还可以再化简一步变成：

\[B(x)\equiv \frac{A(x)}{2C(x)}+\frac{C(x)}2\;\bmod x^n\]

Code

void solvesqr(int len,int *a,int *b){
    if(len==1) return b[0]=1,void();
    solvesqr(len>>1,a,b);
    solveinv(len,b,tmpb);
    get(len);
    for(int i=0;i<len;i++) tmpa[i]=a[i];
    ntt(tmpb,1),ntt(tmpa,1);
    for(int i=0;i<lim;i++) tmpa[i]=1ll*tmpa[i]*tmpb[i]%mod;
    ntt(tmpa,-1);
    for(int i=0,inv2=mod+1>>1;i<lim;i++) b[i]=1ll*(tmpa[i]+b[i])%mod*inv2%mod;
    for(int i=len;i<lim;i++) b[i]=0;
    for(int i=0;i<lim;i++) tmpa[i]=tmpb[i]=0;
}

三模数NTT

Description

求 \(F(x)G(x)\)，对 \(p\) 取模。\(p\) 不保证有原根。

Sol

拆成对三个模数 \(469762049,998244353,1004535809\) 分别做 \(NTT\)，它们的原根都是 \(3\)。

先合并前两个式子，得到:

\[ans\equiv C\;\operatorname{mod}M\]

\[ans\equiv c_3\operatorname{mod} m_3\]

设:

\[ans=xM+C=ym_3+c_3\]

求出 \(x\) 在 \(\operatorname{mod}m_3\) 意义下的值：

\[xM\equiv c_3-C\;\operatorname{mod}m_3\]

在 \(\operatorname{mod} m_3\) 意义下， \(ym_3\) 被消掉了。所以有：

\[x\equiv (c_3-C)M^{-1}\;\operatorname{mod}m_3\]

算出右边的值为 \(q\),令 \(x=km_3+q\)，代入\(ans=xM+C\):

\[ans=(km_3+q)M+C=km_1m_2m_3+qM+C\]

因为 \(ans\in[0,m_1m_2m_3)\)，所以 \(k=0\)！ 所以 \(ans=qM+C\)！

Code

主要是合并部分

int get(int x,int y,int z,int mod){
    ll M=p[1]*p[0];
    ll tmp1=mul(1ll*x*p[1]%M,inv(p[1]%p[0],p[0]),M); //快速幂之前先取模
    ll tmp2=mul(1ll*y*p[0]%M,inv(p[0]%p[1],p[1]),M);
    ll tmp=(tmp1+tmp2)%M;
    ll xx=1ll*(z-tmp%p[2]+p[2])%p[2]*inv(M%p[2],p[2])%p[2];
    return ((xx%mod)*(M%mod)%mod+tmp%mod)%mod; //注意这里先对mod取模再乘起来
}

多项式exp，ln等学完高数马上填坑

拆系数FFT

Step1 基本介绍

设 \(F_P\) 表示多项式 \(P(x)\) DFT之后的结果。

对两个长度为 \(L\) 的多项式 \(A(x),B(x)\) 进行DFT，可以利用
\[ \begin{align}& P(x)=A(x)+iB(x)\\ &Q(x)=A(x)-iB(x) \end{align} \]
对 \(P(x)\) 进行DFT，得到 \(F_p\)。

那么 \(F_q[k]=!(F_p[2L-k])\) ，（!表示取共轭）（证明见论文）

则
\[ \begin{align}DFT(A[k])&=\frac{F_p[k]+F_q[k]}2 \\DFT(B[k])&=-i\frac{F_p[k]-F_q[k]}2 \end{align} \]
这就是两两合并计算DFT的方法，2次DFT优化成了1次。

主要部分的代码：

for(int i=0;i<lim;i++){
    int to=i?lim-i:0;// emm这里并不知道为什么要这样写
    // 嗯现在知道了 这个本质上就是 (lim-i)%lim <--- 这个是对的
    b[i]=(a[i]*a[i]-!a[to]*!a[to])*cp(0,-0.25);
    // 像下面这样写也可以
    // tmpx[i]=(a[i]+!a[to])*cp(0.5,0);
    // tmpy[i]=(a[i]-!a[to])*cp(0,-0.5);
    // b[i]=tmpx[i]*tmpy[i];
}

这样搞一下我的FFT跑的就比NTT还快了！

Step2 拆系数FFT

设 \(M_0=\sqrt M\)，把多项式 \(A(x)\) 和 \(B(x)\) 拆成 \(k\times M_0+b\) 的形式，其中 \(k,b\) 都小于 \(M_0\)

所以两个多项式相乘就变成了 \(K_a(x)K_b(x)M_0^2+(K_a(x)B_b(x)+K_b(x)B_a(x))M_0+B_a(x)B_b(x)\)

然后可以把DFT从四次优化到两次，IDFT从三次优化到两次

这样只需要做四次DFT就行了！

说一下IDFT的两两合并：

假如说现在知道了多项式 \(A(x),B(x),C(x),D(x)\) 的DFT，要计算的是 \(A(x)B(x)\) 和 \(C(x)D(x)\)。令 \(DFT(A)\) 表示多项式 \(A(x)\) 的点值表达。

那么设：
\[ \begin{align}dfta[k]=DFT(A[k])\cdot DFT(B[k])\\dftb[k]=DFT(C[k])\cdot DFT(D[k])\end{align} \]
我们需要对 \(dfta(x)\) 和 \(dftb(x)\) 进行IDFT。注意到这里的IDFT一定是实数，然后经过一番复杂的推导，再令：
\[ p[k]=dfta[k]+i\cdot dftb[k] \]
那么 \(IDFT(p)\) 的实部除以 \(lim\) 就是 \(A(x)B(x)\)，虚部除以 \(lim\) 就是 \(C(x)D(x)\)。

代码：

void mul(int *a,int *b,int *ans,int n){
    static cp p[N],q[N],c[N],d[N];
    for(int i=0;i<lim;i++)
        p[i]=cp(a[i]>>15,a[i]&32767),q[i]=cp(b[i]>>15,b[i]&32767);
    fft(p,1),fft(q,1);
    for(int i=0;i<lim;i++){
        static cp aaa=cp(0.5,0),bbb=cp(0,-0.5),o=cp(0,1);
        static cp ka,kb,ba,bb;
        int j=(lim-i)%lim;
        ka=(p[i]+!p[j])*aaa,kb=(p[i]-!p[j])*bbb;
        ba=(q[i]+!q[j])*aaa;bb=(q[i]-!q[j])*bbb;
        c[i]=ka*ba+kb*ba*o; d[i]=ka*bb+kb*bb*o;
    } fft(c,-1),fft(d,-1);
    for(int i=0;i<=n;i++){
        ll A=((ll)(c[i].x+0.5))%mod,B=((ll)(c[i].y+0.5))%mod,C=((ll)(d[i].x+0.5))%mod,D=((ll)(d[i].y+0.5))%mod;
        ans[i]=((A<<30)%mod+(B+C<<15)%mod+D)%mod;
    }
}

多项式对数函数

Description

给定 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod x^n\) 下的多项式 \(B(x)\)，满足 \(B(x)\equiv \ln A(x)\;\bmod x^n\)。对 \(998244353\) 取模。

Sol

求ln还是挺傻逼的吧...

根据微积分基本知识，复合函数 \(y=f(g(x))\) 的导数和函数 \(y=f(u),u=g(x)\) 的导数间的关系为:

\[y_x'=y_u'\cdot u_x'\]

然后这题就是求 \(y=\ln(A(x))\)

根据上边的定义转变成 \(y'=\frac{A'(x)}{A(x)}\)

求导求逆再积分就行了

然后发现自己求逆的一点小bug就是设置长度设置为不小于 \(n\) 的 \(2\) 的整次幂就行了，不用大于 \(2n\)。

void ds(int *a,int *b,int n){
    for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod;b[n-1]=0;
}
void jf(int *a,int *b,int n){
    for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*inv(i)%mod;b[0]=0;
}
void solveln(int *a,int *b,int len){
    lim=1;while(lim<len) lim<<=1;
    solveinv(a,fa,lim);
    lim=1;while(lim<len+len) lim<<=1;
    for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
    ds(a,da,len);
    ntt(fa,3),ntt(da,3);
    for(int i=0;i<lim;i++) fa[i]=1ll*fa[i]*da[i]%mod;
    ntt(fa,g); jf(fa,b,len);
    for(int i=0;i<lim;i++) fa[i]=da[i]=0;
}

多项式指数函数 (多项式exp)

Description

给定 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod x^n\) 下的多项式 \(B(x)\)，满足 \(B(x)\equiv e^{A(x)}\)。对 \(998244353\) 取模。

Sol

首先要知道牛顿迭代。 在针对一个函数的时候，牛顿迭代可以求函数零点，而针对一个多项式的时候求的就是多项式的零点。即对于一个多项式 \(G(x)\)，求满足条件的 \(G(F(x))\equiv 0\;(\bmod x^n)\) 的多项式 \(F(x)\)。

然后直接放公式就是 \(F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\;(\bmod x^n)\)

其中 \(F_0(x)\) 满足 \(G(F_0(x))\equiv 0\;(\bmod x^{\frac n2})\)

回到多项式exp，那我们现在就是要计算 \(F(x)=e^{A(x)}\)

变形一下得
\[ \ln F(x)-A(x)=0 \]
我们设 \(G(F(x))=\ln F(x)-A(x)\)，那么就是要求这个多项式的零点。如果我们把 \(F(x)\) 看做自变量，那 \(A(x)\) 就是常数，对 \(G(F(x))\) 进行求导，则 \(G'(F(x))=\frac1{F(x)}\)。

代入牛顿迭代的公式得：
\[ F(x)\equiv F_0(x)-F_0(x)\cdot (\ln F(x)+A(x)) \]

\[ F(x)\equiv F_0(x)\cdot (1-\ln F(x)+A(x)) \]

Code

    void solveexp(int *a,int *b,int len){
    if(len==1) return b[0]=1,void();
    solveexp(a,b,len>>1); solveln(b,tmpb,len); lim=len+len;
    for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
    tmpb[0]=(a[0]+1-tmpb[0]+mod)%mod;
    for(int i=1;i<len;i++) tmpb[i]=(a[i]-tmpb[i]+mod)%mod;
    ntt(b,3),ntt(tmpb,3);
    for(int i=0;i<lim;i++) b[i]=1ll*b[i]*tmpb[i]%mod;
    ntt(b,g);
    for(int i=len;i<lim;i++) b[i]=tmpb[i]=0;
    for(int i=0;i<len;i++) tmpb[i]=0;
}

转载于:https://www.cnblogs.com/YoungNeal/p/10209424.html

weixin_30355437
关注
多项式全家桶
Recalling_Clouds的博客
01-23 752
下面的所有知识均属多项式范畴。 前言 本文略去多项式的定义，毕竟初中数学都讲过。 本文会出现一些可能没有见过的表达： 对于一个关于 xxx 的多项式 f(x)f(x)f(x)，有时也称其最高次项次数为 f(x)f(x)f(x) 的度数。 多项式同样定义了四则运算，包括常数与多项式的四则运算和多项式之间的四则运算，同时也衍生出了同余符号 ≡\equiv≡、对数函数 ln⁡(f(x))\ln(f(x))ln(f(x))，指数函数 exp⁡(f(x))\exp(f(x))exp(f(x)) 等等。 下文一般采
[学习笔记] 多项式全家桶 学习笔记
STJqwq的博客
08-20 2649
写在前面：建议初一就学懂这个东西，这是多项式题目的基础。 0.复数类的实现 定义一个复数为 a+bia+bia+bi 。 加法： a+bi+c+di=(a+c)+(b+d)ia+bi+c+di=(a+c)+(b+d)ia+bi+c+di=(a+c)+(b+d)i 减法：同上，系数改成 −1-1−1 乘法：(a+bi)(c+di)=ac+adi+bci+bdi2=(ac−bd)+(ad+bc)i(a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i(a+bi)(c+di)
【模板】多项式全家桶
weixin_33912638的博客
03-21 149
FFT struct Z { double x,y; Z(double _x=0,double _y=0):x(_x),y(_y){}; }; Z operator +(Z x,Z y) {return Z(x.x+y.x,x.y+y.y);} Z operator -(Z x,Z y) {return Z(x.x-y.x,x.y-y.y);} Z operator *(Z x,Z...
【学习笔记】超简单的多项式开方
繁凡さん的博客
01-01 1260
整理的算法模板合集： ACM模板 目录P5205 【模板】多项式开根 点我看多项式全家桶(●^◡_◡◡​^●) P5205 【模板】多项式开根 #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> #include <cmath> using namespace std; typedef long long ll; const int N
【学习笔记】超简单的多项式对数（多项式求导 + 积分）
繁凡さん的博客
01-01 3070
整理的算法模板合集： ACM模板 目录多项式对数P4725 【模板】多项式对数函数（多项式 ln） 点我看多项式全家桶(●^◡_◡◡​^●) 多项式对数 其实非常简单，就是直接求导再积分即可。 P4725 【模板】多项式对数函数（多项式 ln） 直接简单实现一下就行了。 我们输入多项式 fff ，先求fff的导数 AAA ，fff 的逆元 BBB，然后把他们乘起来（NTTNTTNTT），最后求积分 ggg 既是答案，输出即可。 至于求导和积分，因为我们的多项式f(x)=a0+a1×x+a2×x2…an
【学习笔记】超简单的多项式快速幂
繁凡さん的博客
01-03 1629
整理的算法模板合集： ACM模板 目录P5245 【模板】多项式快速幂 点我看多项式全家桶(●^◡_◡◡​^●) 我们知道指数有一个性质： lnxn=n lnxlnx^n=n\ lnxlnxn=n lnx 我们要求的是f(x)kf(x)^kf(x)k 可以很自然地构造出： f(x)k=elnf(x)k=ek lnf(x)f(x)^k=e^{lnf(x)^k}=e^{k\ lnf(x)}f(x)k=elnf(x)k=ek lnf(x) 所以我们只需要求f(x)f(x)
【学习笔记】超简单的多项式求逆（含全套证明）
繁凡さん的博客
12-30 3327
整理的算法模板合集： ACM模板 目录多项式求逆一、分治FFT二、倍增法及其证明三、多项式求逆例题[P4238 【模板】多项式乘法逆](https://www.luogu.com.cn/problem/P4238) 点我看多项式全家桶(●^◡_◡◡​^●) 多项式求逆 一、分治FFT 二、倍增法及其证明 三、多项式求逆例题 P4238 【模板】多项式乘法逆 直接按照上述的思路实现一下模板即可。注意要中间 h(x)h(x)h(x) 要取的是模 x⌈n⌉2x^{\frac{\lceil n \rceil}
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值