原根,BSGS,扩展BSGS,miller_rabbin算法,高斯消元 证明及模板

原根

如果两个整数\(a,b\)互质,则有\(a^{\phi(b)}\%b=1\)
定义模\(b\)意义下的\(a\)的阶为使\(a^d\%b=1\)的最小正整数\(d\)
显然,模\(b\)的阶\(d|\phi(b)\)
如果模\(b\)意义下\(a\)的阶为\(\phi(b)\),则称\(a\)\(b\)原根

欧拉证明了(我证明不了)
在这里插入图片描述

\(b\)存在原根的充要条件为:\(b=2,4,p^n,2p^n\),其中\(p\)为奇素数,\(n∈N^*\)
当模\(b\)有原根时,其原根个数为\(\phi(\phi(b))\)

那么如何求一个数的原根?
首先判断它是否有原根
如果有,一般最小原根都很小,可以选择美丽的暴力

预处理出\(\phi(b)\)的所有因子,存放在数组\(D\)
\([2,\phi(b))\)区间中,从小到大米枚举变量\(i\)
如果存在\(j∈D\),使得\(i^{\frac{\phi(b)}{j}}\%b=1\),则\(i\)不是原根,否则\(i\)是原根

为什么这样是正确的呢?
因为如果存在一个数\(1\le x<\phi(b)\),使得\(i^x\%b=1\)
\(x|\phi(b)\),并且一定存在一个\(\phi(b)\)的质因子\(j\),使得\(x|\frac{\phi(b)}{j}\)


BSGS大步小步算法

BSGS是用来解决离散对数问题
\(a^x\equiv b\ (mod\ \ p)\),其中\(a,b,p\)已知,且\((a,p)=1\),求\(x\)

因为\(a^{\phi(p)}\equiv 1\ (mod\ \ p),\phi(p)<p\)
所以可以采取枚举法,在\(O(p)\)的时间复杂度求出\(x\)

\(BSGS\)可以在\(O(\sqrt{p})\)的时间复杂度内求出\(x\)

\(m=\lceil \sqrt{p} \rceil,r=x\mod m\),则\(x=k*m+r\),其中\(0\le k<m,0\le r<k\),有

\[a^{k*m+r}\equiv b\ \ (mod\ \ p) \]

因为\((a,p)=1\),方程两边同乘\(a^{-r}\)

\[a^{k*m}\equiv b*a^{-r}\ \ (mod\ \ p) \]

对于\(0\le i<r\),求出所有的\(b*a^{-1}\mod\ \ p\),将其值以及\(i\)存入\(map\)
然后再求左边的\(a^{j*m}\mod p,(0\le j<k)\),并在\(map\)里寻找是否出现过相同的值
如果有,代表着同余,已经找到了答案,如果没有则是无解

然而。。。
可以稍微转变一下算法过程,避免求逆元的操作
在这里插入图片描述

\[x=k*m+r\ ——>x=k*m-r,(1\le k\le m,0\le r<m) \]

则$$a^{k*m-r}\equiv b\ \ (\mod p)$$
两边同时乘以\(a^r\)

\[a^{k*m}\equiv b*a^r\ \ (\mod p) \]

先求出右边所有的\(b*a^i(\mod p)(1\le i\le m)\)保存在\(map\)
然后再求左边的\(a^{j*m}\mod\ p\),并在\(map\)中查找是否出现过

如果出现过,左边枚举的是\(j\),右边枚举的是\(i\),则答案为\(x=j*m-i\),这样就避免了求逆元的操作,却仍然用到了逆元
因为推导必须是等价推导,只有当\((a,p)=1\),即\(a^r\)逆元存在时,才可以大胆两边同乘\(a^r\)等价,因为如果\((a,p)≠1\),式子倒推不回去

#include <map>
#include <cmath>
#include <cstdio>
using namespace std;
#define int long long
map < int, int > mp;
int a, mod, r;

int bsgs() {
	mp.clear();
	int m = ceil( sqrt( mod ) );
	int tmp = 1;
	for( int i = 1;i <= m;i ++ ) {
		tmp = tmp * a % mod;
		mp[tmp * r % mod] = i;
	}
	int res = 1;
	for( int i = 1;i <= m;i ++ ) {
		res = res * tmp % mod;
		if( mp[res] ) return m * i - mp[res];
	}
	return -1;
}

signed main() {
	int ans;
	while( ~ scanf( "%lld %lld %lld", &mod, &a, &r ) ) {
		r %= mod, a %= mod;
		if( r == 1 ) ans = 0;
		else if( ! a ) {
			if( r ) ans = -1;
			else ans = 1;
		}
		else ans = bsgs();
		if( ans == -1 ) printf( "no solution\n" );
		else printf( "%lld\n", ans );
	}
	return 0;
}

