ユークリッドの互除法から中国剰余定理ライブラリの実装と、それを用いる問題2つ

1. 内容

1-1. 概要

・ユークリッドの互除法とは？
・中国剰余定理とは？
ユークリッドの互除法を用いて、中国剰余定理ライブラリを作成する過程を説明。
実際にライブラリを使って解いた問題とそのコードも最後に記載。

1-2. 使用言語, 環境

C++(GCC 9.2.1)

2. ユークリッドの互除法、中国剰余定理

以下記事に詳しい解説があったので、省略。

ユークリッドの互除法

中国剰余定理

3. 実装

コード内コメントアウトの形で説明。
流れとしては、

crt1で、ax-by=gcd(a,b)を満たす最小の正のxと、それに対応するyを求める
→crt2でa mod bかつ、c mod dを満たす、val mod (lcm(b,d))となるval, lcmを返す
→crtで、初期値をa=0, b=1とし、cとdに値を順に入れていって再帰でn個の場合に対応させる

という感じ。

ll mod(ll x, ll y) { // x%yをxが負の場合でも使えるように拡張
    if (x<0) x+=(abs(x)/y)*y+y;
    return x%y;
}

vector<ll> crt1(ll x, ll y) {
    //x>=yで入れる！！(x<yでも使えるかもしれないがそこのエラーチェックをしていないので)
    // x*s-y*t=gcd(x,y)となるs,t,gcdを返す
    if (x%y==0) return {1,x/y-1,y};

    // xs-yt=v[2](gcd)を満たす
    // →x=(x/y)*y+(x%y)なので、
    // (x%y)=x-(x/y)*y・・・(1)
    // y*v[0]-(x%y)*v[1]=v[2]に(1)を代入して、
    // y*v[0]-(x-(x/y)*y)*v[1]=v[2]となる。これを整理して、
    // -v[1]*x+(v[0]+(x/y)*v[1])*y=v[2]
    // となるので、s=-v[1]のときにs,tの整数解が存在することがわかる。
    auto v = crt1(y,x%y);
    ll s = mod(-v[1],y/v[2]);
    ll t = (x*s-v[2])/y;
    return {s,t,v[2]};
}

P crt2(ll a, ll b, ll c, ll d) {
    // a mod bかつ、c mod dとなるval mod(lcm(b,d))

    // c-aがgcd(b,d)で割り切れない場合は、条件を満たすものは存在しない
    if (abs(c-a)%__gcd(b,d)!=0) return {INF,-1};

    // b>dの形に統一しておく
    if (b<d) {swap(b,d); swap(a,c);}
    
    // まず、crt1をもちいてbx-dy=gcd(b,d)となるものを求め、
    auto v = crt1(b,d);

    // それを((c-a)/gcd(b,d))倍することで、bx-dy=c-aを満たすxが求まる
    // (ここでmod dすることでこれを満たす最小のxになる)
    ll val = mod(v[0]*((c-a)/v[2]),d);

    // 今、xを求めたので、この値からxb+aをもとめる。
    val = val*b+a;
    ll lcm = b*d/v[2];

    // 多分必要ないが、一応modをとっている(いらなさそうだが、これを消してテストしていないので)
    val = mod(val,lcm);

    return {val,lcm};
}

P crt(vector<ll> rem, vector<ll> div) {
    // CRT(crt2のn個拡張版) remに余り、divに法をindexが一致するように入れる
    ll x = 0; ll y = 1; P v;
    // x,yにたいしてrem, divにあるものを使って、crt2で行った操作を加えていく
    REP(i,rem.size()) {
        if (y<div[i]) {swap(x,rem[i]); swap(y,div[i]);}
        v = crt2(x,y,rem[i],div[i]);
        if (v.second==-1) return {-1,-1};
        x = v.first; y = v.second;
    }
    return v;
}

4. 実際に使った問題

1. yukicoder No.186　中華風

こちらは互除法そのままである。

https://yukicoder.me/problems/no/186
#include <bits/stdc++.h>
using namespace std;
typedef long long ll; //int:2*10**9
typedef long double ld;
typedef pair<ll,ll> P;
#define REP(i,n) for(ll i = 0; i<(ll)(n); i++)
#define FOR(i,a,b) for(ll i=(a);i<=(b);i++)
#define FORD(i,a,b) for(ll i=(a);i>=(b);i--)
#define vec2(name,i,j,k) vector<vector<ll>> name(i,vector<ll>(j,k))
#define vec3(name,i,j,k,l) vector<vector<vector<ll>>> name(i,vector<vector<ll>>(j,vector<ll>(k,l)))
#define pb push_back
#define MOD 1000000007 //998244353
#define PI 3.141592653
#define INF 8000000000000000000 //14
//cin.tie(0);cout.tie(0);ios::sync_with_stdio(false);
ll mod(ll x, ll y) {
    if (x<0) x+=(abs(x)/y)*y+y;
    return x%y;
}

