算法之数论应用篇(二)
算法之數論應用篇二
- 最大公約數
- 線性篩
- Hankson的趣味題
- 歐拉函數
- 前言
- 可見的點(數學知識+歐拉函數)
- 最大公約數(可見的點擴展)
- 同余
- 取模的性質
- 定義
- 基本性質
- 運算規則
- 重要定理
- 重要定理
- 費馬小定理
- 歐拉定理
- 擴展歐幾里得算法
- 乘法逆元
- 線性同余方程式
- 同余方程(數學知識+線性同余方程式)
- 青蛙的約會(同余+數學方程式)
- 曹沖養豬(中國剩余定理)
最大公約數
題目轉送門
線性篩
void get_primes(int n ){for(int i = 2 ; i <= n ; ++i){if(!st[i])primes[cnt ++] = i;for(int j = 0 ; primes[j] <= n / i; ++j){st[i*primes[j]] = true;if(i % primes[j] == 0)break;}} }利用了每個合數必有一個最小素因子。每個合數僅被它的最小素因子篩去正好一次。所以為線性時間。
代碼中體現在:
if(i%prime[j]==0)break;
prime數組 中的素數是遞增的,當 i 能整除 prime[j],那么 i*prime[j+1] 這個合數肯定被 prime[j] 乘以某個數篩掉。
因為i中含有prime[j], prime[j] 比 prime[j+1] 小。接下去的素數同理。所以不用篩下去了。
在滿足i%prme[j]==0這個條件之前以及第一次滿足改條件時,pr[j]必定是pr[j]*i的最小因子.
Hankson的趣味題
題目轉送門
解題思路
這題是求滿足條件的x分個數。條件是 gcd(a,x) = b , lcm(x,c) = d;
由lcm(c,x) = d ,我們知道x 一定為d的約數。那么我們就可以從d的約數中找符合條件的x。不過因為在找約數的時間復雜度是d\sqrtze8trgl8bvbqd?的,我們需要進行一些優化。我們知道一個數的余數可以有其分解質因數的質因子組合而成。因此我們就可以將d先進行分解質因數,之后通過dfs進行組合形成約數。因為2?109內的數含有的不同質因子的個數不會超過10個(上一篇提到),因此速度會得到很大的優化。
代碼:
#include<iostream> #include<algorithm> #include<cmath> #include<cstring> using namespace std; const int N = 2000 , M = 50010 ; typedef long long LL; int n ; int primes[M] , cnt; bool st[M];struct Factor{int p , s ; }factor[M]; int fcnt ,dcnt; int divtor[M];void get_primes(int n ){for(int i = 2 ; i <= n ; ++i){if(!st[i])primes[cnt ++] = i;for(int j = 0 ; primes[j] <= n / i; ++j){st[i*primes[j]] = true;if(i % primes[j] == 0)break;}} }int gcd(int a , int b){if(b == 0 )return a;else return gcd(b , a % b); }void divite(int n ){fcnt = 0;for(int i = 0 ; i < cnt && primes[i] <= n / primes[i] ; ++i){int p = primes[i];if(n % p == 0){int s = 0;while(n % p == 0) n /= p , s ++;factor[++fcnt] = {p,s};}}if(n > 1)factor[++fcnt] = { n ,1}; }void dfs(int u , int p){if( u > fcnt){divtor[dcnt++] = p;return ;}for(int i = 0 ; i <= factor[u].s ; ++i){dfs(u+1 , p);p *= factor[u].p;} }int main(){get_primes(M - 1);cin>>n;while(n--){int a ,b , c ,d;cin>>a>>b>>c>>d;divite(d);dcnt = 0;dfs(1,1);int res = 0;for(int i = 0 ; i < dcnt ; ++i){int x = divtor[i];if(gcd(a,x) == b && (LL)c * x / gcd(x,c) == d) res++;}cout<<res<<endl;}return 0; }歐拉函數
前言
歐拉函數是一個函數嘛?不它不是的!對于一個數N的歐拉函數我們常記作ψ(N)\psi(N)ψ(N),它的含義是從1~N中與N互質的個數。
其中
ψ(N)=N?Π質數p∣Np?1p\psi(N) = N * \Pi_{質數p|N}\frac{p-1}{p}ψ(N)=N?Π質數p∣N?pp?1?
可見的點(數學知識+歐拉函數)
題目轉送門
解題思路
分析題目容易發現,除了(1,0),(0,1)和(1,1)這三個釘子外,一個釘子(x,y)能被看到,當且僅當1≤x,y≤N.x+y并且ged(x,y)=1。
在1≤x,y≤N,x≠y1≤x,y≤N,x \ne y1≤x,y≤N,x?=y中能看到的釘子關于過(0,0)和(N,N)的直線對稱。我們可以考慮其中的一半,即1≤x<y≤N。換言之,對于每個2≤y≤N,我們需要統計有多少個x滿足1≤x<y并且gcd(x,y)=1。這樣的x的數量恰好就是(y)。
綜上所述,本題的答案就是 3+2?Σi=2Nψ(i)。3+2* \Sigma_{i=2}^N\psi(i)。3+2?Σi=2N?ψ(i)。我們要做的事情就是求出2~N中每個數的歐拉函數.
代碼:
#include<iostream> #include<algorithm> using namespace std; const int N = 1010; int n ; int primes[N]; int phi[N];void init(int n){for(int i = 1 ; i <= n ; ++i)phi[i] = i;for(int i = 2 ; i <= n; ++i){if(phi[i] == i)for(int j = 1 ; j <= n / i ; ++j)phi[j*i] = phi[j*i]/i * (i-1);} }int main(){init(N-1);int T ;cin>>T;for(int i = 1; i <= T ; ++i){cin>>n;int res = 0;for(int i = 2 ; i <= n ; ++i) res += phi[i];cout<<i<<" "<<n<<" "<< res * 2 + 3 <<endl;}return 0; }最大公約數(可見的點擴展)
題目轉送門
解題思路
之所以說是上一題的擴展是有原因的,對于這個題目是求 1 - N ,中 滿足 gcd(x,y)=p,p是質數gcd(x,y) = p ,p是質數gcd(x,y)=p,p是質數的點對數。兩個數x,y之間的最大公約數為p,那么其除以p之后也就是gcd(x/p,y/p)=1gcd(x/p,y/p) = 1gcd(x/p,y/p)=1因為p是它們的最大公約數。而知就轉換成了上一題。有一些出入的就是這里是從 1 開始計數的,同時對于每一個p我們加上 1 ~ N/p 范圍內對于上一題的解。
解題步驟
代碼:
#include<iostream> #include<algorithm> #include<cstring> using namespace std; typedef long long LL; const int N = 1E7 + 10; int n ; int primes[N] ,phi[N]; int st[N] , cnt; LL sum[N];void init(int n ){for(int i = 2 ; i <= n ; ++i){if(st[i] == 0){primes[cnt++] = i;phi[i] = i-1;}for(int j = 0; primes[j] <= n / i ; ++j){st[primes[j] * i] = 1;if( i % primes[j] == 0){phi[primes[j] * i] = phi[i] * primes[j];break;}phi[primes[j] * i] = phi[i] * ( primes[j] - 1) ;}}for(int i = 1 ; i <= n ; ++i)sum[i] = sum[i-1] + phi[i]; }int main(){cin>>n;init(n);LL ans = 0;for(int i = 0 ; i < cnt ; ++i)ans += sum[n/primes[i]]*2 + 1 ;cout<<ans <<endl;return 0; }同余
取模的性質
在介紹有關同余的一系列性質之前有必要列舉一下取模的一些性質
定義
給定一個正整數p,任意一個整數n,一定存在等式 :
n = kp + r ;
其中 k、r 是整數,且 0 ≤ r < p,則稱 k 為 n 除以 p 的商,r 為 n 除以 p 的余數。
對于正整數 p 和整數 a,b,定義如下運算:取模運算:a % p(或a mod p),表示a除以p的余數。
說明:
基本性質
自反性:(a % p)=(a % p)意味a≡a (% p)
對稱性:a≡b (% p)等價于b≡a (% p)
傳遞性:若a≡b (% p)且b≡c (% p) ,則a≡c (% p)
運算規則
模運算與基本四則運算有些相似,但是除法例外。其規則如下:
因此對于除法取模,我們就會利用乘法逆元來解決。下面我們會提到如何求乘法逆元。
重要定理
(a * c) ≡ (b * d) (%p),(a / c) ≡ (b / d) (%p); 我們證明加法,其余同理。由前提我們令
a=k1?p+x,b=k2?p+x,c=k3?p+y,d=k4?p+ya = k_1*p + x , b = k_2*p+x , c = k_3*p+y , d = k_4*p+ya=k1??p+x,b=k2??p+x,c=k3??p+y,d=k4??p+y那么(a+c)=(k1+k3)?p+x+y(b+d)=(k2+k4)?p+x+y(a+c) = (k_1+k_3)*p+x+y \quad (b+d) =(k_2+k_4)*p+x+y(a+c)=(k1?+k3?)?p+x+y(b+d)=(k2?+k4?)?p+x+y它們對p取模的結果都是(x+y)(x+y)%p(x+y)
重要定理
費馬小定理
費馬小定理:a(p?1)≡1(modp)a^(p-1) ≡ 1(mod p)a(p?1)≡1(modp)
前提:p為質數,且a,p互質
在證明之前我們先證明兩個引理
∵a,p互質
∴a,2a,3a,4a…(p-1)a也是mod p的完全剩余系
∴123…(p-1)a ≡ a2a3a…(p-1)*a (mod p)
∴ (p-1)! ≡ (p-1)!*a^(p-1) (mod p)
兩邊同時約去(p-1)!
a^(p-1) ≡ 1 (mod p)
歐拉定理
對于費馬定理,我們知道 p(質數)的歐拉函數就是 p - 1
因此可以成
若正整數a,n互質,則aψ(n)≡1(mod)n若正整數a,n互質,則a^{\psi(n)}≡1 (mod)n若正整數a,n互質,則aψ(n)≡1(mod)n
推論
若正整數a,n互質,則對于任意正整數b,有ab≡abmodψ(n)(mod)n若正整數a,n互質,則對于任意正整數b,有a^≡a^{b mod \psi(n)}(mod)n若正整數a,n互質,則對于任意正整數b,有ab≡abmodψ(n)(mod)n
而這個就是快速冪取??尚械囊粋€依據
擴展歐幾里得算法
定理
對于任意整數a,b,存在一對整數x和y,滿足ax+by=gcd(a,b)對于任意整數a,b,存在一對整數x和y,滿足ax + by = gcd(a,b)對于任意整數a,b,存在一對整數x和y,滿足ax+by=gcd(a,b)
證明
在歐幾里得算法的最后一步,即b=0時,顯然有一對整數x=1,y=0,使得a?1+0?0=gcd(a,0)。a*1+0*0=gcd(a,0)。a?1+0?0=gcd(a,0)。
若b>0,則gcd(a,b)=gcd(b,amodb)。gcd(a,b)=gcd(b,a mod b)。gcd(a,b)=gcd(b,amodb)。假設存在一對整數x,y,滿足b?x+(amodb)?y=gcd(b,amodb)b*x+(amodb)*y=gcd(b,a mod b)b?x+(amodb)?y=gcd(b,amodb),因為bx+(amodb)y=bx+(a??a/b?b])?y=ay+b(x?[a/b]y)bx+(amodb)y=bx+(a-\lfloor a/b\rfloor b])*y=ay+b(x-[a/b]y)bx+(amodb)y=bx+(a??a/b?b])?y=ay+b(x?[a/b]y),所以令 x′=y,y′=x?∣a/b∣y,x′=y,y′=x?∣a/b∣y,x′=y,y′=x?∣a/b∣y, 就得到了 。ax′+by′=gcd(a,b)。 對歐幾里得算法的遞歸過程應用數學歸納法,可知Bézout定理成立。
證畢。
代碼:
int exgcd(int a , int b ,int & x ,int & y){if(b == 0){x = 1 , y = 0;return a;}int d = exgcd(b , a % b , y ,x);y -= a / b * x;對于更為一般的方程ax+by=c,它有解當且僅當d|c。(d = gcd(a,b))我們可以先求出ax+by=d的一組特解xo,yo,然后令xo,yo同時乘上c/d,就得到了ax+by=c的一組特解(c/d)?xo,(c/d)?yo。ax+by=d的一組特解x_o,y_o,然后令x_o,y_o同時乘上c/d,就得到了ax+by=c的一組特解(c/d)* x_o,(c/d)* y_o。ax+by=d的一組特解xo?,yo?,然后令xo?,yo?同時乘上c/d,就得到了ax+by=c的一組特解(c/d)?xo?,(c/d)?yo?。
事實上,方程ax+by=c的通解可以表示為:
x=cd?x0+(kbd)y=cd?y0?(k?ad)(k∈Z)x=\frac {c}ze8trgl8bvbq*x_0+(k\frac ze8trgl8bvbq) \quad y=\frac {c}ze8trgl8bvbq*y_0 -(k*\frac {a}ze8trgl8bvbq) (k∈Z)x=dc??x0?+(kdb?)y=dc??y0??(k?da?)(k∈Z)
其中k取遍整數集合.
乘法逆元
定理
若b,p互質,并且b∣a,則存在一個整數x,使得ab≡a?x(modp)若b,p互質,并且b|a,則存在一個整數x,使得\frac{a} ≡ a*x (mod p)若b,p互質,并且b∣a,則存在一個整數x,使得ba?≡a?x(modp),稱x為b的模p的乘法逆元,記作b?1(modp)稱x為b的模p的乘法逆元,記作b^{-1}(mod p)稱x為b的模p的乘法逆元,記作b?1(modp)
因為ab≡a?b?1≡ab?b?b?1\frac{a} ≡ a*b^{-1}≡\frac{a} *b*b^{-1}ba?≡a?b?1≡ba??b?b?1,所以b?b?1≡1(modp)b*b^{-1}≡1(mod p)b?b?1≡1(modp)如果m是質數,且b < p ,那么由費馬小定理有 ,bp?2b^{p-2}bp?2是b的乘法逆元。
如果只是b,p互質 , 我們也可以通過求同余方程式 b?x≡1(modp)b*x≡1(mod p)b?x≡1(modp)
現在有了乘法逆元對于除法的取模運算我們也可以適當的進行加速運算了。
線性同余方程式
定義:給定整數a,b,m,求一個整數x滿足a?x≡b(modm)a*x ≡b(modm)a?x≡b(modm),或者給出無解。因為未知數的指數為1,所以我們稱之為一次同余方程,也稱線性同余方程。
a?x≡b(modm)a*x ≡b(modm)a?x≡b(modm)等價于a?x?b是ma*x-b是ma?x?b是m的倍數,不妨設為?y-y?y倍。于是,該方程可以改寫為a?x+m?y=ba*x+m*y=ba?x+m?y=b。
根據Bézout定理及其證明過程,線性同余方程有解當且僅當gcd(a,m)∣b線性同余方程有解當且僅當gcd(a,m)|b線性同余方程有解當且僅當gcd(a,m)∣b。而這也證明了上面情況乘法逆元的存在性。b?x≡1(modp)b*x≡1(mod p)b?x≡1(modp)??梢詫懗?span id="ze8trgl8bvbq" class="katex--inline">bx+py=1bx + py = 1bx+py=1因為b,p互質,所以 gcd(b,p)=1gcd(b,p) = 1gcd(b,p)=1那么解一定存在。
在有解時,先用歐幾里得算法求出一組整數x0,y0,滿足a?x0+m?y0=gdd(a,m)x_0,y_0,滿足 a?x_0+m?y_0= gdd(a,m)x0?,y0?,滿足a?x0?+m?y0?=gdd(a,m)。然后x=x0?b/gcd(a,m)x=x_0*b/gcd(a,m)x=x0??b/gcd(a,m)就是原線性同余方程的一個解。
方程的通解則是所有模m/gcd(a,m)與x同余的整數。
同余方程(數學知識+線性同余方程式)
題目轉送門
解題思路
題目要求的是滿足條件的最小正整數。而通解的形式是x0+kbx_0 + kbx0?+kb表示的是所有模b與x0x_0x0?的余數相同的最小正整數。 循環來求就是
不過我們發現其實我們是不用循環的,這個(x % b +b ) % b)就是所求解。
代碼:
#include<iostream> #include<algorithm> #include<cstring> using namespace std;int a , b ,x ,y;int exgcd(int a , int b ,int & x ,int & y){if(b == 0){x = 1 , y = 0;return a;}int d = exgcd(b , a % b , y ,x);y -= a / b * x;return d; }int main(){cin>>a>>b;exgcd(a,b,x,y);cout<< (x % b + b)%b<<endl;return 0; }青蛙的約會(同余+數學方程式)
題目轉送門
解題思路
那么它們相遇的時候有方程式
a+mx≡b+nx(modL)a + mx \equiv b+nx \quad(mod L)a+mx≡b+nx(modL)等價于(a?b)+(m?n)x是L的倍數(a-b)+(m-n)x是L的倍數(a?b)+(m?n)x是L的倍數我們假設倍數是-y倍,那么有(m?n)x+Ly=b?a(m-n)x + Ly = b - a(m?n)x+Ly=b?a這就是一次同余方程。那么有解當且僅當 gcd(L,m?n)∣(a?b)gcd(L,m-n) | (a-b)gcd(L,m?n)∣(a?b)。如果有解我么用擴展歐幾里得算法解出x0x_0x0?后x=x0?(b/gcd(L,m?n))x = x_0 * (b/gcd(L,m-n))x=x0??(b/gcd(L,m?n))就是方程的一個特解。通解就是X=x+Ld,d=gcd(L,m?n)X = x + \frac{L}ze8trgl8bvbq ,d = gcd(L,m-n)X=x+dL?,d=gcd(L,m?n),至于最小的正整數解就類似上一題所說為(xmodLd+Ld)modLd(x \quad mod \quad \frac{L}ze8trgl8bvbq + \frac{L}ze8trgl8bvbq)mod \frac{L}ze8trgl8bvbq(xmoddL?+dL?)moddL?
代碼:
#include<iostream> #include<algorithm> #include<cmath> using namespace std;typedef long long LL; LL n , m , a ,b ,L , x , y;LL exgcd(LL a , LL b ,LL & x , LL & y){if(b == 0 ){x = 1 , y = 0;return a;}LL d = exgcd(b , a%b , y , x);y -= a/b * x;return d; }int main(){cin>>a>>b>>m>>n>>L;LL d = exgcd(m - n , L ,x ,y);if((b - a) % d)cout<<"Impossible\n";else{x *= (b - a) / d ;LL t = abs(L/d);cout<< ( x % t + t) % t <<endl;}return 0; }曹沖養豬(中國剩余定理)
題目轉送門
解題思路
這題是一個板子題,對于中國剩余定理是
中國剩余定理是構造了一組解,使得滿足條件。
代碼:
#include<iostream> #include<algorithm> using namespace std; typedef long long LL; const int N = 11; int n ; int a[N] , b[N];int exgcd(LL a , LL b ,LL & x , LL & y){if(!b){x = 1 , y = 0;return a;}LL d = exgcd(b , a % b , y ,x);y -= a/b * x;return d; }int main(){cin>>n;LL m = 1;for(int i = 1 ; i <= n ; ++i){cin>>a[i]>>b[i];m *= a[i];}LL ans = 0;for(int i = 1 ; i <= n ; ++i){LL x , y ;LL M = m / a[i];exgcd( M, a[i] , x , y);ans += b[i] * M * x; }cout<<(ans % m+ m ) % m<<endl;return 0; }總結
以上是生活随笔為你收集整理的算法之数论应用篇(二)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 算法之基础数论应用篇(一)
- 下一篇: 算法之组合数学及其算法篇(一) ----