【算法竞赛学习笔记】数论（三）-数学提升计划

---

title : 数论笔记（三）
date : 2021-8-12
tags ： 数论,ACM

---

快速乘

龟速乘

事实上快速乘是为了防止溢出，又不想写高精度，所以我们模仿二进制加法来完成两数的取模乘积。复杂度$O(logn)$

ll mul(ll x,ll y,ll mod){ //=>x*y%mod
    ll res=0;
    while(y){
        if(y&1) res=(res+x)%mod;
        x=(x+x)%mod;
        y>>=1;
    }
    return res;
}//其中ll是long long类型

优秀的long double快速乘

先上代码，复杂度只有$O(1)$噢

ll mul(ll a,ll b,ll mod){
    ull c=(ull)a*b-(ull)((ld)a/mod*b+0.5L)*mod;
} //其中ull是unsigned long long，ld是long double类型

仔细一看，这直接用了乘法操作不就爆掉了吗？

其实a*b和(a*b/p)*p两部分都是会溢出的，但是unsigned保证了他们溢出后的差值不变，因此不会影响最终结果。（反正很巧妙）

判断素数

范围较小的时候，可以从2~sqrt(n)中判断有无n的因子，

范围较大且要找出所有素数时，可以使用素数筛。

Miller Rabin 素数测试

前置知识：唯一分解定理、威尔逊定理、费马定理

威尔逊定理：$若p为素数，则(p-1)!\equiv -1 \mod p $

威尔逊逆定理：$若(p-1)!\equiv -1\mathrm{m}\mathrm{o}\mathrm{d}p，则p一定为素数$

费马定理：$若p为素数，a为正整数，且a和p互质，则a^{p-1}\equiv 1\mathrm{m}\mathrm{o}\mathrm{d}p$

测试过程

$$\begin{gathered} (1)计算奇数M，使得N=2^r\ast M+1 \\ (2)选择随机数A<N \\ (3)对于任意i<r，若A^{2^i\ast M}\mathrm{m}\mathrm{o}\mathrm{d}N=N-1,则N通过随机数A的测试 \\ (4)或者，若A^M\mathrm{m}\mathrm{o}\mathrm{d}N=1,则N通过随机数A的测试 \\ (5)让A取不同的值对N进行5次测试，若全部通过则判定N为素数 \end{gathered}$$

概括

不断选取不超过n-1的基b（共s次），计算是否每次都有$b^{n-1}\equiv 1\mathrm{m}\mathrm{o}\mathrm{d}n$，若每次都成立则n是素数，否则为合数。

二次探测定理

$如果p是奇素数，则x^2\equiv 1\mathrm{m}\mathrm{o}\mathrm{d}p的解为x=1||x=p-1\mathrm{m}\mathrm{o}\mathrm{d}p$

这个定理可以提高测试效率。

#include <bits/stdc++.h>
#define int long long
#define ull unsigned long long
#define ld long double
using  namespace std;

const int p[]={2,3,5,7,11,13,17,19,61}; //越多越准确 

int mul(int a, int b, int p){ //快速乘 
	ull c=(ull)a*b-(ull)((ld)a/p*b+0.5L)*p; //自动溢出 
	return c<p?c:c+p; //转化为正数 
}

int fpow(int a, int b, int mod){ //快速幂里的乘法用快速幂替换 
    int res = 1;
    while(b){
    	if(b&1) res=mul(res,a,mod);
    	a=mul(a,a,mod);
    	b>>=1;
	}
    return res;
}

bool miller_rabin(int x){ //米勒拉宾质数判定 
	if (x < 2) return 0;
	int i,j,y=x-1,a,b;
	while(y&1^1) y>>=1; //相当于先化为奇数 
	for(i=0;i<9;i++){
		if (x%p[i]<1)	return x == p[i]; //不可以是质数p的倍数 
	}
	for (i=2;i<9;i++){
		for (a=b=fpow(p[i],y,x),j=y;j<x&&a>1;j+=j,a=b){
			if (b=mul(b,b,x),b==1&&a!=x-1)	return 0; //二次检测 
			//b和x互质，b^2modx=1，要b=1,要么b=x-1,那么b==1时a应该为x-1 
		}
		if(a!=1) return 0; //最终a不为1，没有通过二次检验 
	}
	return 1;
}

signed main(){
	int n;
	while(cin>>n){ 
		puts(miller_rabin(n)?"Y":"N");
	}
	return 0;
}

素数筛

暴力枚举求素数

枚举每个数，看是否有正整数能整除这个数。$复杂度O(n\sqrt{n})$​

for(int i=2;i<=n;i++){
    bool g=0;
    for(int j=2;j*j<=i;j++)
        if(i%j==0){g=1;break;}
    if(!g){
    	tot++;
        p[tot]=i;
    }
}

埃氏筛（普通筛法）

每个合数可以分解为素数的乘积，那么在搜索到一个数为素数的时候，就把它的倍数标记为合数。$复杂度O(nlognlogn)$

for(int i=2;i<=n;i++){
    if(prime[i]==0){
        p[++tot]=i;
        for(int j=2;j*i<=n;j++){
            prime[j*i]=1;
        }
    }
}

欧拉筛 （线性筛法）

每次用最小质因子来筛素数，避免重复筛选。确保每个合数只被最小质因子p筛一次。

for(int i=2;i<=n;i++){
    if(prime[i]==0) p[++tot]=i;
    for(int j=1;j<=tot&&i*p[j]<=n;j++){
        prime[i*p[j]]=1;
        if(i%p[j]==0) break;
    }
}

