확장 유클리드 알고리즘, 페르마 소정리, 중국인 나머지 정리

확장 유클리드 정리는 $ax+by=c$ ($a,b,c$는 상수 $ x,y$는 정수인 미지수)를 만족시키는 $x,y$를 찾는 알고리즘이다. 예를 들어서 $a=3,b=4,c=10$이라면 $x=2,y=1$을 찾을 수 있게 해주는 것이다. 여기서 반드시 충족되야할 조건은, $a$와 $b$의 최대공약수가 $c$의 약수여야 한다는 것이다. 수식으로 표현하자면 $gcd(a,b)=d$일 때, $d | c$를 만족해야 한다는 것이다. 위 식에선 d=1이고 c=10이므로 충족한다. 이를 충족시키지 못하면 $x,y$는 존재하지 않는다. 

그러면 본격적으로 $ax+by=c$를 만족하는 $x,y$는 어떻게 구할까? 바로 기존의 유클리드 알고리즘을 거꾸로 거슬러 올라가기만 하면 된다. 유클리드 알고리즘은 $gcd(a,b)=gcd(b,a\%b)$이다. 예를 들어서 a=20,b=12를 가정해보자. 유클리드 알고리즘에 따라서

 

$gcd(20,12)=gcd(12,20\%12=8)$

$gcd(12,8)=gcd(8,12\%8=4)$

$gcd(8,4)=4$이다.

 

이를 거꾸로 올라가면

 4=12-8

 4=12-(20-1*12)=2*12-20이 나온다.

확장 유클리드 알고리즘을 따라서 나오는 결과는 $a$와 $b$의 최대공약수이다.  즉 $20x+12y=4$를 만족하는 $x,y$는 $x=2,y=-1$이라는 말이다. 하지만 이는 유일한 해가 아니며 실제로는 해가 무수히 많다. 이 해의 일반항 방정식은

$x=x_0+(b/d)*t, y=y_0-(a/d)*t$ ($x_0,y_0$는 방정식을 만족하는 임의의 정수)이다.   ($x=2+3t, y=-1-5t$)

따라서 상수 $c$값을 맞춰주기 위해서는 $x_0,y_0$에 적절한 정수배를 해주면 되는데, 이를 디오판테인 방정식이라고 한다. 예를 들어서 위와 같이 $a=20,b=12$이고 $c=20$이라면 최대공약수 4에 5배를 해줘야 하므로 $x_0=10,y_0=-5$가 된다. 따라서 일반항 방정식은 

$x=10+3t,y=-5-5t$ 이다.

 

이를 코드로 구현하면 다음과 같다. (이 블로그를 참고함)

 

tuple<int,int,int> ex_euclid(int a,int b){
    if(b==0) return {a,1,0};
    int g,x,y;
    tie(g,x,y)=ex_euclid(b,a%b);
    return {g,y,x-(a/b)*y};
}

 

함수의 인자로 a,b를 넣어주면 반환값으로 첫번째는 a,b의 최대공약수(d), 두번째와 세번째는 $ax+by=d$를 만족시키는 a,b를 반환해준다.

 

지금까지 배운 내용을 토대로 문제를 풀어보자.

 

[BOJ 14565] 역원(Inverse) 구하기

 

 

이 문제를 풀기 위해서는 역원의 개념에 대해서 숙지하고 있어야 한다. 역원이란 연산했을 때, 항등원이 나오도록 하는 수이다. 항등원은 연산해서 자기 자신이 나오는 수이다. 예를 들어 덧셈에 대한 항등원은 '0'이다. $a+?=a$라는 식이 있을 때 ?는 0이기 때문이다. 그리고 덧셈에 대한 항등원은 $a+??=0$을 만족해야 하기 때문에 ??는 $-a$이다. 그럼 비슷하게 곱셈에 대한 항등원은 1, 역원은 $1/a$로 정의할 수 있을 것이다. 그런데 이는 모듈러 연산이 없을 때의 얘기이고, 이 문제는 모듈러 연산이 존재하기 때문에 곱셈에 대한 항등원,역원이 성립하지 않고 새롭게 정의해야 한다. 

 

