多项式基础学习笔记(2)
FWT
我们平时说的多项式卷积（就是 FFT 那个）是加法卷积，也就是 \(\sum\limits_{i+j=k}f_ig_j\) ，而 FWT 是用来解决位运算卷积的，比如与、或和异或。
其实思路是和 FFT 类似的，先求出另一个多项式，然后将对应位置直接乘起来，最后复原。模板题
或卷积
显然这个东西满足交换律，仔细想想发现也满足结合律。那么对于一个最高次项是 \(2^n\) 的多项式 \(A\) ，类似 FFT 的做法分成两部分：前 \(2^{n-1}\) 次项和后 \(2^{n-1}\) 次项，记为 \(A_0\) 和 \(A_1\) 。有
那么 \(IFWT\) 就是
证明的话直接挂链接了 Orz yyb
与卷积
异或卷积
模板
简洁易懂。不过一开始因为没想到 \(\times opt\) 会变成负数然后 WA0 了两发……
//Author: RingweEH
void FWT_Or( int *f,int opt )
{
int i,j,k,len;
for ( i=2; i<lim; i<<=1 )
for ( len=(i>>1),j=0; j<lim; j+=i )
for ( k=j; k<j+len; k++ )
(f[k+len]+=bmod(1ll*f[k]*opt%Mod))%=Mod;
}
void FWT_And( int *f,int opt )
{
int i,j,k,len;
for ( i=2; i<lim; i<<=1 )
for ( len=(i>>1),j=0; j<lim; j+=i )
for ( k=j; k<j+len; k++ )
(f[k]+=bmod(1ll*f[k+len]*opt))%=Mod;
}
void FWT_Xor( int *f,int opt )
{
int i,j,k,len;
for ( i=2; i<lim; i<<=1 )
for ( len=i>>1,j=0; j<lim; j+=i )
for ( k=j; k<j+len; k++ )
{
int t1=f[k],t2=f[k+len];
f[k]=(t1+t2)%Mod; f[k+len]=(t1+Mod-t2)%Mod;
if ( opt==-1 ) f[k]=1ll*f[k]*Inv2%Mod,f[k+len]=1ll*f[k+len]*Inv2%Mod;
}
}
多项式多点求值
给定一个多项式 \(f(x)\) 和 \(m\) 个值 \(a_i\) ，求 \(i\in[1,m],f(a_i)\) . 模板题
令 \(\displaystyle g_0(x)=\prod\limits_{i=1}^{\lfloor \frac{m}{2}\rfloor}(x-a_i),g_1(x)=\prod\limits_{i=\lfloor \frac{m}{2}\rfloor+1}^m(x-a_i)\) 。显然， \(g_0(a_1)=g_0(a_2)=\cdots g_0(a_{\frac{m}{2}})=g_1(a_{\frac{m}{2}+1})=\cdots g_1(a_m)=0\) 。
令 \(f(x)=g_0(x)*q_0(x)+r_0(x)=g_1(x)*q_1(x)+r_1(x)\) 。将 \(a_{1\cdots \frac{m}{2}}\) 带入前一个式子，\(f(x)=r_0(x)\) ，可以继续分治下去。
具体实现上，先分治NTT求出 \(g(x)\) ，求 \(r(x)\) 直接多项式除法即可。
#define lc (pos<<1)
#define rc (pos<<1|1)
int GMod[N<<6],*arr=GMod,*Evg[M];
void Eval_Mul( int n,int m,int *f,int *g,int *res ) //Evaluation
{
static int tf[M],tg[M]; Poly_Init(n);
Copy(tf,f,n); Clear(tf+n,lim-n); NTT( tf,1 );
Copy(tg,g,m); Clear(tg+m,lim-m); NTT( tg,1 );
for ( int i=0; i<lim; i++ ) tf[i]=1ll*tf[i]*tg[i]%Mod;
reverse( tf+1,tf+lim ); NTT( tf,0 ); Copy( res,tf+m-1,n-m+1 );
}
void Eval_Pre( int l,int r,int pos,int *f )
{
Evg[pos]=arr; arr+=r-l+1;
if ( r-l==1 ) { Evg[pos][0]=Mod-f[l],Evg[pos][1]=1; return; }
int mid=(l+r)>>1;
Eval_Pre( l,mid,lc,f ); Eval_Pre( mid,r,rc,f );
Poly_Mul( mid-l+1,r-mid+1,Evg[lc],Evg[rc],Evg[pos] );
}
void Eval_Sol( int l,int r,int pos,int *f,int *g )
{
if ( r-l==1 ) { g[l]=f[0]; return; }
int mid=(l+r)>>1,*fl,*fr;
fl=arr; arr+=mid-l; fr=arr; arr+=r-mid;
Eval_Mul( r-l,r-mid+1,f,Evg[rc],fl );
Eval_Mul( r-l,mid-l+1,f,Evg[lc],fr );
Eval_Sol( l,mid,lc,fl,g ); Eval_Sol( mid,r,rc,fr,g );
}
void Poly_Eval( int n,int m,int *a,int *f,int *g )
{
static int t[M]; n=max( n,m );
Eval_Pre( 0,n,1,a ); reverse( Evg[1],Evg[1]+n+1 );
Poly_Inv( n,Evg[1],t ); reverse( t,t+n );
Poly_Mul( n,n,t,f,t ); Copy( t,t+n,n );
Eval_Sol( 0,n,1,t,g );
for ( int i=0; i<m; i++ ) bmod( g[i]=1ll*g[i]*a[i]%Mod+f[0] );
for ( int i=m; i<n; i++ ) g[i]=0;
}
多项式快速插值
拉格朗日插值复杂度是 \(\mathcal{O}(n^2)\) 的，重心也不能处理 \(k=x[i]\) 的情况，所以需要进行一些奇妙优化。
前置：洛必达法则
如果函数 \(f(x),g(x)\) 满足：
- \(\lim\limits_{x\to x_0}f(x)=\lim\limits_{x\to x_0}g(x)=0\)
- 两个函数在 \(x_0\) 的某去心邻域均可导。（点 \(a\) 的邻域就是以 \(a\) 为中心点的任何开区间，去心就是去掉这个点）
- \(g'(x)\neq 0\)
那么有：
关于邻域和去心邻域可以参考 这篇
回归拉格朗日插值
我们最开始的式子长这样：
把它拆掉
令 \(g(x)=\prod\limits_{i=1}^n(x-x_i)\) ，那么 \(\prod\limits_{j\neq i}(x_i-x_j)=\dfrac{g(x_i)}{x_i-x_i}\) 发现全是0 怎么办呢，当然是上洛必达了！代入之后上式就是 \(g'(x_i)\) 。原式就变成
分治NTT求出 \(g'(x)\) ，然后直接多点求值就能得到 \(\dfrac{y_i}{g'(x_i)}\) 。
现在来求 \(f_{l,r}\) 表示 \((x_{l\sim r},y_{l\sim r})\) 的答案，有
于是就做完了。
代码实现上一个坑点就是，我习惯开 static ，但是递归过程中这个数组是要往下一层递归传参的……于是就暴毙了。还有就是略有卡常，开空间的时候按照所需大小开就可以过了 反正我过了 。
我代码是卡着时限过的……写太丑就不贴了。这里给个 Froggy 的代码 ，我足足跑了 13s ，Froggy 就 9s ……松松怪TAT
二次剩余
定义及基本性质
定义：设 \(m\) 是正整数，\(a\) 是整数，如果 \((a,m)=1\) ，且同余方程 \(x^2\equiv a\pmod m\) 有解，那么称 \(a\) 为模 \(m\) 的 二次剩余 ，否则称为 二次非剩余 。
定理1 ：对于奇素数 \(p\) ，在整数 \(1\sim p-1\) 中，\(p\) 的二次剩余与二次非剩余个数相同。
引理：对于奇素数 \(p\) 和不被 \(p\) 整除的 \(a\) ，同余方程 \(x^2\equiv a\pmod p\) 只可能是无解或者恰有两个模 \(p\) 不同余的解。
引理证明：设 \(x=x_0\) 是其中一个解，那么 \(x=-x_0\) 是不同余的解。而如果有两个解 \(x_0,x_1\) ，那么有 \((x_0+x_1)(x_0-x_1)\equiv 0\) ，所以一定有 \(x_1\equiv -x_0\) 或者相等。
要在 \(1\sim p-1\) 中找出 \(p\) 的所有二次剩余，考虑这 \(p-1\) 个平方，要么无解要么有两个解，所以在 \(1\sim p-1\) 中，\(p\) 的二次剩余恰有 \((p-1)/2\) 个，剩下 \((p-1)/2\) 个是二次非剩余。
定理2 ：设 \(p\) 是素数，\(r\) 是其原根，\(a\) 是不被 \(p\) 整除的整数。如果 \(\text{ind}_r(a)\) （ \(\text{ind}_r(a)\) 表示以 \(r\) 为底 \(a\) 的指数）是偶数，那么 \(a\) 是 \(p\) 的二次剩余，如果是奇数，那么是二次非剩余。
设 \(p\) 是奇素数，整数 \(a\) 不被 \(p\) 整除。定义 勒让德符号 \(\left(\dfrac{a}{p}\right)\) 为：
于是有 欧拉判别法 （或者叫欧拉准则）：设 \(p\) 是奇素数，\(a\) 是不被 \(p\) 整除的正整数，则
对于 \(=1\) 的情况，由费马小定理有 \(\sqrt a^{p-1}\equiv 1\pmod p\) ；对于 \(=0\) 显然。
考虑 \(x^2\equiv a\pmod p\) 无解的情况。对于每个 \((i,p)=1\) 的整数 \(i\) ，存在整数 \(j\) 使得 \(ij\equiv a\pmod p\) ，且显然 \(i\neq j\) . 于是 \(1\sim p-1\) 可以分成 \((p-1)/2\) 对，每一对乘积为 \(a\) ，相乘有 \((p-1)!\equiv a^{(p-1)/2}\pmod p\) ，由威尔逊定理，有 \(a^{(p-1)/2}\equiv -1\) 。
因此有推论： \(\left(\dfrac ap\right)\left(\dfrac bp\right)=\left(\dfrac{ab}{p}\right)\) ，一个重要的特殊情形是 \(\left(\dfrac{a^2}{p}\right)\equiv a^{p-1}\) .
求二次剩余：Cipolla
模板题 求解方程 \(x^2\equiv N\pmod p\) ，保证 \(p\) 是奇素数。
做法如下：在 \([0,p-1]\) 中随机产生一个数 \(a\) ，令 \(w=a^2-n\) ，如果 \(\left(\dfrac wp\right)=-1\) ，那么 \((a+\sqrt w)^{(p+1)/2}\) 是一个二次剩余。
（虽然 \(w\) 是个非二次剩余并不能模意义开方，但是可以定义 \(\sqrt w\) 为一个类似实数 \(\to\) 虚数中产生的 \(i\) 的东西，这样就能计算了。这个好像叫 扩展数域 ）
下面来简单证明一下：
首先，根据二项式定理不难得出模意义下有结论：\((a+b)^p\equiv a^p+b^p\pmod p\) .
令 \(t=(a+\sqrt w)^{(p+1)/2}\) ，那么根据上面关于勒让德的推论和费马小定理， \(t^2\equiv(a^p+\sqrt w^p)(a+\sqrt w)\equiv(a-\sqrt w)(a+\sqrt w)\) 。
再看看 \(w=a^2-n\) ，那么 \(t^2\equiv a^2-w\equiv n\pmod p\) .
证毕。
实现上，由于 \(w\) 是已知的，所以直接把 \(\sqrt w\) 当成虚数中类似 \(i\) 的东西，重载运算符即可。
//Author: RingweEH
//P5491 【模板】二次剩余
ll Mod,w,n;
struct Complex
{
ll x,y;
Complex( ll _x=0,ll _y=0 ) : x(_x),y(_y) {}
Complex operator * ( Complex t ) { return Complex(x*t.x+y*t.y%Mod*w,x*t.y+y*t.x); }
Complex operator % ( ll Mod ) { return Complex( (x%Mod+Mod)%Mod,(y%Mod+Mod)%Mod ); }
};
ll power( ll a,ll b,ll Mod ) { ll res=1; for ( ; b; b>>=1,a=a*a%Mod ) if ( b&1 ) res=res*a%Mod; return res; }
Complex power( Complex a,ll b,ll Mod ) { Complex res(1,0); for ( ; b; b>>=1,a=a*a%Mod ) if ( b&1 ) res=res*a%Mod; return res; }
ll Solve( ll n,ll P )
{
n%=P;
if ( power(n,(P-1)>>1,P)==P-1 ) return -1; //判断无解
ll a=0;
while ( 1 )
{
a=rand()%P; w=(a*a+P-n)%P;
if ( power(w,(P-1)>>1,P)==P-1 ) break;
}
return power( Complex(a,1),(P+1)>>1,P ).x%P;
}
多项式三角函数
复变函数中的欧拉公式将复指数函数与三角函数联系在一起：
具体证明和解析可以看 这篇专栏 。在上式中用 \(x=-x\) 带入，有 \(e^{i(-x)}=\cos x-i\sin x\) ，那么我们就有三角函数的另一个表达式：
我们现在要求多项式 \(f(x)\) 的三角函数，也就是：
这样就可以用之前所学的 \(\exp\) 来完成三角函数求解。注意，由于我们是在模意义下做 NTT ，那么虚数单位 \(i\) 也要换成模意义下的 \(-1\) 的平方根，也就是 \(\sqrt{998244352}\equiv 86583718\pmod{998244353}\) . （虽然这里直接给出了答案，但结合上一节你会知道这就是二次剩余求解的结果而已）
const int TriI=86583718,Inv2=power(2,Mod-2),InvI=power(TriI,Mod-2);
void Poly_SinCos( int n,int *f,int *g,int opt ) //0：sin，1：cos
{
Poly_Init(n); static int tf[M];
Copy(tf,f,n); Clear(tf+n,lim-n);
for ( int i=0; i<n; i++ ) tf[i]=1ll*tf[i]*TriI%Mod;
Poly_Exp( n,tf,g ); Clear(tf,lim); Poly_Inv( n,g,tf ); //g:exp(if(x)),tf:exp(-if(x))
if ( !opt )
{
for ( int i=0; i<n; i++ ) bmod(g[i]+=Mod-tf[i]);
int cons=1ll*Inv2*InvI%Mod;
for ( int i=0; i<n; i++ ) g[i]=1ll*g[i]*cons%Mod;
Clear(g+n,lim-n);
}
else
{
for ( int i=0; i<n; i++ ) bmod(g[i]+=tf[i]);
for ( int i=0; i<n; i++ ) g[i]=1ll*g[i]*Inv2%Mod;
Clear(g+n,lim-n);
}
}
参考资料
Froggy ：多项式大杂烩
《初等数论及其应用》 ：第十一章 二次剩余
C3H5ClO ：二次剩余与高斯整数
OI-Wiki ：多项式三角函数
洛谷题解
