前言
以下代码多项式系数均对$998244353(7\times 17\times 2^{23}+1)$取膜,如果想对于$1004535809$等其他模数通用的话,请将原根,原根的逆和虚数单位,$\texttt{Cipolla}$中的$\texttt{exponent}$,以及三角函数代码里的$\texttt{inv2}$更改即可。
如果有什么比如说像链接有问题或者是别的问题,请在下面发评论。
多项式的求导与不定积分
$\frac{\mathrm{d}}{\mathrm{d}x}x^n=nx^{n-1},\int x^n=\frac{x^{n+1}}{n+1}\,\mathrm{d}x$,所以直接套公式即可。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16inline void deriv(ll fd,ll *f,ll *res)
{
for(register int i=1;i<fd;++i)
{
res[i-1]=(li)f[i]*i%MOD;
}
res[fd-1]=0;
}
inline void integ(ll fd,ll *f,ll *res)
{
for(register int i=1;i<fd;++i)
{
res[i]=(li)f[i-1]*qpow(i,MOD-2)%MOD;
}
res[0]=0;
}
多项式乘法
略
多项式求逆
如果$F(x)$只有一项,那么答案是它的逆元。
剩下的情况考虑递归做,假设已经求出了$\bmod{x^{\lceil\frac{n}{2}\rceil}}$意义下的答案$H(x)$,考虑推出$\bmod{x^n}$意义下的答案$G(x)$,于是有:
于是可以快活地上板子啦qwq1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30inline void inv(ll fd,ll *f,ll *res)
{
static ll tmp[MAXN];
if(fd==1)
{
res[0]=qpow(f[0],MOD-2);
return;
}
inv((fd+1)>>1,f,res);
ll cnt=1,limit=-1;
while(cnt<(fd<<1))
{
cnt<<=1,limit++;
}
for(register int i=0;i<cnt;++i)
{
tmp[i]=i<fd?f[i]:0;
rev[i]=(rev[i>>1]>>1)|((i&1)<<limit);
}
NTT(tmp,cnt,1),NTT(res,cnt,1);
for(register int i=0;i<cnt;++i)
{
res[i]=(li)(2-(li)tmp[i]*res[i]%MOD+MOD)%MOD*res[i]%MOD;
}
NTT(res,cnt,-1);
for(register int i=fd;i<cnt;++i)
{
res[i]=0;
}
}
多项式除法
首先定义$F_R(x)$是指$F(x)$系数反转后所得的多项式。
很显然,对于$n$次多项式$F(x),F_R(x)=x^nF(\frac{1}{x})$
还有,对于一个$n$次的多项式$F(x)$和一个$m$次的多项式$G(x)$,它们的商$Q(x)$是$n-m$次的,余数$R(x)$是$m-1$次的。
所以,有了这些就可以愉快的推公式了qwq
把$\frac{1}{x}$代入
两边乘个$x^n$
把它换成$\bmod {x^{n-m+1}}$意义下的式子
通过这个式子求得$Q(x)$,再代入原来的式子求出$R(x)$即可。
精彩的上代码环节
1 | inline void div(ll fd,ll gd,ll *f,ll *g,ll *q,ll *r) |
多项式对数函数
大力推一波公式:
根据链式求导法则$\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{\mathrm{d}y}{\mathrm{d}u}\frac{\mathrm{d}u}{\mathrm{d}x}$和基本公式$\frac{\mathrm{d}}{\mathrm{d}x}\ln x=\frac{1}{x}$,有
再积分回来
于是又可以快活地上板子啦qwq1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24inline void ln(ll fd,ll *f,ll *res)
{
static ll pinv[MAXN],der[MAXN];
ll cnt=1,limit=-1;
while(cnt<(fd<<1))
{
cnt<<=1,limit++;
}
inv(fd,f,pinv),deriv(fd,f,der);
for(register int i=0;i<cnt;i++)
{
rev[i]=(rev[i>>1]>>1)|((i&1)<<limit);
}
NTT(pinv,cnt,1),NTT(der,cnt,1);
for(register int i=0;i<cnt;i++)
{
der[i]=(li)der[i]*pinv[i]%MOD;
}
NTT(der,cnt,-1),integ(fd,der,res);
for(register int i=0;i<cnt;i++)
{
der[i]=pinv[i]=0;
}
}
多项式指数函数
还是大力推公式:
先考虑求导,有
指数函数还是没有消去,看来这种方法布星。
两边取个对数,再把$F(x)$移过去
强行解这个方程布星,考虑牛顿迭代。在膜$x^n$意义下的牛顿迭代式子是这样的(这里不证明)
其中,$F_0(x)$是$\bmod x^{\lceil\frac{n}{2}\rceil}$意义下的答案
所以,设$H(x)=\ln x-F(x)$,有
又因为保证了$F(0)=0$,所以$G(x)$常数项为$1$
于是可以愉快的上板子了qwq
1 | inline void exp(ll fd,ll *f,ll *res) |
多项式开根
Luogu P5277(不要脸的宣传一发)
这里考虑两种做法。
第一种做法,两边取个对数
再取指数
于是就完工啦qwq
另一种做法是用另一种方式推公式
设$H^2(x)\equiv F(x)\pmod {x^{\lceil\frac{n}{2}\rceil}}$,有
接着考虑$F_0\neq 0$的情形。
由于第一种做法和下面的分析,可以很明显的知道$G(x)$的每一项都要乘$F_0$的二次剩余。
使用$\texttt{Cipolla}$求解即可。(这种好东西,两段代码都会写)
但是,考虑$\texttt{Cipolla}$的随机性,有些时候可能不会得到想要的结果,可以多试几次,找到那个结果。(我在造$\texttt{P5277}$的数据时出现了好多个常数项是$1$的情况,然后被迫改了十几个测试点的常数项,$\texttt{mmp}$)
由于一些原因,这里先放上第一种做法的代码qwq(很容易被卡)1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59typedef pair<ll,ll> pii;
inline ll rd()
{
return rand()%MOD;
}
inline bool chkx(ll num)
{
return qpow(num,(MOD-1)>>1)==1;
}
inline pii mul(pii x,pii y,ll k)
{
ll x1=x.first,x2=x.second,y1=y.first,y2=y.second;
ll res1=((li)x1*y1%MOD+(li)x2*y2%MOD*k)%MOD;
ll res2=((li)x2*y1%MOD+(li)x1*y2%MOD)%MOD;
return make_pair(res1,res2);
}
inline ll qres(ll num)
{
if(!chkx(num))
{
return 1;
}
ll k=rd();
while(chkx(((li)k*k-num+MOD)%MOD))
{
k=rd();
}
ll kk=((li)k*k-num+MOD)%MOD,exponent=499122177;
pii res=make_pair(1,0),base=make_pair(k,1);
while(exponent)
{
if(exponent&1)
{
res=mul(res,base,kk);
}
base=mul(base,base,kk),exponent>>=1;
}
return min(res.first,MOD-res.first);
}
inline void sqrt(ll fd,ll *f,ll *res)
{
static ll tsqrt[MAXN];
ll k=f[0],cnt=1;
ln(fd,f,tsqrt);
while(cnt<(fd<<1))
{
cnt<<=1;
}
for(register int i=0;i<=cnt;i++)
{
tsqrt[i]=tsqrt[i]&1?(tsqrt[i]+MOD)>>1:tsqrt[i]>>1;
}
exp(fd,tsqrt,res),k=qres(k);
for(register int i=0;i<cnt;i++)
{
res[i]=(li)res[i]*k%MOD;
tsqrt[i]=0;
}
}
多项式快速幂
最好不要用倍增快速幂,尽管大概率卡不掉,但是码量极大。
先推一波公式
两边取对数
再取指数
但是,如果$F(x)$的常数项不为$1$,那么代码中求$\ln$的地方会自动将$F(x)$转换成常数项为$1$的多项式,也就是除掉常数项。
这时,要考虑把结果乘上一个$F_0^k$
可是$F_0$可以等于$0$,这时采用先把$0$删掉再补$0$的方法计算。
所以可以愉快的上代码了qwq
1 | inline void qpow(ll fd,ll *f,ll *res,ll exponent) |
多项式三角函数
像求导积分之类的肯定布星,因为
然后这里试一下求导。
于是,求导很明显解决不了问题。
牛顿迭代也会把你劝退,因为式子也会很复杂,这里就不推导了。
$\texttt{So?}$
考虑使用欧拉公式,将$F(x)$代入
再把$-F(x)$代入
诱导公式搞一下
两式相加相减
把那些系数除过去
等等,这个$\mathrm{i}$是怎么回事,乱搞一下
很明显蒯一个$\texttt{Cipolla}$的板子,算出来答案是$86583718$($911660635$也是可以滴)
于是就可以考虑上代码啦qwq
这里窝把$\sin$和$\cos$放在一起,为了节省空间1
2
3
4
5
6
7
8
9
10
11
12
13
14
15inline void trig(ll fd,ll *f,ll *res,ll type)
{
ll inv2=499122177,inv2i=qpow(I<<1,MOD-2);
static ll tmp[MAXN],tmp2[MAXN],texp[MAXN],texp2[MAXN];
for(register int i=0;i<fd;i++)
{
tmp[i]=(li)I*f[i]%MOD,tmp2[i]=MOD-tmp[i];
}
exp(fd,tmp,texp),exp(fd,tmp2,texp2);
for(register int i=0;i<fd;i++)
{
res[i]=type?(texp[i]+texp2[i])%MOD:(texp[i]-texp2[i]+MOD)%MOD;
res[i]=(li)res[i]*(type?inv2:inv2i)%MOD;
}
}
多项式反三角函数
先讨论$G(x)=\sin^{-1}F(x)$的情形。
两边同时求个导
再积分回来
同理,$G(x)=\tan^{-1}F(x)$就会有
积分回来
所以上个板子qwq(这里$\sin^{-1}$和$\tan^{-1}$分开了)
1 | inline void asin(ll fd,ll *f,ll *res) |
多项式多点求值
把求值的点均分成两份,并构造多项式
显然,对于$\forall i\in [1,\lfloor\frac{n}{2}\rfloor]$,$G_0(x_i)=0$,$G_1$同理。
于是拿$F(x)$(也就是给定的多项式)暴力除$G_0(x)$,有
这样,对于$\forall i\in [1,\lfloor\frac{n}{2}\rfloor]$,就有$F(x_i)=R_0(x_i)$。
同理,构造
对于$\forall i\in [\lfloor\frac{n}{2}\rfloor+1,n]$,就有$F(x_i)=R_1(x_i)$。
所以,问题转化为前$\lfloor\frac{n}{2}\rfloor$个点对$G_0(x)$求值和后$n-\lfloor\frac{n}{2}\rfloor$个点对$G_1(x)$求值的问题,然后就可以愉快的分治啦qwq,分治的过程可以用类似线段树的方法维护(一定要开$4$倍内存)
这里有一个细节:如果$F(x)$的最高次数大于$\prod\limits_{i=0}^N(x-x_i)$要先对$\prod\limits_{i=0}^N(x-x_i)$模一下。
但是,毒瘤的出题人把时间给卡成了$\texttt{1s}$,这样做是没有分的。于是可以在分治到一个地方时循环展开秦九韶暴力算($\texttt{wys}$在召唤)。经人肉二分可以算出来是$1024$。
可是,由于评测姬的不稳定性,交此题要避开评测高峰期,在一个夜深人静的时候交上去,才有可能$\texttt{AC}$啦(尽管我是在评测高峰期交的)
终于可以上代码啦qwq(写的又长又臭,包括了优化后的多项式取模)
1 | inline void mod(ll fd,ll gd,ll *f,ll *g,ll *r) |
多项式快速插值(凑到$800$行啦。。。)
代码全家福
1 | inline ll qpow(ll base,ll exponent) |