一、素数
Miller-Rabin
首先介绍一下伪素数:若n是一个正整数,且存在正整数a满足$a^{n-1}equiv1;(mod;n)$
(费马小定理,但n不一定为素数)
那么我们说n是基于a的伪素数
如果一个数是伪素数,它很大概率是素数
但一个数不是伪素数,它一定不是素数
那么对于要判断的数$n$,我们只需要多次选取$a$来判断是否是伪素数即可
然而是否伪素数通过了所有底数的测试就很大概率是素数了呢?
其实有一种叫$Carmichael$数的东西,十亿内有$600$多个,最小的如$561$,可以通过所有底数的测试。
那咋整?我们可以继续探测啊QwQ
二次探测定理:若$a^2~mod~p=1$,且$a!=1$,$a!=-1$则$p$必为合数。
设$p-1=2^ld$,可以先测$a^d$,然后平方$l$次。
当两倍两倍往上加的时候,判断是否合法,如果出现不合法就$return$。判定方法结合定理看代码吧QwQ……
至于伪素数的判定就两倍两倍加回原数之后在弄吧。
代码和$Pollard~Rho$的放一起吧QAQ
Pollard Rho
大数分解……
暂时还不是很懂就不瞎讲了……先坑着吧,不过我写的是第一篇博客里的第二种判环方法QAQ
1 #include<iostream> 2 #include<cstring> 3 #include<cstdio> 4 #include<algorithm> 5 #define LL long long 6 using namespace std; 7 8 LL T,maxn,x; 9 LL prime[10]={2,3,5,7,11,13,17,19,23}; 10 11 LL Mul(LL a,LL b,LL MOD) 12 { 13 LL tmp=a*b-(LL)((long double)a*b/MOD+0.1)*MOD; 14 return tmp<0?tmp+MOD:tmp; 15 } 16 17 LL Qpow(LL a,LL b,LL MOD) 18 { 19 LL ans=1; 20 while (b) 21 { 22 if (b&1) ans=Mul(ans,a,MOD); 23 a=Mul(a,a,MOD); b>>=1; 24 } 25 return ans; 26 } 27 28 LL gcd(LL a,LL b) {return b==0?a:gcd(b,a%b);} 29 30 bool Miller_Rabin(LL n) 31 { 32 if (n==2) return 1; 33 if (n<2 || n%2==0) return 0; 34 LL m=n-1,l=0; 35 while (m%2==0) ++l, m>>=1; 36 for (int i=0; i<9; ++i) 37 { 38 LL p=prime[i],w=Qpow(p,m,n); 39 if (w==1 || w==n-1 || p==n) continue; 40 for (int j=1; j<=l; ++j) 41 { 42 LL u=Mul(w,w,n); 43 if (u==1 && w!=1 && w!=n-1) return 0; 44 w=u; 45 } 46 if (w!=1) return 0; 47 } 48 return 1; 49 } 50 51 LL Pollard_Rho(LL n,LL c) 52 { 53 LL x=rand()%n,y=x,p=1,k=2; 54 for (LL i=1; p==1; ++i) 55 { 56 x=(Mul(x,x,n)+c)%n; 57 p=x>y?x-y:y-x; 58 p=gcd(p,n); 59 if (i==k) y=x,k+=k; 60 } 61 return p; 62 } 63 64 void Solve(LL n) 65 { 66 if (n==1) return; 67 if (Miller_Rabin(n)) {maxn=max(maxn,n); return;} 68 LL t=n; 69 while (t==n) t=Pollard_Rho(n,rand()%(n-1)+1); 70 Solve(t); Solve(n/t); 71 } 72 73 int main() 74 { 75 scanf("%lld",&T); 76 while (T--) 77 { 78 scanf("%lld",&x); 79 maxn=0; 80 Solve(x); 81 if (maxn==x) puts("Prime"); 82 else printf("%lld ",maxn); 83 } 84 }
埃筛,欧拉筛
没什么好说的,讲解网上的很好,直接上code吧
埃筛
1 void work(int n) 2 { 3 for (int i=2;i<=n;++i) 4 if (!vis[i]) 5 { 6 prime[++cnt]=i; 7 for (int j=i*2;j<=n;j+=i) 8 vis[j]=1; 9 } 10 }
欧拉筛
1 void work(int n) 2 { 3 for (int i=2;i<=n;++i) 4 { 5 if (!vis[i]) 6 prime[++cnt]=i; 7 for (int j=1;j<=cnt&&prime[j]*i<=n;++j) 8 { 9 vis[prime[j]*i]=1; 10 if (i%prime[j]==0) 11 break; 12 } 13 } 14 }
二、欧拉函数
欧拉函数$φ(n)$表示1~n中与n互质的个数(注意1和任何数互质)
线筛欧拉函数
1 void Euler(long long n) 2 { 3 phi[1]=1; 4 for (int i=2;i<=n;++i) 5 if (!phi[i]) 6 for (int j=i;j<=n;j+=i) 7 { 8 if (!phi[j]) phi[j]=j; 9 phi[j]=phi[j]/i*(i-1); 10 } 11 }
欧拉定理:$a^{φ(n)}equiv1;(mod;n)$
三、中国剩余定理(CRT)
CRT
中国剩余定理,就是用来求解一组同余方程组的一个解的,其中模数$m_{i}$两两互质
引入一个经典问题
“今有物不知其数,三三数之剩二(除以3余2),五五数之剩三(除以5余3),七七数之剩二(除以7余2),问物几何?”这个问题称为“孙子问题” (一个很妙的证明)
设$q_{1}$为$3*k+2$,也就是除3余2的一个数
同理设$q_{2}$为$5*k+3$,$q_{3}$为$7*k+2$
那么可以导出以下三点:
1、要使$q_{1}+q_{2}+q_{3}$的和满足除3余2,$q_{2}和q_{3}$必须是3的倍数
2、要使$q_{1}+q_{2}+q_{3}$的和满足除5余3,$q_{1}和q_{3}$必须是5的倍数
3、要使$q_{1}+q_{2}+q_{3}$的和满足除7余2,$q_{1}和q_{2}$必须是7的倍数
整理可得
1、$q_{1}$除以3余2,且是5和7的公倍数。
2、$q_{2}$除以5余3,且是3和7的公倍数。
3、$q_{3}$除以7余2,且是3和5的公倍数。
显然$q_{1}+q_{2}+q_{3}$是符合条件的一组解,可以用exgcd进行求解了
这里有一个小技巧,比如求$q_{1}$的时候可以求出除3余1的数然后再乘2,
也就是先求出5和7的公倍数模3下的逆元,再用逆元去乘余数。
看代码可能更好懂一点
1 long long CRT() 2 { 3 long long x,y,ans=0,lcm=1; 4 for(i=1; i<=n; ++i) lcm=lcm*m[i];//因为互质 5 for(i=1; i<=n; i++) 6 { 7 long long kl=lcm/m[i]; 8 exgcd(kl,m[i],x,y); 9 x=(x%m[i]+m[i])%m[i]; 10 ans=(ans+a[i]*x*kl)%lcm; 11 } 12 return ans; 13 }
exCRT
用来解决模数不两两互质的同余方程组
假设我们这里有两个方程
$x = a_1 * x_1 + b_1$
$x = a_2 * x_2 + b_2$
$a_1,a_2$是模数,$b_1,b_2$是余数,那么我们可以合并这两个方程:
因为正负号对变量$x_1 , x_2$无影响,所以稍加变形可得:$a_1 * x_1 + a_2 * x_2 = b_2 - b_1$
然后扩欧就可以求解出最小正整数解$x_1$,那么令$k=(a_1*x_1+b_1)$
因为$k=a_1*x_1+b_1=a_2*x_2+b_2$,显然 $k$ 和$a_1 , a_2$分别同余 ,那么$ k $ 和 $ lcm( a_1 , a_2 ) $一定同余
所以 $x equiv k ; (mod ; lcm( a_1 , a_2 ) )$
不停的合并就可以得到最终答案了
1 long long exCRT()//M是模数,A是余数 2 { 3 long long M=m[1],A=a[1],d,x,y; 4 for (int i=2; i<=n; ++i) 5 { 6 exgcd(M,m[i],d,x,y); 7 if ((a[i]-A)%d) return -1; 8 x*=(a[i]-A)/d; t=m[i]/d; x=(x%t+t)%t; 9 A=M*x+A; M=M/d*m[i]; A%=M; 10 } 11 A=(A%M+M)%M; 12 return A; 13 }
四、Lucas
Lucas
证明不会,会用就行
Lucas用于解决当模数p为质数且比较小的时候,可以用$log$的复杂度解决求解$C(n,m)$
$Lucas(n,m) = Lucas(n/p,m/p)*C(n\%p,m\%p)$
其中$C(n\%p,m\%p)$可以通过预处理阶乘和阶乘的逆元$O(1)$进行求解
1 void Init(long long p) 2 { 3 fac[0]=1; inv[1]=1; facInv[0]=1; 4 for (int i=1; i<=p; ++i) 5 { 6 if (i!=1) inv[i]=(p-p/i)*inv[p%i]%p; 7 fac[i]=fac[i-1]*i%p; facInv[i]=facInv[i-1]*inv[i]%p; 8 } 9 } 10 11 long long C(long long n,long long m) 12 { 13 if (m>n) return 0; 14 return fac[n]*facInv[m]%MOD*facInv[n-m]%MOD; 15 } 16 17 long long Lucas(long long n,long long m) 18 { 19 if (m>n) return 0; 20 long long ans=1; 21 for (; m; n/=MOD,m/=MOD) 22 ans=(ans*C(n%MOD,m%MOD))%MOD; 23 return ans; 24 }
exLucas
可以解决$C(n,m)\%p$中模数不是质数的Lucas问题
首先我们将模数$p$分解,$p={p_1}^{k_1}*{p_2}^{k_2}*...*{p_q}^{k_q}$
很容易得出一组同余方程组
$left{ egin{array}{c} ansequiv c_1pmod {{p_1}^{k_1}}\ ansequiv c_2pmod {{p_2}^{k_2}}\ ...\ ansequiv c_qpmod {{p_q}^{k_q}}\ end{array} ight.$
很显然模数是两两互质的,只要求出每一个$c_i$那么我们就可以用中国剩余定理求解出$C(n,m)$了。
那么问题来了,怎么求解每一个$c_i$呢?
因为$C_n^m={n!over m!(n-m)!}$,所以只要能求出 $n!~\% {p_i}^{k_i},~m!~\% {p_i}^{k_i},~(n-m)!~\% {p_i}^{k_i}$,我们就可以根据逆元求解每一个$c_i$了。
举个栗子
$n=22,p_i=3,k_i=2$
$22!=(1*2*3*4*5*6*7*8*9)*(10*11*12*13*14*15*16*17*18)*(19*20*21*22)$
剔除因子3之后
$22!=(1*2*4*5*7*8)*(10*11*13*14*16*17)*(19*20*22)*3^6*(1*2*3*4*5*6*7)$
发现按照如下分组后$3^6$前面部分的前$frac{n}{{p_i}^{k_i}}$组在$mod~{p_i}^{k_i}$ 后结果相同,因此只需要算出来一组然后快速幂就好了。
对于冗余部分$(19*20*22)$暴力计算。
最后的$(1*2*3*4*5*6*7)$递归计算即可
3因子被剔除了就不考虑了
代码还写不来先坑着
未完待续目前懒得写了