原根，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,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\)裡尋找是否出現過相同的值
如果有，代表著同餘，已經找到了答案，如果沒有則是無解
然而。。。
可以稍微轉變一下演算法過程，避免求逆元的操作
則$$a^{k*m-r}\equiv b\ \ (\mod p)$$
兩邊同時乘以\(a^r\)
先求出右邊所有的\(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\)喲(＾Ｕ＾)ノ~ＹＯ
#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;
}
高斯消元
指路
已經有一篇模板合集了，不重複敲寫，反正高斯消元很簡單不是嗎？？