위 문제에서 물어보는 것은 $Ax≡1 (mod\;n)$에 대해서 물어보고 있다. 이를 유클리드 알고리즘으로 풀어서 쓰면 

$Ax+ny=1$ 형태로 바꿔서 쓸 수 있다. 따라서 위에서 적었던 ex_euclid()함수를 이용해서 반환되는 값(두번째,세번째)가 답이 된다. 여기서 확장 유클리드 정리를 충족시키기 위해서는 언급했다시피 A와 n의 최대공약수가 1의 약수여야 한다는 것이다. 1의 약수는 1 하나밖에 없으므로 A와 n이 서로소이어야만 성립한다는 뜻이다. 서로소가 아니면 역원은 존재하지 않는다. 

그리고 또 주의해야할 점은 ex_euclid()함수를 돌고 나온 값이 음수일 수도 있다는 점이다. 이 문제에서는 이럴 때는 디오판테인 방정식에서처럼 적당히 $b/d$값을 더해줘서 양수로 맞춰주면 된다.

tuple<ll,ll,ll> ex_euclid(ll a,ll b){
    if(b==0) return {a,1,0};
    ll g,x,y;
    tie(g,x,y)=ex_euclid(b,a%b);
    return {g,y,x-(a/b)*y};
}

int main(){
    ll a,b;   
    scanf("%lld %lld",&a,&b);

    printf("%lld ",a-b);

    ll GCD,x,y;
    tie(GCD,x,y)=ex_euclid(a,b);

    if(GCD!=1) {printf("-1\n"); return 0;}
    if(y<0){
        y+=a;
    }
    printf("%lld\n",y);

    return 0;
}

 

 [BOJ 3955] 캔디 분배

 

문제 상황을 요약해보면 

$$CX=KY+1$$ ($X$=총 봉지 수, $Y$=한 사람당 사탕 수)

$CX-KY=1$ 꼴의 식으로 표현된다. 이 식에 확장 유클리드 알고리즘을 적용하면 된다.

당연히 위에서 언급했다시피 $C$와 $K$의 최대공약수가 결과로 나오기 때문에 1이 나오기 위해서는 $C$와 $K$가 서로소여야만 답이 존재한다.

그리고 주의해야할 점은 이 사탕은 최대 1봉지 이상 10억봉지 이하만 살 수 있끼 때문에 $0<X<10^9$이고 한 사람당 최소 하나 이상의 사탕은 받아야하므로 $Y>0$인 조건을 만족하는 $X,Y$값을 찾아야한다. 

디오판테인 방정식의 형태로 $X=x0+kt, Y=y0-ct$를 구하고 위의 부등식에 대입해서 적절한 t값을 구해야 한다. 분수의 부등식에 따른 구현은 myceil, myfloor함수를 사용했다.

tuple<ll,ll,ll> ex_euclid(ll a,ll b){
    if(b==0) return {a,1,0};
    ll g,x,y;
    tie(g,x,y)=ex_euclid(b,a%b);
    return {g,y,x-(a/b)*y};
}

ll myceil(ll a,ll b){
    if(b<0){a=-a; b=-b;}
    if(a>=0) return (a+b-1)/b;
    return a/b;
}

ll myfloor(ll a,ll b){
    if(b<0){a=-a; b=-b;}
    if(a>=0) return a/b;
    return (a-b+1)/b; 
}

int main(){
    int t;  scanf("%d",&t);

    while(t--){
        ll k,c; scanf("%lld %lld",&k,&c);
        ll G,S,T;
        tie(G,S,T)=ex_euclid(c,k);
        
        if(G!=1) {printf("IMPOSSIBLE\n"); continue;}
        
        ll left=max(myfloor(T,c)+1,myfloor(-S,k)+1), right= (1e9-S)/k;
        if(left>right){printf("IMPOSSIBLE\n"); continue;}
        ll ans=S+left*k;
        printf("%lld\n",ans);
    }
    return 0;
}

 

 

 