欧拉函数

对正整数n，欧拉函数是小于等于n的数中与n互质的数的数目。

欧拉函数又称为$\varphi$函数，例如$\varphi (8)=4$，因为1，3，5，7均与8互质。

引理
$$\begin{gathered} ①如果n为某一个素数p，则：\varphi (p)=p-1 \\ ②如果n为某一个素数p的幂次p^a，则\varphi (p^a)=(p-1)\ast p^{a-1} \\ ③如果n为任意两个互质的数a,b的积，则\varphi (a\ast b)=\varphi (a)\ast \varphi (b) \end{gathered}$$

证明

$$\begin{gathered} ②共p^a-1个数比p^a小，其中有p^{a-1}-1个数能被p整除，表示为p\ast t(t=1,2,...,p^{a-1}-1),相减可得 \\ ③在比a\ast b小的a\ast b-1个整数中，有\varphi (a)个与a互质的数，有\varphi (b)个与b互质的数，必须既与a互质，又与b互质，才会与a\ast b互质，满足条件的数共有\varphi (a)\ast \varphi (b)个。 \end{gathered}$$

暴力

int Eular(int m){
	int ret=m;
	for(int i=2;i<m;i++){
		if(m%i==0) ret-=ret/i;
		while(m%i==0){
			m/=i;
		}
	}
	if(m>1) ret-=ret/m;
	return ret;
} 

筛选法

void init(){ //筛选法打欧拉函数表 
	euler[1]=1;
	for(int i=2;i<maxn;i++) euler[i]=i;
	for(int i=2;i<maxn;i++)
		if(euler[i]==i)
			for(int j=i;j<maxn;j+=i)
				euler[j]=euler[j]/i*(i-1);
}

直接求解欧拉函数

int eu(int n){  //直接求解欧拉函数 
        int res=n,a=n;  
        for(int i=2;i*i<=a;i++){
                if(a%i==0){  //分解质因数 
                        res=res/i*(i-1);
                        while(a%i==0) a/=i; //把质因数彻底分解 
                }
        }
        if(a>1) res=res/a*(a-1);
        return res;
}

线性（重要）

void Get_phi(){
	cnt=0;
	memset(flag,1,sizeof(flag));
	phi[1]=1;
	for(int i=2;i<maxn;i++){
		if(flag[i]){
			p[cnt++]=i;
			phi[i]=i-1;
		}
		for(int j=0;j<cnt;j++){
			if(i*p[j]>maxn) break;
			flag[i*p[j]]=false;
			if(i%p[j]==0){
				phi[i*p[j]]=p[j]*phi[i];
				break;
			}else phi[i*p[j]]=(p[j]-1)*phi[i];
		}
	}
}	

需要注意的性质：

①phi§==p-1是因为素数p除了1以外的因子只有p，所以与p互质的个数是p-1个；

②$phi(p^k)==p^k-p^{k-1}==(p-1)\ast p^{k-1}$

Pollard Rho算法求大数因子

时间效率要求较高时，用来分解一个合数n。

luoguP4718 【模板】Pollard-Rho算法

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define lll __int128
ll max_factor,t,n;
inline ll read(){	ll x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-') f=f*-1;ch=getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}
inline ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
inline ll qp(ll x,ll p,ll mod){
    ll ans=1;
    while(p){
        if(p&1) ans=(lll)ans*x%mod;
        x=(lll)x*x%mod;
        p>>=1;
    }
    return ans;
}

inline bool mr(ll x,ll b){ //米勒罗宾素数检验 
    ll k=x-1;
    while(k){
        ll cur=qp(b,k,x);
        if(cur!=1&&cur!=x-1) return false;
        if((k&1)==1||cur==x-1) return true;
        k>>=1;
    }
    return true;
}

inline bool prime(ll x){
    if(x==46856248255981ll|| x<2) return false;
    if(x==2||x==3||x==7||x==61||x==24251) return true;
    return mr(x,2)&&mr(x,61);
}

inline ll fc(ll x,ll c,ll n){return ((lll)x*x+c)%n;}

inline ll PR(ll x){ //Pollard Rho算法求素因子 
    ll s=0,t=0,c=1ll*rand()%(x-1)+1;
    int stp=0,goal=1;
    ll val=1;
    for(goal=1;;goal<<=1,s=t,val=1){
        for(stp=1;stp<=goal;++stp){
            t=fc(t,c,x);
            val=(lll)val*abs(t-s)%x;
            if((stp%127)==0){
                ll d=gcd(val,x);
                if(d>1) return d;
            }
        }
        ll d=gcd(val,x);
        if(d>1) return d;
    }
}

inline void fac(ll x){ //求最大因子 
    if(x<=max_factor||x<2) return;
    if(prime(x)){
        max_factor=max_factor>x?max_factor:x;
        return;		
    }
    ll p=x;
    while(p>=x) p=PR(x);
    while((x%p)==0) x/=p;
    fac(x),fac(p);
}

signed main(){
	cin>>t;
    srand((unsigned)time(NULL));
	while(t--){
        n=read();
        max_factor=0;
        fac(n);
        if(max_factor==n) puts("Prime");
        else printf("%lld\n",max_factor);
    }
    return 0;
}

参考资料

信息学奥赛之数学一本通 （林厚从）

https://www.cnblogs.com/zjp-shadow/p/9267675.html#autoid-3-3-0