扩展BSGS

对于\(a^x\equiv b(\mod p)\)
如果\((a,p)>1\),则无法直接套用BSGS,此时就要用到扩展BSGS

将要求解的式子变形

\[a^x+k*p=b \]

\((a,p)=g\),若\(b\)不是\(g\)的倍数,则方程无解
不过有一个例外\(b=1,x=0\),这个情况特判即可

式子左右两边除以\(g\)

\[a^{x-1}\frac{a}{g}+k\frac{p}{g}=\frac{b}{g} \]

\(a'=\frac{a}{g},p'=\frac{p}{g},b'=\frac{b}{g}\)

\[a^{x-1}a'+kp'=b' \]

\(a,\frac{p}{g}\)仍然不互质,则继续以上操作找出\(a,p'\)的最大公约数\(g'\)
最新式子两边继续除以\(g'\),直到\((a,p')=1\)为止

在此过程中,假设取出来了\(cnt\)\(a\),出去公约数后剩下的乘积为\(a'\)
此时\((a',p')=1\),于是可以转化为$$a^{x-cnt}\equiv b'(a')^{-1}\ (\mod p)$$
其中\(cnt\)表示两边除以最大公约数\(g\)的次数

此处右边有逆元,为了避免求逆元,将\(a'\)保留在左边
在BSGS枚举左边时,初始值设为\(a'\)即可

如果在除以\(g\)的过程中,发现\(b'(a')^{-1}=1\),则立马可以得到答案
\(x-cnt=0\ ——>x=cnt\)

接下来,直接套用基础BSGS即可,记得最后的答案不要忘记加上\(cnt\)哟(^U^)ノ~YO
在这里插入图片描述

#include <map>
#include <cmath>
#include <cstdio>
using namespace std;
#define int long long
map < int, int > mp;
int a, mod, r;

int gcd( int x, int y ) {
	if( ! y ) return x;
	else return gcd( y, x % y );
}

int qkpow( int x, int y, int mod ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

int exbsgs() {
	if( r == 1 ) return 0;
	int cnt = 0, tmp = 1, d;
	while( 1 ) {
		d = gcd( a, mod );
		if( d == 1 ) break;
		if( r % d ) return -1;
		r /= d, mod /= d;
		tmp = tmp * ( a / d ) % mod;
		++ cnt;
		if( tmp == r ) return cnt;
	}
	mp.clear();
	int res = r;
	mp[r] = 0;
	int m = ceil( sqrt( 1.0 * mod ) );
	for( int i = 1;i <= m;i ++ ) {
		res = res * a % mod;
		mp[res] = i;
	}
	int temp = qkpow( a, m, mod );
	res = tmp;
	for( int i = 1;i <= m;i ++ ) {
		res = res * temp % mod;
		if( mp.count( res ) )
			return i * m - mp[res] + cnt;
	}
	return -1;
}

signed main() {
	while( ~ scanf( "%lld %lld %lld", &a, &mod, &r ) ) {
		if( ! a && ! mod && ! r ) return 0;
		int ans = exbsgs();
		if( ans == -1 ) printf( "No Solution\n" );
		else printf( "%lld\n", ans );
	}
	return 0;
}

miller_rabbin

miller_rabbin是一个质数判断算法,采用随机算法能够在很短的时间里判断一个数是否是质数

首先,根据费马定理,若\(p\)为质数,取\(a<p\),则\(a^{p-1}\%p=1\),但反过来则不一定成立
因为有一类数——卡迈克尔数,它不是质数,但也满足\(a^{p-1}\%p=1\)

如何将卡迈克尔数甄别出来,这里要用到二次探测方法

\(p\)是要判断的数,将\(p-1\)分解为\(2^r*k\),则有\(a^{p-1}=(a^k)^{2^r}\)
可以先求出\(a^k\),然后将其不断平方,并取模\(p\)
如果第一次出现模\(p\)的值为\(1\),并且检测上一次的值是否为\(-1\)
如果是,则\(p\)是质数,反之不是质数

证明:
\(y^2\%p=1\),即\((y-1)*(y+1)-1=k*p\)
如果\(p\)为质数,则\(\begin{cases}y-1=0\\y+1=p\end{cases}\)
所以\(y\)的值只能为\(1\)\(p-1\)
那么如果\(y≠1\)\(y≠p-1\),就可以断定\(p\)不是质数
其次如果\(r\)次平方过程中,模\(p\)的值都不为\(1\),肯定也不是质数

如果\(r\)次平方的过程中,结果为\(1\)了,而且第一次结果为\(1\),而之前的结果不为\(-1\),则判断其不是质数,这样的方法称之为二次探测

通过了二次检测,是不是就一定能够断定\(p\)是质数呢?也不一定
但我们认为它逼近正确,是质数的概率很大
在这里插入图片描述

于是,我们可以多次选择 ,如果选了许多\(a\)来做上述的操作,都不能判断出\(p\)是合数,则\(p\)多半就是一个质数了

miller_rabbin算法是一个随机算法,并不能完全保证结果准确,但是出错的概率非常小
我们取\(a\)的次数越多,出错的概率越小
一般情况下,取20~30次\(a\),正确率接近100%了

如何分析这个概率呢?
对于一个质数,它无疑可以百分百地通过检测
而一个合数,可能有极少量的\(a\)值,能让它通过费马定理和二次检测
假设这些a值的个数,占\(\frac{1}{10}\)
那么通过\(1\)次检测的为\(\frac{1}{10}\),通过两次检测的概率为\(\frac{1}{100}\),通过十次检测的概率就是\(\frac{1}{10^{10}}\)

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <ctime>
using namespace std;
#define int long long

int qkmul( int x, int y, int mod ) {
	int ans = 0;
	while( y ) {
		if( y & 1 ) ans = ( ans + x ) % mod;
		x = ( x + x ) % mod;
		y >>= 1;
	}
	return ans;
}

int qkpow( int x, int y, int mod ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = qkmul( ans, x, mod );
		x = qkmul( x, x, mod );
		y >>= 1;
	}
	return ans;
}

signed main() {
	srand( ( unsigned ) time( NULL ) );
	int n;
	while( cin >> n ) {
		if( n == 2 || n == 3 ) {
			printf( "It is a prime number.\n" );
			continue;
		}
		else if( n == 1 || n & 1 == 0 ) {
				printf( "It is not a prime number.\n" );
				continue;
			}
		int r = 0, k = n - 1;
		while( ! ( k & 1 ) ) r ++, k >>= 1;
		int ans = 1, last = 0;
		bool flag = 0;
		for( int i = 0;i <= 10;i ++ ) {
			ans = rand() % ( n - 1 ) + 1;
			ans = qkpow( ans, k, n );
			if( ans == 1 ) continue;
			last = ans;
			for( int j = 0;j < r;j ++ ) {
				ans = qkmul( ans, ans, n );
				if( ans == 1 ) {
					if( last != n - 1 ) flag = 1;
					break;
				}
				last = ans;
			}
			if( flag ) break;
			if( ans != 1 ) { 
				flag = 1;
				break;
			}
		}
		if( flag ) printf( "It is not a prime number.\n" );
		else printf( "It is a prime number.\n" );
	}
	return 0;
}

but!wait a minute...
刚开始看上面那份代码的原始模样,真的又臭又长,没有看的兴趣,结果后来仔细敲了一遍发现也就那样嘛
果然还是爷码风美丽
在这里插入图片描述

在网上看到的miller_rabbin基本上都没有考虑卡迈克尔数
且此时的概率不再是\((\frac{1}{10})^x\)而是\((\frac{1}{4})^x\),但我也不造为什么,别问我>﹏<

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iostream>
using namespace std;
#define int long long
#define MAXN 50

int qkmul( int x, int y, int mod ) {
	int ans = 0;
	while( y ) {
		if( y & 1 ) ans = ( ans + x ) % mod;
		x = ( x + x ) % mod;
		y >>= 1;
	}
	return ans;
}

int qkpow( int x, int y, int mod ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = qkmul( ans, x, mod );
		x = qkmul( x, x, mod );
		y >>= 1;
	}
	return ans;
}

bool miller_rabbin( int n ) {
	if( n == 2 ) return 1;
	if( n < 2 || ! ( n & 1 ) ) return 0;
	srand( ( unsigned ) time ( NULL ) );
	for( int i = 1;i <= MAXN;i ++ ) {
		int x = rand() % ( n - 2 ) + 1;
		if( qkpow( x, n - 1, n ) != 1 ) return 0;
	}
	return 1;
}

signed main() {
	int n;
	while( cin >> n ) {
		if( miller_rabbin( n ) )
			printf( "It is a prime number.\n" );
		else
			printf( "It is not a prime number.\n" );
	}
	return 0;
}

高斯消元

指路
已经有一篇模板合集了,不重复敲写,反正高斯消元很简单不是吗??
在这里插入图片描述

posted @ 2020-09-23 21:35  读来过倒字名把才鱼咸  阅读(259)  评论(0编辑  收藏  举报