다음 문제를 풀기 위해서는 페르마 소정리에 대한 이해가 필요하다. 페르마 소정리는 $p$가 소수이고 $p$가 $a$의 약수가 아닐 때, $a^{p-1}≡1 (mod p)$가 성립한다. 페르마 소정리에 대한 이론은 [BO 13172] Σ 우연스럽게 이 문제의 디스크립션에서 친절하게 설명을 해주고 있다. 

디스크립션을 요약하자면 페르마소정리는 모듈러 곱셈에 대한 역원을 정의하는 방법을 제시한다. 곱셈에 대한 역원이란, 어떤 수를 곱해서 1이 나오도록 하는 수이다. 예를 들어 mod11일 때 3의 곱셈에 대한 모듈러 역원은 4이다. (3*4 (mod11) =1 ). 

곱셈에 대한 모듈러 역원이 필요한 이유는 모듈러의 성질 중에서 덧셈, 뺄셈, 곱셈은 모두 분배법칙이 성립하지만 나눗셈(분수)에 대해서는 성립하지 않기 때문이다. 예를 들어

$(a+b)%m = ( (a%m) + (b%m) ) % m$ 은 성립하고,

$(a*b)%m = ( (a%m) * (b%m) ) % m$ 은 성립하지만,

$(a/b)%m = ( (a%m) / (b%m) ) % m$ 은 성립하지 않는다. 

이를 해결하기 위한 수단이 모듈러 역원이고 여기서 페르마 소정리의 힘을 빌린다고 생각하면 되겠다. 페르마 소정리에 따르면 곱셈에 대한 모듈러 역원은 a^(p-2)이다. 

예를 들어 분수 7/3에 대한 모듈러 역원을 구한다고 하자. 즉 $7/3 * x (mod 11) = 1$이 되게하는 정수 x를 찾아보자. 

$7/3$은 $7*3^(-1)$로 쓸 수 있고 나는 결국 3의 역원을 구해야 한다는 뜻이다. mod 11에 대한3의 역원을 위에서처럼 임의의 수를 대입해봐서 구할 수도 있겠지만 이는 시간이 얼마나 걸릴지 장담 못하므로 페르마 소정리를 이용하면 단번에 구할 수 있다. mod 11에서 11은 소수이고, 11은 3의 약수가 아니므로 $3^(11-2)=3^9$는 3의 역원이라고 할 수 있다. 즉 $3*3^9 = 3^10 ≡ 1 (mod 11)$이다. 따라서 $7 * (3^9) ≡ 1 (mod 11)$이 된다. 그럼 나는 여기서 $3^9$를 어떻게 해석할 수 있는가? 내가 생각했을 때 $3^9$는 페르마 소정리의 힘을 빌려서 원래 있었던 $3^(-1)$의 형태를 띠고 있는 환상(illusion)이라고 생각한다. 이 illusion은 영원할 수 없고 $a$값이 $p$의 약수가 아닐 때만 그 효력이 유지된다. 따라서 $a^(p-2)$의 힘은 오직 소수 $p$값의 따라 운명이 결정된다. 만약 $p$값이 큰 소수라면 $[1,p-1]$까지의 $a$는 모두 $p$와 서로소일 것이므로 상당히 넓은 범위에 대해서 페르마소정리를 보장해줄 것이다. 그러나 $p$값이 작은 소수라면 $p$가 $a$의 약수가 될 가능성이 커지므로 페르마소정리를 보장해줄 수 있는 범위도 그만큼 줄어드는 것이다.

 

 

[BOJ 11401] 이항계수 3

 

이 문제에서는 n과 r값이 주어기고 $nCr$ (mod 1e9+7)을 구하는 것이 목표이다. 

이 문제를 풀기 위해서는 이항 계수 식을 분수 형태로 표현해야 한다. 이항 계수 $nCr = n! / (r! * (n-r)! )$으로 표현할 수 있다. 즉 $n! / (r! * (n-r)! ) (mod 1e9+7)$을 구하면 된다. 위에서 설명했다시피 분수 형태는 모듈러 계산이 성립하지 않으므로 역원을 찾아서 곱해야 한다. 

