• 多项式基础学习笔记(2)


    书接上回

    FWT

    我们平时说的多项式卷积(就是 FFT 那个)是加法卷积,也就是 (sumlimits_{i+j=k}f_ig_j) ,而 FWT 是用来解决位运算卷积的,比如与、或和异或。

    其实思路是和 FFT 类似的,先求出另一个多项式,然后将对应位置直接乘起来,最后复原。模板题

    或卷积

    [C_k=sum_{i|j=k}A_i imes B_j ]

    显然这个东西满足交换律,仔细想想发现也满足结合律。那么对于一个最高次项是 (2^n) 的多项式 (A) ,类似 FFT 的做法分成两部分:前 (2^{n-1}) 次项和后 (2^{n-1}) 次项,记为 (A_0)(A_1) 。有

    [FWT(A)= egin{cases} (FWT(A_0),FWT(A_0)+FWT(A_1)) & n>0\\ A & n=0 end{cases} ]

    那么 (IFWT) 就是

    [IFWT(A')=(IFWT(A_0'),IFWT(A_1')-IFWT(A_0')) ]

    证明的话直接挂链接了 Orz yyb

    与卷积

    [FWT(A)= egin{cases} (FWT(A_0)+FWT(A_1),FWT(A_1)) & ngt0 \\ A & n=0 end{cases}\\ IFWT(A')=(IFWT(A_0')-IFWT(A_1'),IFWT(A_1')) ]

    异或卷积

    [FWT(A)= egin{cases} (FWT(A_0)+FWT(A_1),FWT(A_0)-FWT(A_1)) & n>0\\ A & n=0 end{cases}\\ IFWT(A')=left(dfrac{IFWT(A_0')+IFWT(A_1')}{2},dfrac{IFWT(A_0')-IFWT(A_1')}{2} ight) ]

    模板

    简洁易懂。不过一开始因为没想到 ( 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)

    那么有:

    [limlimits_{x o x_0}dfrac{f(x)}{g(x)}=limlimits_{x o x_0}dfrac{f'(x)}{g'(x)} ]

    关于邻域和去心邻域可以参考 这篇

    回归拉格朗日插值

    我们最开始的式子长这样:

    [f(x)=sum_{i=1}^ny_iprod_{i eq j}dfrac{x-x_j}{x_i-x_j},给定点值为(x_i,y_i) ]

    把它拆掉

    [f(x)=sumlimits_{i=1}^nleft(dfrac{y_i}{prodlimits_{j eq i}(x_i-x_j)}cdot prodlimits_{j eq i}(x-x_j) ight) ]

    (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)) 。原式就变成

    [f(x)=sumlimits_{i=1}^nleft(dfrac{y_i}{g'(x_i)}cdot prodlimits_{j eq i}(x-x_j) ight) ]

    分治NTT求出 (g'(x)) ,然后直接多点求值就能得到 (dfrac{y_i}{g'(x_i)})

    现在来求 (f_{l,r}) 表示 ((x_{lsim r},y_{lsim r})) 的答案,有

    [egin{aligned} f_{l,r} &=sumlimits_{i=l}^rdfrac{y_i}{g'(x_i)}prodlimits_{j=l,i eq j}^r(x-x_j)\\ &=prodlimits_{j=mid+1}^r(x-x_j)f_{l,mid}+prodlimits_{j=l}^{mid}f_{mid+1,r}\\ &=g_{mid+1,r}f_{l,mid}+g_{l,mid}f_{mid+1,r} end{aligned} ]

    于是就做完了。

    代码实现上一个坑点就是,我习惯开 static ,但是递归过程中这个数组是要往下一层递归传参的……于是就暴毙了。还有就是略有卡常,开空间的时候按照所需大小开就可以过了 反正我过了

    我代码是卡着时限过的……写太丑就不贴了。这里给个 Froggy 的代码 ,我足足跑了 13sFroggy9s ……松松怪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)) 为:

    [left(dfrac ap ight)= egin{cases} 1 & a是p的二次剩余\\ -1 & a是p的二次非剩余\\ 0 & p|n end{cases} ]

    于是有 欧拉判别法 (或者叫欧拉准则):设 (p) 是奇素数,(a) 是不被 (p) 整除的正整数,则

    [left(dfrac ap ight)equiv a^{(p-1)/2}pmod 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;
    }
    

    多项式三角函数

    复变函数中的欧拉公式将复指数函数与三角函数联系在一起:

    [e^{ix}=cos x+isin x ]

    具体证明和解析可以看 这篇专栏 。在上式中用 (x=-x) 带入,有 (e^{i(-x)}=cos x-isin x) ,那么我们就有三角函数的另一个表达式:

    [sin x=dfrac{e^{ix}+e^{-ix}}{2i}\\ cos x=dfrac{e^{ix}-e^{-ix}}{2} ]

    我们现在要求多项式 (f(x)) 的三角函数,也就是:

    [sin f(x)=dfrac{exp(if(x))-exp(-if(x))}{2i}\\ cos f(x)=dfrac{exp(if(x))+exp(-if(x))}{2} ]

    这样就可以用之前所学的 (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多项式三角函数

    洛谷题解

  • 相关阅读:
    【纪中集训】2019.08.01【NOIP提高组】模拟 A 组TJ
    有上下界的网络流 学习笔记
    平衡规划后的华丽暴力——莫队算法
    【纪中集训】2019.07.11【NOIP提高组】模拟 B 组TJ
    hdu 2063 过山车 二分匹配(匈牙利算法)
    hdu 5373 The shortest problem(杭电多校赛第七场)
    hdu 3729 I'm Telling the Truth(二分匹配_ 匈牙利算法)
    hdu 2119 Matrix(二分匹配)
    hdu 1498 50 years, 50 colors(二分匹配_匈牙利算法)
    HNU Joke with permutation (深搜dfs)
  • 原文地址:https://www.cnblogs.com/UntitledCpp/p/Poly_AllInOne_2.html
Copyright © 2020-2023  润新知