math Algorithm 정리
목차
| 제목 | 날짜 |
|---|---|
| 소수 (prime number) | 24.09.07 |
| 유클리드 알고리즘 (Euclidean Algorithm) | 24.09.07 |
| 빠른 거듭제곱 연산(exponentiation-by-squaring) | 24.09.07 |
| 행렬에서의-연산 (caculation-by-matrix) | 24.09.07 |
| 피보나치 수열 | 24.09.07 |
| 모듈러 연산(Modular Arthimetic) | 24.09.08 |
| 진법변환(base conversion | 24.09.08 |
소수 (prime number)
소수판정 - 밀러라빈 알고리즘
소수를 판별하려면 2부터 $sqrt(n)$까지의 모든 수에 대해서 $n$을 나누어 떨어지는 수가 있는지를 확인 하면 되고, 그런 수가 없다면 해당 수를 소수라고 판별할 수 있다.
하지만 이는 $n$이 작을때 사용하는 것이고, 만약 $n$이 크다면, 밀러-라빈 알고리즘을 사용할 수 있다. 밀러라빈 알고리즘은 확률 알고리즘으로, 소수의 성질을 이용하여 만들어낸 알고리즘이다. 간단히 설명을 하면 홀수인 $n$이 소수인지 알아보면, $n-1$은 짝수가 되고 이는 $n-1 = 2^s d$로 표기를 할수 있다. n이 소수라고 가정해보았을때, 페르마의 소정리를 이용하면 \(a^{n-1} = a^{2^s d} \equiv 1 \: (\text{mod}\: n) \; \exists \, a \in \mathbb{N}, a < n\) 다음 처럼 나타낼 수 있다. 이를 정리하면 $a^{2^r d} = -1 : (\text{mod} : n) \;or\; a^d = 1 : (\text{mod} : n) \;\; some \;r \in [0,s-1]$이 되는데, 다음 두 식을 만족한다면 n은 아마도 소수라고 유추를 할 수 있다.
만약 잘못 판단하는 경우를 줄이고 싶다면 많은 a를 이용하면 된다.
실제 PS 에서 사용할 때에는 int 범위에서 $a = {2,7,61}$ 정도만 확인 하면 되고, long long 범위에서는 $a = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}$로 테스트하는 게 가장 정확하다고 한다.
이때 시간복잡도는 $O(k \log^3 n)$ 이라고 한다.
유클리드 알고리즘 (Euclidean Algorithm)
$ lcm(a,b) = \frac{ab}{gcd(a,b)}$
다음은 유클리드 알고리즘으로 최소공배수 또는 최대공약수 둘중 하나만 알고 있다면 나머지 하나를 알수 있음
다음 사이트에서 효율적 계산 식을 알수 있음
$(a,b) = (b, a \mod b)$ 이 동작하는 이유는 다음 처럼 설명이 가능하다.
$x = \text{gcd} (a,b) $ 일때, $x$는 $a$와 $b$ 모두의 약수이다. 또한 $ a \mod b$의 약수도 된다. 그러므로 수의 크기를 줄여가면서 효율적으로 구현이 가능하다.
Binary Euclid Algorihtm
Binary Euclid Algorithm 이 알고리즘은 현대 컴퓨터는 shift 연산이 더 효율적이기 때문에, 이 shift 연산을 이용하여 조금더 성능이 좋아진 알고리즘이다.
int binary_gcd(int a, int b){
int k;
if(a == 0) return b;
if(b == 0) return a;
for(k = 0; ((a|b)&1)== 0; ++k){
a >>= 1;
b >>= 1;
}
while(!(a&1)) a>>=1;
do{
while(!(b&1)) b>>=1;
if(a > b) swap(a,b);
b -= a;
b >>=1;
} while(b);
return a << k;
}
확장된 유클리드 알고리즘
정수 $a,b$가 주어졌을 때, $ax + by = gcd(a,b)$를 만족하는 $x,y$를 구하는 알고리즘.
이를 구현하기 전 베주 항등식에 대해 알면 이해하기 편하다.
베주 항등식
두 정수의 최대공약수를 원래 두수의 배수의 합으로 나타낼수 있음.
$a,b$가 주어지고, $d$는 $a,b$의 최대공약수일때, 다음 등식을 성립시키는 $m,n$이 존재한다. $ma + nb = d$
이를 통해서 만약 $d \neq k\times \text{gcd}(a,b) \quad k \in \mathbb{\mathbb{N}}$이라면 이는 x,y가 정수해가 아니라는 뜻이다.
다시 확장 유클리드 알고리즘으로 돌아와서…
베주 항등식에 의해 정수해 $x,y$가 존재한다는 것을 알수 있고, 유클리드 알고리즘을 이용하여 $x,y$를 구해낼 수 있다.
확장된 유클리드 알고리즘에서는 정수해 $x,y$를 찾는 것을 중점으로 두기 때문에 반환하는 값도 $x,y$여야 한다.
$\text{gcd}(a,b) = \text{gcd}(b, a \mod b)$
$a \mod b = a-\left \lfloor a/b\right \rfloor \cdot b$
이 두가지 식을 이용하여
$bx’ + (a \mod b)y’ = \text{gcd}(a,b)$ 라고 나타내고, 위의 식을 대입해준다면
$ay’ + b(x’-\left \lfloor a/b\right \rfloor \cdot y’) = gcd(a,b) $ $x = y’ , y = x’-\left \lfloor a/b\right \rfloor \cdot y’$
임을 구할 수 있다.
이를 코드로 구현할때에는 tuple을 이용하여 구할수 있다.
tuple<int,int,int> gcd(int a, int b){
if (b == 0) return {1,0,a};
else {
int x,y, g;
tie(x,y,g) = gcd(b, a%b);
return {y, x-(a/b)*y , g};
}
}
선형 디오판토스 방정식
부정방정식 $ax+by = c$에서 정수해 만을 구하는 방정식을 말한다. 만약 $c = \text{gcd}(a,b)$ 라면 베주 항등식과 동일하다.
디오판토스 방정식에서 정수 해를 구하려면 c가 $\text{gcd}(a,b)$의 배수여야 한다. $c = k \cdot \text{gcd}(a,b)$ 인경우, $ax+by = \text{gcd}(a,b)$를 통해 $x,y$를 구해주고 거기에 $k$를 곱해주면 된다.
또한 디오판토스 방정식의 해는 유일하지 않고, 하나의 해만 알고 있다면 무한히 많은 해를 만들 수 있다.
다음은 그 조합에 해당된다.
\[\left ( x+\frac{pb}{\text{gcd}(a,b)}, y-\frac{pa}{\text{gcd}(a,b)} \right ) , \exists p \in \mathbb{N}\]빠른 거듭제곱 연산 (Exponentiation by squaring)
정수론 문제에서는 $x^n \mod m$을 빠른 시간에 계산해야하는 경우가 생긴다. 이때 다음 성질을 이용하여 거듭제곱을 빠르게 풀수 있다.
int fast_pow(int x, int n, int m){
if(n == 0) return 1 % m;
long long u = fast_pow(x,n>>1, m);
u = (u*u)%m;
if(n&1) u = (u*x)%m;
return u;
}
위의 아이디어를 이용한다면 행렬의 거듭제곱에서 빠르게 사용할 수 있다.
행렬에서의 연산(caculation by matrix)
행렬에서의 빠른 거듭제곱
크기가 $n \times n$인 행렬 $A,B$가 있을 때 $A \times B$ 은 naive 하게 구현하면 $O(n^3)$을 가지게 된다.
strassen’s algorithm
스트라센 알고리즘은 위의 naive 구현보다 조금 더 빠른 $O(n^{2.81})$ 에 작동하는 알고리즘이다. strassen’s algorithm 을 참고할 수 있다.
참고로 현존하는 행렬곱셈 알고리즘 중 가장 빠른 알고리즘은 Le gall의 알고리즘이라고 한다. le gall’s algorithm
다시 거듭제곱으로 돌아온다면..
이번엔 크기가 $n \times n$인 행렬 $A$가 있을 떄 $A^k$를 naive 하게 구현하면 시간이 오래 걸릴 것이다. 이를 보완하기 위해 위의 빠른 거듭제곱의 방식을 이용하여 $O(n^3 \log k)$정도로 줄일 수 있다.
matrix mat_pow(matrix A, int k){
if( k== 1) return A;
if( k &1) return A* mat_pow(A,k-1);
else{
matrix c = mat_pow(A,k/2);
return c*c;
}
}
행렬 제곱의 합
행렬 제곱의 합을 구하는 경우도 있을 텐데, 이도, 빠른 거듭제곱 발상을 이용하여 구할 수 있다. \(S^B=\begin{cases}(A^{B/2}+I)\cdot S^{B/2} &\text{if}~2 \mid B \\(A^{\left\lfloor B/2 \right\rfloor+1}+(A^{\left\lfloor B/2 \right\rfloor+1}+I)\cdot S^{\left\lfloor B/2 \right\rfloor} & \text{if}~2 \nmid B\end{cases}\)
matrix sum_pow_matrix(int k, matrix A){
if (k == 1) return A;
if(k%2 == 1){
return sum_pow_matrix(k/2,A) * (I + power(k/2, A)) + power(k,A);
}
else return sum_pow_matrix(k/2,A) * (I + power(k/2,A));
}
피보나치 수열
피보나치 수열은 다양한 형태로 구현을 할수 있으며, 그 방식에는 recursion, iterative, matrix, dp 등 다양하게 구현이 가능하다. 피보나치 수는 식에서 유도가능한 다양한 성질을 가지고 있으므로, 위의 링크를 통해 잘 살펴보는 것이 좋다.
recursion, iterative에서는 $O(N)$ 의 시간이 걸리게 되고, matrix, dp 에서는 $O(\log N)$ 정도 까지 시간을 줄일 수 있다. 그중 map을 이용한 dp 방식의 코드하나만 첨부한다.
#define MOD 1000000007
unordered_map<int64_t, int64_t> Fb;
int64_t Fibo(int64_t N) {
if (N == 0) return 0;
if (N <= 2) return 1;
if (Fb.find(N) != Fb.end()) return Fb[N];
int64_t k;
if (N & 1) {
k = (N + 1) / 2;
return Fb[N] = (Fibo(k) * Fibo(k) + Fibo(k - 1) * Fibo(k - 1)) % MOD;
}
else {
k = N / 2;
return Fb[N] = (2 * Fibo(k - 1) + Fibo(k)) * Fibo(k) % MOD;
}
}
제켄도르프의 정리
모든 양의 정수를 한개 이상의 서로 다른 피보나치 수의 합으로, 그 합에 연속된 두 피보나치 수가 포함되지 않도록 유일히 나타낼 수 있다.
모듈러 연산 (Modular arithmetic)
일반적인 연산 (+, -, $\times$) 에서는 연산할때마다 모듈러 연산을 취해주면 된다.
이때 나눗셈인 경우 조금 달라지게 되는데, $(a/b) \mod \;m$인 경우에 $\text{gcd}(b,m) = 1$ 즉 $b,m$이 서로소라면 모듈러 역원 $\text{inv} (b)$를 구한후에 이를 곱해주는 방식으로 연산하게 된다.
오일러 정리
오일러 피 함수(euler’s totient function) $ \varphi (n)$은 1이상 $n$ 이하의 정수중 $n$과 서로소인 정수의 개수를 나타내는 함수이다. 이 값은 $n$을 소인수 분해한 후에 얻어낼 수 있다.
이때 오일러 정리란 모든 서로소인 정수 $x,m$에 대해 다음 식이 성립한다. \(x^{\varphi(m)}\, \text{mod} \,m = 1\) 만약 이식에서 m이 소수인 경우에는 $\varphi(m) = m-1$과 같으므로 \(x^{m-1}\, \text{mod} \,m = 1\) 로 나타낼 수 있고, 이는 페르마 소정리(Fermat’s little theorem)이라고 불린다.
모듈러 역원 (Modular multiplicative inverse)
$a,m$이 서로소일때, $a x \equiv 1\;(\text{mod} \,m)$인 경우 $x$를 모듈러 역원이라고 한다.
이를 통해 모듈러 역원이 존재하는 경우는 $x,m$이 서로소인 경우와 동치가 된다. 그리고 오일러 정리를 이용하면 $\text{inv}(x) = x^{\varphi (m)-1}$과 같다라는 것을 알수 있다. $m$이 소수인 경우에는 $\text{inv}(x) = x^{m-2}$가 된다.
수 $x$의 역원을 구하려면 기본적으로 공식에 의해 오일러 피함수를 구하면 된다. 하지만 이를 구하는 것은 소인수분해를 하고, 그 개수를 세야하므로, 시간이 오래 걸리게 된다. 그러므로 Extened Euclidean Algorithm 을 이용하여 구하게 된다. $ax + my \equiv 1$ 을 이용하면 된다.
Multiple inverses
n개의 정수로 이루어진 수열 $A$의 원소를 $a_i$라고 할때, $a_1 … a_n$까지의 역원을 구하는 방법이다. 기본적으로 역원을 구하는 것은 $O(\log(m))$과 동일하지만 n개의 수의 역원을 구하려면 $O(n\log(m))$이 걸리게 된다. 이를 $O(n+log(m))$안에 처리시키는 알고리즘이다. 하지만 이는 n이 상당히 클때 이득이 있다.
다중역원의 방식을 이용한다면, Factorial 수들에 대한 역원도 전처리로 미리 구해낼 수 있다. 이는 이항계수에서 사용되므로, 한번 계산해보면 좋다.
vector<int> m_inv(vector<int> A, int m){
int n = A.size();
vector<int> B(n);
B[0] = A[0];
for(int i = 1; i<n;i++){
B[i] = A[i]*B[i-1];
}
int inv,y,g;
tie(inv,y,g)= gcd(B[n-1],m);
vector<int> ans(n);
ans[n-1] = inv;
for(int i = n-2; i>=1;i--){
ans[i] = A[i+1]*ans[i+1];
}
tie(inv,y,g) = gcd(A[0],m);
ans[0] = inv;
return ans;
}
다음처럼 식을 구성해볼수 있다.
중국인의 나머지 정리 (Chinese Remainder Theorem)
$x=a_i \, \text{mod}\;m_i \;,\; \forall i\in [1,n] $ 인 문제가 나온다면 사용할 수 있는데, 이때 $m_i$의 모든 조합은 서로소이다. 중국인의 나머지 정리 증명에서 자세하게 증명을 할 수 있다.
일단 $n_i$가 모두 서로소여야 하는 것이 전제이다. N을 $n_i$의 곱이라고 생간한다. 그리고 $N_i = \frac{N}{n_i}$ 라고 생각한다면 $N_i, n_i$는 서로소가 된다. 이는 $N$에서 $n_i$만 없어지는 것이기 때문에 무조건 서로소가 된다. 그럼 베주 항등식에 의해 $M_i N_i + m_i n_i = 1$인 정수해 $M_i, m_i$를 찾을 수 있다. 그럼 우리가 구해야하는 값 x는 $\sum a_iM_iN_i$가 된다. 이 식을 다음 합동식으로 바꿀 수 있다. 이 때 $M_iN_i \equiv 1 \mod n_i$ 이 식이 나오게 된다. 즉 $M_i$는 $N_i$의 역원인것을 알수 있고, 확장 유클리드 알고리즘을 이용하여 역원을 구할 수 있다. 그리고 위의 x를 구하고, N으로 모듈러연산을 취해주면 된다.
int mod(int x, int n){
x %= n;
if(x<0) x+= n;
return x;
}
tuple<int,int,int> gcd(int a, int b){
if (b == 0) return {1,0,a};
else {
int x,y, g;
tie(x,y,g) = gcd(b, a%b);
return {y, x-(a/b)*y , g};
}
}
int crt(vector <int> a, vector<int> n){
int NN = 1;
vector <int> N(n.size());
for(auto i : n) NN *= i;
for(int i = 0;i<N.size();i++){
N[i] = NN/n[i];
}
vector<int> M(N.size());
for(int i = 0; i<M.size();i++){
int y, g;
tie(M[i],y,g) = gcd(N[i],n[i]);
}
int x = 0;
for(int i = 0;i<M.size();i++){
x += a[i]*N[i]*M[i];
}
x = mod(x,NN);
return x;
}
서로소가 아니라면?
서로소가 아닌 경우도 존재하는데, 이때는 이 링크에서 나온 공식을 토대로 연쇄적으로 $x$를 풀어나가면 된다.
int mod(int x, int n){
return (x+n)%n;
}
int _gcd(int a, int b){
if(b == 0) return a;
return gcd(b, a%b);
}
tuple<int,int,int> gcd(int a, int b){
if (b == 0) return {1,0,a};
else {
int x,y, g;
tie(x,y,g) = gcd(b, a%b);
return {y, x-(a/b)*y , g};
}
}
int crt(vector <int> a, vector<int> m){
int n = m.size();
int a1 = a[0], m1 = m[0];
a1 %= m1;
for(int i = 1; i<n;i++){
int a2 = a[i], m2 = m[i];
int g = _gcd(m1,m2);
if(a1%g != a2%g) return -1;
auto [x,y,_] = gcd(m1/g,m2/g);
int M = m1 / g * m2;
a1 = (a1*(m2/g)%M) * y % M + (a2*(m1/g)%M)*x % M;
a1 = mod(a1, M);
m1 = M;
}
return a1;
}
꼭 문제를 풀때, 조건이 서로소인지, 아닌지를 생각하고 두개의 코드를 잘 이용해야한다.
진법변환 (base conversion)
BOJ 1112 이 문제를 보면 일반적인 진법변환이 아닌, 음수 진법변환도 포함되어 있는 경우가 있다. 이런 경우에는 어떻게 해야할지 수학적으로 접근을 해보아야한다.
\[x = a_nb^n+ a_{n-1}b^{n-1}+...+a_1b+a_0\]우리는 어떤 10진수 $x$를 위에 처럼 표기할 수 있다. 일단 $x,b$가 양수일때 생각해보자.
모듈러 연산을 이용하면 다음을 알아낼 수 있다.
\[m_0 = x \equiv a_0 \mod b\]그리고 위의 식을 다음처럼 변형할 수 있다.
\[\frac{m_0-a_0}{b} = a_nb^{n-1}+ a_{n-1}b^{n-2}+...+a_1\]이 식은 다음 합동식이 된다. \(m_1 = \frac{m_0-a_0}{b} \equiv a_1 \mod b\)
즉 일반화를 시켜보면 다음처럼 만들 수 있다.
\[m_i = \frac{m_{i-1}-a_{i-1}}{b} \equiv a_i \mod b\]이를 이용하면 $a_i$를 모두 구할 수 있고, $a$를 뒤집으면 진법이 변환이 된다.
이제 $x$가 음수고 $b$가 양수인 경우 생각해보면, 위에서 그대로 구하고, 앞에 마이너스 부호만 붙여주면 된다.
그럼 b가 음수인 경우를 생각해보자. b가 음수인경우 모듈러가 음수인 경우는 불가능 하기 때문에, $|b|$로 모듈러를 취해주면 된다. 나머지 과정은 위의 증명을 통해 똑같다는 것을 알 수 있다.
pair<bool,vector<int>> crt(int n, int b){
vector <int> ans;
bool p = b>0 && n<0;
n = b>0 ? abs(n):n;
while(n!=0){
ans.push_back( mod(n,abs(b)));
n =((n-ans.back())/b);
}
if(ans.empty()) ans.push_back(0);
reverse(ans.begin(), ans.end());
return {p, ans};
}