$n! = A , ( r! * (n-r)! ) = B$라고 두면 $A*B^(-1) (mod m)$을 구해야 한다. $B^(-1)$을 페르마 소정리에 의해 $B^(p-2)$로 구할 수 있다. 여기서 p가 10억이 넘으므로 거듭제곱 계산을 할 때 $O(N)$시간복잡도로 하면 시간이 터진다. 따라서 $O(logN)$만에 할 수 있는 방법으로 해야한다.

 

#include<stdio.h>
typedef long long ll;
const int MOD = 1e9 + 7;
const int SIZE = 4000009;
ll fact[SIZE];

ll modpow(ll a,ll x){
    ll res = 1;
    while(x){
        if(x % 2) res = res * a % MOD;
        a = a * a % MOD;
        x /= 2;
    }
    return res;
}

ll getInv(ll x){
    return modpow(x,MOD - 2);
}

ll comb(ll N,ll R){
    ll res =  fact[N] * getInv(fact[R]) % MOD;
    res = res * getInv(fact[N-R]) % MOD;
    return res;
}

int main(){
	int n,r;	scanf("%d %d",&n,&r);
	fact[0] = 1;
	for(int i=1 ; i<=n ; i++) fact[i] = fact[i-1] * i % MOD;
	printf("%lld\n",comb(n,r));
}

 

 

마지막으로 살펴볼 개념은 중국인의 나머지 정리이다.중국인의 나머지 정리에 대한 이론은 이 블로그를 참고했다. 

중국인의 나머지 정리는 아래와 같이 두개 이상의 연립합동식의 해를 찾는 방법이다. 

 

x≡a1(modm1)x≡a2(modm2)

 

여기서 m1,m2라 서로소라는 보장이 없다면 위 식은 다음과 같이 풀 수 있다.

 

x≡m1k0+a1(modlcm(m1,m2))

 

[BOJ 6064] 카잉 달력

 

전형적인 중국인의 나머지 정리 유형 문제이다. 사실 이 문제는 최대공약수,최소공배수로 풀 수 있는 방법도 존재하긴 하지만 지금은 그렇게는 풀지 말자. 문제의 디스크립션은 좀 복잡하게 써놓았지만 실제로는 입력이 m,n,a,b로 주어졌다면 

x≡a (mod m)

x≡b (mod n)

을 풀기만 하면 된다. 여기서 m,n이 서로소가 아닐 수도 있으므로 LCM을 이용하도록 하자.

 

이 문제는 0보다 크면서 가장 작은 x값을 찾기 위해서 적절한 k값을 찾는게 중요한데, 여기서 나는 부등식끼리 꼬여있는게 너무 복잡하다고 생각해서 바이너리 서치의 lower bound로 찾았다. 부등식을 풀기 귀찮을 때는 이것도 좋은 방법일 것 같다.

tuple<int,int,int> exEuclid(int a,int b){
    if(b==0) return {a,1,0};
    int g,x,y;
    tie(g,x,y)=exEuclid(b,a%b);
    return {g,y,x-(a/b)*y};
}

int main(){
    int t;  scanf("%d",&t);

    while(t--){
        int m1,m2,a1,a2;    scanf("%d %d %d %d",&m1,&m2,&a1,&a2);

        int G,x,y;
        tie(G,x,y)=exEuclid(m1,m2);

        if((a2-a1)%G) {printf("-1\n"); continue;}
        x*=(a2-a1)/G;
        y*=(a2-a1)/G;

        ll l=-2e9,r=2e9;

        //lower bound
        while(l+1<r){
            ll mid=(l+r)/2;

            ll k=1LL*x+m2/G*mid;
            ll cur=k*m1+a1;
            if(cur>0) r=mid;
            else l=mid;
        }
        ll k=1LL*x+m2/G*r;
        printf("%lld\n",k*m1+a1);
        continue;

    }

    return 0;
}

 

[BOJ 11661] 해의 개수