zoukankan      html  css  js  c++  java
  • Noip前的大抱佛脚----数论

    数论

    Tags:Noip前的大抱佛脚

    知识点

    Exgcd

    (O(logn))求解(Ax+By=C)的问题
    1、若(C\%gcd(A,B)!=0)则无解
    2、(Gcd=gcd(A,B);A/=Gcd,B/=Gcd,C/=Gcd)
    3、代入下面代码求(Ax+By=1)
    4、(x*C),得到一组特解
    5、通解为(egin{cases}x=x_0+k*B \y=y_0+k*Aend{cases})

    void Exgcd(ll a,ll b,ll &x,ll &y)
    {
    	if(!b){x=1;y=0;return;}
    	Exgcd(b,a%b,y,x);y-=a/b*x;
    }
    

    逆元

    (gcd(A,P)==1)时,(A)在模(P)意义下存在逆元(证明可以符合Exgcd有解的证明),其余情况不存在逆元

    通常(P)为质数时就是(A^{p-2})作为逆元(费马小定理)

    • 欧拉定理(费马小定理) (x^{-1}equiv x^{p-2}(mod p)),要求(p)为质数

    • 解方程(Exgcd)(Px+Ay=1,A=frac{1}{y}(mod P))

    • 线性递推逆元:

      (a=lfloorfrac{p}{i} floor,b=p\%i),则有(ai+b=p,iequiv-frac{b}{a}(mod p),frac{1}{i}equiv-frac{a}{b}(mod p))

      由于(b<i),所以(b)的逆元已经求出,直接可以得到:

      inv[i]=p-(p/i*inv[p%i])%p;
      

    gcd

    具有一些奇妙的性质,如可合并

    • 往数集中加入一个数,要么gcd不变,要么至少变为原来的(frac{1}{2})(所以可以用来分块了(同样的xor也每次只会增加log次))

    欧拉函数(varphi(x))

    表示小于x且与x互质的数的个数

    计算公式

    [varphi(n)=n*prod(1-frac{1}{p_i})$$,其中$p_i$表示$n$的不相同的质因子 #### **欧拉公式** $$a^{varphi(p)}equiv 1 (mod p)]

    降幂公式

    • (x,p)互质:(x^kequiv x^{k\%varphi(p)} (mod p))
    • (x,p)不互质:(x^kequiv egin{cases}x^{k} (mod p),klevarphi(p) \x^{k\%varphi(p)+varphi(p)} (mod p),k>varphi(p)end{cases})

    CRT&EXCRT

    (O(nlogn))求解一系列同余方程,如

    [egin{cases}xequiv A_1(mod P_1) \xequiv A_2(mod P_2) \...\xequiv A_n(mod P_n) \end{cases} ]

    代码如下(未判无解,无解即exgcd无解)

    ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
    void Exgcd(ll a,ll b,ll &x,ll &y)
    {
        if(!b) {x=1;y=0;return;}
        Exgcd(b,a%b,y,x);y-=a/b*x;
    }
    int EXCRT()
    {
        for(int i=2;i<=n;i++)
        {
            ll g=gcd(P[i-1],P[i]),C=(A[i]-A[i-1])/g,x,y;
            Exgcd(P[i-1]/g,P[i]/g,x,y);
            P[i]=P[i-1]*P[i]/g;
            A[i]=A[i-1]+P[i-1]*x*C;
            A[i]=(A[i]%P[i]+P[i])%P[i];
        }
        return A[n];
    }
    

    方法是每次合并两个方程,手推式子(P[i]x+A[i]=P[i+1]y+A[i+1] -> P[i]x+P[i+1]y=A[i+1]-A[i])

    可以发现模数为原来的(lcm)(先除以(gcd)再乘,很容易爆(long long)),余数为原来余数加上模数的(x)

    注意经常x算出来是负数所以时刻(+mod)\%mod !)

    有的时候题目并没有那么裸

    • (x)前带系数(屠龙勇士)

    (Axequiv B(mod p) -> Ax+pk=B(当且仅当B\%gcd(A,p)==0时有解) -> x=x_0+tp -> xequiv x_0(mod p))

    • CRT合并答案

    对于一些公式,只适用于模数是质数的情况,而题目要求的模数要是任意正数,所以可以求得

    [egin{cases}Ansequiv A_1(mod P_1) \Ansequiv A_2(mod P_2) \...\Ansequiv A_n(mod P_n) \end{cases} ]

    然后用CRT合并答案,求得(Ansequiv OUTANS (mod Mod))

    但是就所见到的题目来说(任意模数NTT),只有在答案不太大的时候适用,例如该题答案不会超过(10^9*10^9*10^5=10^{23}),方法是选取三个乘积大于(10^{23})的NTT模数,得到

    [egin{cases}Ansequiv A_1(mod P_1×P_2) \Ansequiv A_2(mod P_3) \end{cases} ]

    (M=P_1*P_2)于是(Ans=kM+A_1=k_3P_3+A_2),即(kM+A_1equiv A_2(mod P_3))

    由于(Ansle10^{23}<P_1P_2P_3=MP_3),所以(k<P_3),根据此同余方程可以求得(k)(mod P_3)的解也就是(k)的实际值

    然后就可以带入求得答案(Ans=kM+A_1)了,不过会爆(long long),可以采用这个

    ll mul(ll x,ll y,ll p) {return (x*y-(ll)((long double)x/p*y+0.5)*p+p)%p;}
    

    BSGS&EXBSGS

    快速((O(sqrt n)))求解(A^x=B(mod P))(x)的解,板子题:[SDOI2013]随机数生成器

    普通的BSGS要求P是质数,拓展的可以不用是质数

    通过降幂公式,能够知道(xle P-1),令(M=sqrt P),则(x=iM+j),可以枚举(i)的值,得到(t=A^{(i+1)M}),查表看是否存在(T=A^{M-j}B)满足条件,得到解就返回,这样就能保证解是最小的了

    int BSGS(int A,int B,int P)
    {
      	if(!A%P) return -1;
      	int M=sqrt(P)+1;Hash.reset();
      	for(int i=0,t=B;i<M;i++,t=1ll*t*A%P) Hash.Add(t%Mo,t,i);//存入B*A^i的哈希值
      	for(int i=1,bs=ksm(A,M,P),t=bs;i<=M;i++,t=1ll*t*bs%P)//t=A^{(i+1)M}
        	if(Hash.Query(t)!=-1) return i*M-Hash.Query(t);
      	return -1;
    }
    

    FFT/NTT/MTT/FWT

    板子见最下方

    原根相关

    定义:

    原根类似于FFT的单位复根,使得能够进行取模操作

    NTT模数是指有原根且是(r×2^k+1)的形式的模数,如(998244353(3))(1004535809(3))

    存在性及判定

    一个数有原根当且仅当它为(2,4,p,2p,p^r)(p为奇质数)

    (g^1,g^2...g^{phi(p)})(mod p)下各不相同,且(g^{phi (p)}=1),当p为质数时,也就是说(g^1...g^{p-1})分别对应(1...p-1)

    (p)为质数时,令(p-1=p_1^{a1}p_2^{a_2}...p_n^{a_n}),若存在(g^{frac{p-1}{p_i}}=1)(g)不是原根,否则是原根(证明

    注意点

    • 两个多项式最高次项分别为(l1,l2),那么需要开的长度(包括0)是(>l1+l2)的第一个形如(2^k)的数

    组合公式

    • ((0,0))走到((n,m)),不碰到直线(y=x+b)的方案数:(C(n+m,n)-C(n+m,n-b))

      意义为沿着直线翻折,越过直线的路径对称域从翻着顶点到终点的路径

    斯特林数

    第一类斯特林数

    第一类斯特林数:(n)个人坐(m)张圆桌的方案数(人不同,圆桌相同,n个元素构成m个圆排列)

    [s[n,m]=s[n-1,m-1]+s[n-1,m]*(n-1) ]

    第二类斯特林数

    第二类斯特林数:(n)个球放入(m)个盒子的方案数(球不同,盒子相同,n个元素分成m个非空集合)

    [S{n,m}=S{n-1,m-1}+S{n-1,m}*m ]

    卡塔兰数

    [C_n=frac{C_{n-1}*(4n-2)}{n+1}$$$$C_n=C(2n,n)-C(2n,n-1) ]

    常用数学公式

    • 平方和的前缀和 (sum_{i=1}^{n}i^2=frac{n(n+1)(2n+1)}{6})
    • 立方和的前缀和 (sum_{i=1}^{n}i^3=(frac{n(n+1)}{2})^2)

    技巧经验

    容斥

    • ((0,0))((n,m))且不经过一些禁点的方案数(BZOJ两双手)

      (f[i]=sum_{j=0}^{j<i}-f[j]*way(j,i),f[0]=-1,way(j,i))表示从j点到i点的方案数,按照从((0,0))走到该点的所需步数排序

    组合计数

    • 可以画成二维平面上的点去移动,辅助理解、推式子

    区间筛

    • 区间筛质数个数

      luogu1835素数密度:筛出长度为(10^6)的一段区间内的质数个数,左右端点在int以内

      方法:枚举到(sqrt r)然后用(j=(l-1)/i*i+i)暴力在这个区间内筛掉合数

    • 区间筛(varphi())之和

      luogu3601签到题:筛出长度为(10^6)的一段区间的(varphi(i))之和,左右端点在(10^{12})以内

      方法:枚举到(sqrt r)然后(ans[i]=ans[i]/p*(p-1))

    • 区间筛(mu())之和

      luogu3653小清新数学题:筛出长度为(10^6)的一段区间的(mu(i))之和,左右端点在(10^{18})以内

      方法:枚举到(10^6)然后用质数暴力搞,除剩下的数如果不是(1)那么只有三种情况:

    • 质数的平方:用(sqrt)判掉

    • 一个质数

    • 两个质数之积
      上述两种可以用Miller—Rabin素性判断判掉,不过更方便的是(check):(2^{p-1}\%p==1)(3^{p-1}\%p==1)。这题判两次可以过,当然是有可能不准确的。。

    博弈

    • Nim游戏:异或和为0则先手必败,否则必胜
    • 威佐夫博弈:约定(a<b)(a=frac{sqrt 5+1}{2}(a-b))时是奇异局势,先手必败。

    有趣的式子

    gcd有关

    • (sum_{i=1}^{n}gcd(i,n))

      (=sum_{d=1}^{n}(dsum_{i=1}^{n}[gcd(i,n)==d]))

      (=sum_{d=1}^{n}dvarphi(frac{n}{d}))(当然到这步就可以(O(nsqrt n))做了)

      (=sum_{d=1}^{n}frac{n}{d}varphi(d))

      (=sum_{d=1}^{n}frac{n}{d}dprod_{p_i|d}frac{p_i-1}{p_i})

      (=nsum_{d=1}^{n}prod_{p_i|d}frac{p_i-1}{p_i})

      考虑每个质数选或者不选的生成函数:((b_ifrac{p_i-1}{p_i}+1)),其中(b_i)表示(p_i)的指数,含义为选了这个质数就会有(b_i)种指数

      所以最后答案就是(nprod_{p_i|n}(b_ifrac{p_i-1}{p_i}+1))

    数论模板库

    黑科技

    (long long)相乘取模

    ll mul(ll x,ll y,ll P) {return (x*y-(ll)((long double)x/P*y+0.5)*P+P)%P;}
    

    子集枚举

    (0)~(2^n)的子集个数之和是(3^n)(一个位置可以有三种情况:没被选到、选到但是没有被子集枚举到、选到并被枚举到)

    以下可以快速枚举子集

    for(int i=s;;i=(i-1)&s)
    {
    	cout<<i<<endl;
    	if(!i) break;
    }
    

    高维前缀和

    统计子集

    for(int p=1;p<=1<<20;p<<=1)
    	for(int i=0;i<1<<21;i++)
    		if(i&p) f[i]+=f[i^p];
    

    统计超集

    for(int p=1;p<=1<<20;p<<=1)
    	for(int i=0;i<1<<21;i++)
    		if(!(i&p)) f[i]+=f[i|p];
    

    各种线性筛

    用接近线性的方法筛取一些函数,首先要先找到其值与质数的关联才能做
    积性函数都能线性筛

    线性筛质数
    【模板】线性筛素数

    void Get_pri()
    {
    	ispri[1]=1;
    	for(int i=2;i<=N;i++)
    	{
    		if(!ispri[i]) pri[++tot]=i;
    		for(int j=1;j<=tot&&i*pri[j]<=N;j++)
    		{
    			ispri[i*pri[j]]=1;
    			if(i%pri[j]==0) break;
    		}
    	}
    }
    

    线性筛(mu)

    void Get_mu()
    {
    	ispri[1]=mu[1]=1;
    	for(int i=2;i<=N;i++)
    	{
    		if(!ispri[i]) pri[++tot]=i,mu[i]=-1;
    		for(int j=1;j<=tot&&i*pri[j]<=N;j++)
    		{
    			ispri[i*pri[j]]=1;
    			if(i%pri[j]) mu[i*pri[j]]=-mu[i];
    			else break;
    		}
    	}
    }
    

    线性筛(varphi)

    void Get_phi()
    {
    	ispri[1]=phi[1]=1;
    	for(int i=2;i<=N;i++)
    	{
    		if(!ispri[i]) pri[++tot]=i,phi[i]=i-1;
    		for(int j=1;j<=tot&&i*pri[j]<=N;j++)
        	{
    			ispri[i*pri[j]]=1;
    			if(i%pri[j]) phi[i*pri[j]]=phi[i]*(pri[j]-1);
    			else {phi[i*pri[j]]=phi[i]*pri[j];break;}
    		}
    	}
    }
    

    高级算法

    Exgcd

    上文有讲用法,注意传入的时候要求都除以了(gcd)

    void Exgcd(int a,int b,int &x,int &y)
    {
        if(!b) {x=1;y=0;return;}
        Exgcd(b,a%b,y,x);y-=a/b*x;
    }
    

    Lucas

    【模板】卢卡斯定理
    用于计算模数比较小的时候的组合数,公式看就是(C(n,k)=C(n/p,k/p)*C(n\%p,k\%p))
    然后继续递归,边界条件是(k)(0)时答案是(1)

    int C(int n,int k) {return 1ll*jc[n]*inv[k]%P*inv[n-k]%P;}
    int Lucas(int n,int k) {return !k?1:1ll*Lucas(n/P,k/P)*C(n%P,k%P)%P;}
    int pre()
    {
        inv[0]=inv[1]=jc[0]=1;
        for(int i=2;i<P;i++) inv[i]=(P-1ll*P/i*inv[P%i]%P)%P;
        for(int i=1;i<P;i++) jc[i]=1ll*jc[i-1]*i%P,inv[i]=1ll*inv[i]*inv[i-1]%P;
    }
    

    EXCRT

    【模板】扩展中国剩余定理(EXCRT)
    合并同余方程,复杂度瓶颈在于exgcd

    ll mul(ll x,ll y,ll p) {return (x*y-(ll)((long double)x/p*y+0.5)*p+p)%p;}
    ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
    void Exgcd(ll a,ll b,ll &x,ll &y) {!b?(x=1,y=0):(Exgcd(b,a%b,y,x),y-=a/b*x);}
    ll n,B1,P1,B2,P2;
    int main()
    {
    	cin>>n>>P1>>B1;n--;
    	while(n--)
    	{
    		cin>>P2>>B2;
    		ll g=gcd(P1,P2),C=(B2-B1)/g,x,y,P;
    		Exgcd(P1/g,P2/g,x,y); P=P1/g*P2;
    		B1=(mul(mul(x,C,P),P1,P)+B1+P)%P; P1=P;
    	}
    	cout<<B1<<endl;
    }
    

    BSGS

    [CQOI2018]破解D-H协议
    (g^a equiv A(mod P)),给定(A),求(a)
    方法:

    • (a=im-j)所以(g^{im}equiv A*g^j),其中(m=sqrt P)能够保证最优复杂度
    • 然后从(0)(m)枚举(i),把(g^{im})存进哈希表
    • 暴力判断对于每一个(j),哈希表里是否有(A*g^j)的值

    复杂度(O(sqrt PlogP))

    const int Mo=100003;
    int g,A,a,b,P,n;
    int ksm(int x,int k)
    struct HashTable
    {
    	struct Line{int next,v1,v2;}a[Mo];
    	int head[Mo],cnt;
    	void reset() {memset(head,0,sizeof(head));cnt=0;}
    	void Add(int p,int v1,int v2) {a[++cnt]=(Line){head[p],v1,v2};head[p]=cnt;}
    	int Query(int x)
    		{
    			for(int i=head[x%Mo];i;i=a[i].next)
    				if(a[i].v1==x) return a[i].v2;return -1;
    		}
    }Hash;
    int BSGS(int g,int A,int P)
    {
    	int M=sqrt(P)+1;Hash.reset();
    	for(int i=0,t=A;i<M;i++,t=1ll*t*g%P) Hash.Add(t%Mo,t,i);
    	for(int i=1,bs=ksm(g,M),t=bs;i<=M;i++,t=1ll*t*bs%P)
    		if(Hash.Query(t)!=-1) return i*M-Hash.Query(t);return -1;
    }
    

    高斯消元

    【模板】高斯消元法
    (O(n^3))求解线性方程组。用线性基的方式实现

    const int N=110;
    const double eps=1e-7;
    int n;
    double F[N][N],A[N];
    int main()
    {
    	cin>>n;
    	for(int w=1;w<=n;w++)
    	{
    		for(int i=1;i<=n+1;i++) cin>>A[i];
    		for(int i=1;i<=n;i++)
    		{
    			if(fabs(A[i])<eps) continue;
    			if(fabs(F[i][i])<eps)
    			{
    				
    				for(int j=i;j<=n+1;j++) F[i][j]=A[j]/A[i];
    				F[i][i]=1;break;
    			}
    			else for(int j=i+1;j<=n+1;j++) A[j]-=F[i][j]*A[i];
    		}
    	}
    	for(int i=1;i<=n;i++)
    		if(fabs(F[i][i])<eps)
    			return puts("No Solution"),0;
    	for(int i=n;i>=1;i--)
    		for(int j=i-1;j>=1;j--)
    			F[j][n+1]-=F[j][i]*F[i][n+1],F[j][i]=0;
    	for(int i=1;i<=n;i++) printf("%.2f
    ",F[i][n+1]);
    }
    
    

    线性基

    【模板】线性基
    很方便地求解一个数集中能够异或出来的数

    ll n,x,a[60],ans;
    int main()
    {
    	for(cin>>n;n;n--)
    	{
    		cin>>x;
    		for(int i=50;i>=0;i--)
    			if(x&(1ll<<i))
    			{
    				if(!a[i]) {a[i]=x;break;}
    				else x^=a[i];
    			}
    	}
    	for(int i=50;i>=0;i--)
    		if((ans^a[i])>ans) ans^=a[i];
    	return cout<<ans<<endl,0;
    }
    

    裴蜀定理

    【模板】裴蜀定理
    问题:给定(A_1,A_2...A_k),求(A_1x_1+A_2x_2+...A_kx_k)所能取到的最小值
    结论:所有数的绝对值的(gcd)

    int g,x,n;
    int main()
    {
    	cin>>n;while(n--) cin>>x,g=__gcd(g,abs(x));
    	cout<<g<<endl;
    }
    

    FFT

    【模板】多项式乘法(FFT)
    (O(nlogn))求解两多项式卷积

    //FFT
    for(int i=1;i<l;i++) r[i]=(r[i>>1]>>1)|((i&1)<<tt);
    for(int i=0;i<l;i++) w[i].rl=cos(Pi/l*i),w[i].im=sin(Pi/l*i);
    void FFT(Complex *P,int op)
    {
        for(int i=1;i<l;i++) if(r[i]>i) swap(P[i],P[r[i]]);
        for(int i=1;i<l;i<<=1)
            for(int j=0,p=i<<1;j<l;j+=p)
                for(int k=0;k<i;k++)
                {
                    Complex W=w[l/i*k];W.im*=op;
                    Complex X=P[j+k],Y=P[j+k+i]*W;
                    P[j+k]=X+Y;P[j+k+i]=X-Y;
                }
    }
    for(int i=0;i<=m+n;i++) printf("%d ",(int)(A[i].rl/l+0.5));
    

    拉格朗日插值

    (n-1)次多项式带入(x=k)求值。

    [f(k)=sum_{i=1}^{n}y[i]prod_{j!=i}frac{k-x[j]}{x[i]-x[j]} ]

    int Lagrange()
    {
    	cin>>n>>k;
    	for(int i=1;i<=n;i++) cin>>x[i]>>y[i];
    	for(int i=1;i<=n;i++)
    	{
    		int fz=1,fm=1;
    		for(int j=1;j<=n;j++)
    			if(i!=j) fz=1ll*fz*(k-x[j])%P,fm=1ll*fm*(x[i]-x[j])%P;
    		(ans+=1ll*y[i]*fz%P*ksm(fm,P-2)%P)%=P;
    	}
    	cout<<(ans+P)%P<<endl;
    }
    

    NTT

    FWT

    这两个一个用NTT代替,一个用高维前缀和。考到了话我直播CS!

  • 相关阅读:
    LinqToXML~读XML文件
    C#~使用FileSystemWatcher来监视文件系统的变化
    HDU4144:Bacon's Cipher
    国产手机的路还很长
    同一个存储过程中,不能多次select into 到同一张表的问题
    IT行业为什么没有进度
    IOT(Index Organized Table)
    模版——容器,迭代器
    shell split分析日志文件
    VMware vSphere 服务器虚拟化之十七 桌面虚拟化之安装View链接服务器
  • 原文地址:https://www.cnblogs.com/xzyxzy/p/9903896.html
Copyright © 2011-2022 走看看