디오판토스 방정식의 해 구하기

ax+by=c (a, b, c는 정수) 형태의 방정식에서 정수해를 구해야 하는 걸 디오판토스 방정식이라고 부른다.

 

이 해는 gcd(a,b)|c일 때에만 존재한다

 

{ax+by|x,y∈ℤ}={gcd(a,b)n|n∈ℤ} 이기 때문인데, 베주 항등식과 관련이 있다. 글주제에서는 조금 벗어난 내용이니까 자세히 서술하지 않음

 

여기서 gcd(a,b)는 유클리드 호제법 Euclidean Algorithm으로 구할 수 있다.

유클리드 호제법을 사용한 예시 : 챗 지피티로 생성함

자.

여기서 유클리드 호제법을 사용하는 동안 등장한 모든 나머지들 (142, 21, 16, 5 ...) 은 전부 2019과 1877의 linear combination으로 나타낼 수 있다.

 

그러니까 이 정수들을 차라리 순서쌍 (x, y)로 나타내도록 하자. (x, y)는 실제로는 정수 2019x+1877y를 표현하는 거다.

 

예를 들어 142 = 2019 - 1877이므로 (1, -1)일 것이다.

그럼 21 = 1877 - 142 * 13 인데, 1877은 (0, 1)이고 142는 (1, -1)이니까 21 = (0,1) - (1, -1) * 13 = (-13, 14)이다.

같은 방식으로 16, 5 ... 쭉 구해서 gcd(2019,1877)도 순서쌍으로 나타낼 수 있다.

그럼 2019x+1877y=c라는 방정식을 보자.

c는 gcd(2019,1877)의 배수이다.

 

그러니 만약 c=N*gcd(2019,1877)이고,

gcd(2019,1877)=(u,v)이라면

 

위 디오판토스 방정식의 해는 (N*u,N*v)가 되는 거다.

 

이 원리로 C++로 디오판토스 방정식의 해를 찾아주는 코드를 짜보았다.

유의할 점!! 디오판토스의 해는 유일하지 않다.

따라서 해를 "전부" 찾아주는 방법론이 아니다!!

해 중 1개를 찾아주는 거다.

//Finding ONE solution to the diophantine equation ax + by = c

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;
typedef long long ll;

//basically, all the integers that appear can be expressed as ax + by for some integers x and y
//so, instead of using the number itself, we can use the pair (x, y) to represent it

class Integer 
{
    public:
    ll x, y;
    static ll a, b;
    Integer(ll x = 0, ll y = 0) : x(x), y(y) {}

    Integer operator+(const Integer &other) const 
    {
        return Integer(x + other.x, y + other.y);
    }

    Integer operator-(const Integer &other) const 
    {
        return Integer(x - other.x, y - other.y);
    }
    Integer operator*(ll scalar) const 
    {
        return Integer(x * scalar, y * scalar);
    }

    ll value()
    {
        return a * x + b * y;
    }
};
ll Integer::a = 0;
ll Integer::b = 0;

int main()
{
    ll a, b, c;
    cin >> a >> b >> c;
    Integer::a = a;
    Integer::b = b;
    // input a, b, c in ax + by = c
    Integer A(1, 0), B(0, 1);

    while(1)
    {
        ll va = A.value(), vb = B.value();
        if (va < vb)
        {
            swap(A, B);
            swap(va, vb);
        }
        ll q=va/vb;
        ll r=va%vb;
        if (r==0)
        {
            break;
        }

        A=A-B*q;
    }
    if (c % B.value() != 0)
    {
        cout << "no solution\n";
        return 0;
    }
    Integer C = B * (c/B.value());
    cout << C.x << " " << C.y << "\n"; 
    cout<< Integer::a <<" * "<< C.x<< " + " << Integer::b <<" * "<< C.y <<" = "<< c <<"\n";
	return 0;
}