FWT
我们平时说的多项式卷积(就是 FFT 那个)是加法卷积,也就是 (sumlimits_{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
与卷积
异或卷积
模板
简洁易懂。不过一开始因为没想到 ( imes 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) ,求 (iin[1,m],f(a_i)) . 模板题
令 (displaystyle g_0(x)=prodlimits_{i=1}^{lfloor frac{m}{2} floor}(x-a_i),g_1(x)=prodlimits_{i=lfloor frac{m}{2} floor+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_{1cdots 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)) 满足:
- (limlimits_{x o x_0}f(x)=limlimits_{x o x_0}g(x)=0)
- 两个函数在 (x_0) 的某去心邻域均可导。(点 (a) 的邻域就是以 (a) 为中心点的任何开区间,去心就是去掉这个点)
- (g'(x) eq 0)
那么有:
关于邻域和去心邻域可以参考 这篇
回归拉格朗日插值
我们最开始的式子长这样:
把它拆掉
令 (g(x)=prodlimits_{i=1}^n(x-x_i)) ,那么 (prodlimits_{j
eq 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_{lsim r},y_{lsim r})) 的答案,有
于是就做完了。
代码实现上一个坑点就是,我习惯开 static
,但是递归过程中这个数组是要往下一层递归传参的……于是就暴毙了。还有就是略有卡常,开空间的时候按照所需大小开就可以过了 反正我过了 。
我代码是卡着时限过的……写太丑就不贴了。这里给个 Froggy 的代码 ,我足足跑了 13s
,Froggy
就 9s
……松松怪TAT
二次剩余
定义及基本性质
定义:设 (m) 是正整数,(a) 是整数,如果 ((a,m)=1) ,且同余方程 (x^2equiv apmod m) 有解,那么称 (a) 为模 (m) 的 二次剩余 ,否则称为 二次非剩余 。
定理1 :对于奇素数 (p) ,在整数 (1sim p-1) 中,(p) 的二次剩余与二次非剩余个数相同。
引理:对于奇素数 (p) 和不被 (p) 整除的 (a) ,同余方程 (x^2equiv apmod p) 只可能是无解或者恰有两个模 (p) 不同余的解。
引理证明:设 (x=x_0) 是其中一个解,那么 (x=-x_0) 是不同余的解。而如果有两个解 (x_0,x_1) ,那么有 ((x_0+x_1)(x_0-x_1)equiv 0) ,所以一定有 (x_1equiv -x_0) 或者相等。
要在 (1sim p-1) 中找出 (p) 的所有二次剩余,考虑这 (p-1) 个平方,要么无解要么有两个解,所以在 (1sim p-1) 中,(p) 的二次剩余恰有 ((p-1)/2) 个,剩下 ((p-1)/2) 个是二次非剩余。
定理2 :设 (p) 是素数,(r) 是其原根,(a) 是不被 (p) 整除的整数。如果 ( ext{ind}_r(a)) ( ( ext{ind}_r(a)) 表示以 (r) 为底 (a) 的指数)是偶数,那么 (a) 是 (p) 的二次剩余,如果是奇数,那么是二次非剩余。
设 (p) 是奇素数,整数 (a) 不被 (p) 整除。定义 勒让德符号 (left(dfrac{a}{p} ight)) 为:
于是有 欧拉判别法 (或者叫欧拉准则):设 (p) 是奇素数,(a) 是不被 (p) 整除的正整数,则
对于 (=1) 的情况,由费马小定理有 (sqrt a^{p-1}equiv 1pmod p) ;对于 (=0) 显然。
考虑 (x^2equiv apmod p) 无解的情况。对于每个 ((i,p)=1) 的整数 (i) ,存在整数 (j) 使得 (ijequiv apmod p) ,且显然 (i eq j) . 于是 (1sim p-1) 可以分成 ((p-1)/2) 对,每一对乘积为 (a) ,相乘有 ((p-1)!equiv a^{(p-1)/2}pmod p) ,由威尔逊定理,有 (a^{(p-1)/2}equiv -1) 。
因此有推论: (left(dfrac ap ight)left(dfrac bp ight)=left(dfrac{ab}{p} ight)) ,一个重要的特殊情形是 (left(dfrac{a^2}{p} ight)equiv a^{p-1}) .
求二次剩余:Cipolla
模板题 求解方程 (x^2equiv Npmod p) ,保证 (p) 是奇素数。
做法如下:在 ([0,p-1]) 中随机产生一个数 (a) ,令 (w=a^2-n) ,如果 (left(dfrac wp ight)=-1) ,那么 ((a+sqrt w)^{(p+1)/2}) 是一个二次剩余。
(虽然 (w) 是个非二次剩余并不能模意义开方,但是可以定义 (sqrt w) 为一个类似实数 ( o) 虚数中产生的 (i) 的东西,这样就能计算了。这个好像叫 扩展数域 )
下面来简单证明一下:
首先,根据二项式定理不难得出模意义下有结论:((a+b)^pequiv a^p+b^ppmod p) .
令 (t=(a+sqrt w)^{(p+1)/2}) ,那么根据上面关于勒让德的推论和费马小定理, (t^2equiv(a^p+sqrt w^p)(a+sqrt w)equiv(a-sqrt w)(a+sqrt w)) 。
再看看 (w=a^2-n) ,那么 (t^2equiv a^2-wequiv npmod 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-isin x) ,那么我们就有三角函数的另一个表达式:
我们现在要求多项式 (f(x)) 的三角函数,也就是:
这样就可以用之前所学的 (exp) 来完成三角函数求解。注意,由于我们是在模意义下做 NTT ,那么虚数单位 (i) 也要换成模意义下的 (-1) 的平方根,也就是 (sqrt{998244352}equiv 86583718pmod{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
:多项式三角函数
洛谷题解