vector<ll> crt1(ll x, ll y) {
    //x>=yで入れる！！
    // x*s-y*t=gcd(x,y)となるs,t,gcdを返す
    if (x%y==0) return {1,x/y-1,y};
    auto v = crt1(y,x%y);
    ll s = mod(-v[1],y/v[2]);
    ll t = (x*s-v[2])/y;
    return {s,t,v[2]};
}

P crt2(ll a, ll b, ll c, ll d) {
    // a mod bかつ、c mod dとなるval mod(lcm(b,d))
    if (abs(c-a)%__gcd(b,d)!=0) return {INF,-1};
    if (b<d) {swap(b,d); swap(a,c);}
    auto v = crt1(b,d);
    ll val = mod(v[0]*((c-a)/v[2]),d);
    val = val*b+a;
    ll lcm = b*d/v[2];
    val = mod(val,lcm);
    return {val,lcm};
}
P crt(vector<ll> rem, vector<ll> div) {
    // CRT(crt2のn個拡張版)
    ll x = 0; ll y = 1; P v;
    REP(i,rem.size()) {
        if (y<div[i]) {swap(x,rem[i]); swap(y,div[i]);}
        v = crt2(x,y,rem[i],div[i]);
        if (v.second==-1) return {-1,-1};
        x = v.first; y = v.second;
    }
    return v;
}

int main(){
    vector<ll> x(3); vector<ll> y(3);
    REP(i,3) cin >> x[i] >> y[i];
    auto v = crt(x,y);
    if (v.first!=0) cout << v.first << endl;
    else cout << v.second << endl;
}

2.キャディプログラミングコンテスト2021(ABC193) E - oversleeping

https://atcoder.jp/contests/abc193/tasks/abc193_e
1つめよりは考察がいるが、主に中国剰余定理。
自分は定理の実装が間に合わず、無限にWAした。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll; //int:2*10**9
typedef long double ld;
typedef pair<ll,ll> P;
#define REP(i,n) for(ll i = 0; i<(ll)(n); i++)
#define FOR(i,a,b) for(ll i=(a);i<=(b);i++)
#define FORD(i,a,b) for(ll i=(a);i>=(b);i--)
#define vec2(name,i,j,k) vector<vector<ll>> name(i,vector<ll>(j,k))
#define vec3(name,i,j,k,l) vector<vector<vector<ll>>> name(i,vector<vector<ll>>(j,vector<ll>(k,l)))
#define pb push_back
#define MOD 1000000007 //998244353
#define PI 3.141592653
#define INF 5000000000000000000 //14
//cin.tie(0);cout.tie(0);ios::sync_with_stdio(false);
ll mod(ll x, ll y) {
    if (x<0) x+=(abs(x)/y)*y+y;
    return x%y;
}

vector<ll> crt1(ll x, ll y) {
    //x>=yで入れる！！
    // x*s-y*t=gcd(x,y)となるs,t,gcdを返す
    if (x%y==0) return {1,x/y-1,y};
    auto v = crt1(y,x%y);
    ll s = mod(-v[1],y/v[2]);
    ll t = (x*s-v[2])/y;
    return {s,t,v[2]};
}

P crt2(ll a, ll b, ll c, ll d) {
    // a mod bかつ、c mod dとなるval mod(lcm(b,d))
    if (abs(c-a)%__gcd(b,d)!=0) return {INF,-1};
    if (b<d) {swap(b,d); swap(a,c);}
    auto v = crt1(b,d);
    ll val = mod(v[0]*((c-a)/v[2]),d);
    val = val*b+a;
    ll lcm = b*d/v[2];
    val = mod(val,lcm);
    return {val,lcm};
}
P crt(vector<ll> rem, vector<ll> div) {
    // CRT(crt2のn個拡張版)
    ll x = 0; ll y = 1; P v;
    REP(i,rem.size()) {
        if (y<div[i]) {swap(x,rem[i]); swap(y,div[i]);}
        v = crt2(x,y,rem[i],div[i]);
        if (v.second==-1) return {-1,-1};
        x = v.first; y = v.second;
    }
    return v;
}

void solve() {
    ll x, y, p, q; cin >> x >> y >> p >> q;
    ll b = 2*(x+y); ll d = p+q;
    ll ans = INF;
    FOR(a,x,x+y-1) FOR(c,p,p+q-1) {
        auto P = crt2(a,b,c,d); ll ansc = P.first;
        ans = min(ans,ansc);
    }
    if (ans==INF) cout << "infinity" << endl;
    else cout << ans << endl;
}

int main(){
    ll t; cin >> t;
    REP(i,t) solve();
    return 0;
}