けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

ACL Contest 1 B - Sum is Multiple (青色, 600 点)

AtCoder AtCoder600点 数学(整数問題) 入力が定数個 テク:約数の個数は少ない 素因数分解 全探索:ビット全探索 全探索 拡張Euclidの互除法 中国剰余定理 枝刈り 互いに素 オーバーフロー ある量を固定して考える 青色diff 【問題集】整数変数の式で表された条件を扱う探索

中国剰余定理使う問題楽しいよね！

問題へのリンク

問題概要

正の整数 N が与えられる。以下の条件を満たす最小の正の整数 K を求めよ。

  • 1 + 2 + \dots + K が N の倍数である

制約

  • 1 \le N \le 10^{15}

考えたこと

まず、 1 + 2 + \dots + K = \frac{1}{2}K(K+1) なので、問題の条件は次と同値になる。

  • K(K+1) が 2N の倍数である

ここで改めて N を 2N で置き換えておく。さて、Euclid の互除法より

{\rm GCD}(K, K+1) = {\rm GCD}(K, 1) = 1

であることから、 K と K+1 は互いに素となる (これは受験数学ではめっちゃよく見る典型)。したがって、 N を素因数分解したとき、各素因子は K と K+1 に別々に振り分けられることになる。つまり

  • K は p_{1}^{e_{1}} \dots p_{a}^{e_{a}} の倍数
  • K+1 は q_{1}^{e_{1}} \dots q_{b}^{e_{b}} の倍数
  • N = p_{1}^{e_{1}} \dots p_{a}^{e_{a}} q_{1}^{e_{1}} \dots q_{b}^{e_{b}}

という風になる。 K と K+1 に素因子をそれぞれどのように振り分けるかについては bit 全探索する！！

なお、

2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29 * 31 * 37 * 41 = 304250263527210

より、登場する素因子の個数は最大でも 13 個になる。よって bit 全探索によって探索する場合の数は最大でも 2^{13} = 8192 通りとなって十分小さい。以後、

  • K は A の倍数
  • K+1 は B の倍数

を満たす最小の正の整数 K を求める問題を考えていく ( AB = N)。

解法 (1)：中国剰余定理

ライブラリさえあればこれが最も楽かなと思う。条件は

  • K ≡ 0 \pmod A
  • K ≡ -1 \pmod B

と言い換えることができる。中国剰余定理より、ある整数 r が存在して、これは

  • K ≡ r \pmod N

と同値になることが保証される。そのような r を求めればよい。この作業に要する計算量は O(\log N)。

#include <bits/stdc++.h>
using namespace std;
template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; }
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }

vector<long long> prime_factorize(long long n) {
    vector<long long> res;
    for (long long p = 2; p * p <= n; ++p) {
        if (n % p != 0) continue;
        int num = 0;
        long long val = 1;
        while (n % p == 0) { ++num; n /= p; val *= p; }
        res.push_back(val);
    }
    if (n != 1) res.push_back(n);
    return res;
}

long long extGcd(long long a, long long b, long long &p, long long &q) {
    if (b == 0) { 
        p = 1, q = 0; 
        return a; 
    }
    long long d = extGcd(b, a % b, q, p);
    q -= a / b * p;
    return d;
}

pair<long long, long long> ChineseRem(const vector<long long> &vr, const vector<long long> &vm) {
    if (vr.empty() || vm.empty()) return make_pair(0, 1);
    long long R = vr[0], M = vm[0];
    for (int i = 1; i < (int)vr.size(); ++i) {
        long long p, q, r = vr[i], m = vm[i];
        if (M < m) swap(M, m), swap(R, r); // prevent overflow
        long long d = extGcd(M, m, p, q); // p is inv of M/d (mod. m/d)
        if ((r - R) % d != 0) return make_pair(0, -1);
        long long md = m / d;
        long long tmp = (r - R) / d % md * p % md;
        R += M * tmp, M *= md;
    }
    R %= M;
    if (R < 0) R += M;
    return make_pair(R, M);
}

int main() {
    long long N;
    cin >> N;
    N *= 2;
    auto pf = prime_factorize(N);
    int M = pf.size();
    long long res = N - 1;
    for (int bit = 0; bit < (1<<M); ++bit) {
        long long A = 1;
        for (int i = 0; i < M; ++i) if (bit & (1<<i)) A *= pf[i];
        long long B = N / A;
        auto p = ChineseRem({0, B-1}, {A, B});
        if (p.first > 0) chmin(res, p.first);
    }
    cout << res << endl;
}

解法 (2)：拡張 Euclid の互除法

K が A の倍数で、 K + 1 が B の倍数であるという条件から、整数 q, r が存在して、

  • K = Aq
  • K + 1 = Br

と表せる。よって、

A(-q) + Br = 1

という関係が成り立つ。これを満たすような q, r はまさに拡張 Euclid の互除法によって求めることができる。この作業に要する計算量は O(\log N)。

#include <bits/stdc++.h>
using namespace std;
template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; }
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }

vector<long long> prime_factorize(long long n) {
    vector<long long> res;
    for (long long p = 2; p * p <= n; ++p) {
        if (n % p != 0) continue;
        int num = 0;
        long long val = 1;
        while (n % p == 0) { ++num; n /= p; val *= p; }
        res.push_back(val);
    }
    if (n != 1) res.push_back(n);
    return res;
}

long long extGcd(long long a, long long b, long long &p, long long &q) {
    if (b == 0) { 
        p = 1, q = 0; 
        return a; 
    }
    long long d = extGcd(b, a % b, q, p);
    q -= a / b * p;
    return d;
}

int main() {
    long long N;
    cin >> N;
    N *= 2;
    auto pf = prime_factorize(N);
    int M = pf.size();
    long long res = N - 1;
    for (int bit = 0; bit < (1<<M); ++bit) {
        long long A = 1;
        for (int i = 0; i < M; ++i) if (bit & (1<<i)) A *= pf[i];
        long long B = N / A;
        long long q, r;
        extGcd(A, B, q, r);
        long long tmp = (A * -q) % N;
        if (tmp <= 0) tmp += N;
        chmin(res, tmp);
    }
    cout << res << endl;
}