hogecoder

つたじろう(Tsuta_J) 競技プログラミングの記録

AtCoder Beginner Contest 020 D: LCM Rush

もう 7 月かあ・・・。

問題概要

 \sum_{i=1}^{N} lcm(i, K) \mod 1,000,000,007 を求めよ。ただし  lcm(a, b) a b の最小公倍数を指す。

制約:  1 \leq N, K \leq 10^{9}

解説

愚直にやるのは無理なので、数学をします。

まず、 lcm(a, b) = a \times b / gcd(a, b) なる関係があります。これを使って与式を変形すると

 \sum_{i=1}^{n} i \times K/gcd(i, K) が求めたいものになります。

ここで、 gcd(i, K) の取りうる値の数は  K の約数の数ぶんしかないことに注目します。この制約下では最大  1344 個であり、比較的少ないです。最大公約数の値が同じものについて、まとめて処理できると嬉しいよなあという気持ちになるため、 gcd の値ごとに分けて考えることにしましょう。

まずは、 gcd(i, K) = 1 を満たす  i について、 lcm(i, K) の総和  S を考えてみましょう。これは、

 S = \sum_{1 \leq i \leq N, gcd(i, K) = 1} lcm(i, K) = K \times \sum_{1 \leq i \leq N, gcd(i, K) = 1} i

となり、条件を満たす  i の総和と  K の積になります。

条件を満たす  i の総和は包除原理によって求めることができます。 K の素因数を列挙し、素因数偶数個の積を約数にもつ数の総和は足し、素因数奇数個の積を約数に持つ数の総和は引き・・・という方針で求められます。(包除原理に関しては他サイトを参照してください)

次に、 gcd(i, K) = x を満たす  i について、 lcm(i, K) の総和  T を考えてみましょう。これは変形すると先ほどの計算と同じ方針でできます。

 gcd(i, K) = x は、  gcd(i/x, K/x) = 1 と同値です。よって、

 T = \sum_{1 \leq i \leq N, gcd(i, K) = x} lcm(i, K) = K \times \sum_{1 \leq i' \leq \lfloor \frac{N}{x} \rfloor, gcd(i', \lfloor \frac{K}{x} \rfloor) = 1} i'

となります。

まとめると、

  1.  K の約数  M それぞれに対して以下を行う
  2.  1 \leq i \leq \lfloor \frac{N}{M} \rfloor かつ  gcd(i, \lfloor \frac{K}{M} \rfloor) = 1 を満たす  i の総和を答えに足す
  3. 2 の操作は包除原理によって高速に処理できる

ということになります。

ソースコード

割と綺麗に書けて個人的には満足です。

int N, K;

// 約数の列挙
vector<int> divisor(int n) {
    vector<int> ret;
    for(int i=1; i*i<=n; i++) {
        if(n % i == 0) {
            ret.push_back(i);
            if(n / i != i) ret.push_back(n/i);
        }
    }
    return ret;
}

// 素因数の列挙
vector<int> pFactorization(int n) {
    vector<int> ret;
    for(int i=2; i*i<=n; i++) {
        if(n % i == 0) {
            ret.push_back(i);
            while(n % i == 0) n /= i;
        }
    }
    if(n != 1) ret.push_back(n);
    return ret;
}

// 1 <= i <= n, gcd(i, m) = 1 なるすべての i の和を求める
int solve(int n, int m) {
    vector<int> vp = pFactorization(m);
    int s = vp.size(), ans = 0;
    rep(bit,0,1<<s) {
        int num = 1;
        rep(i,0,s) if(bit >> i & 1) num *= vp[i];
        int a = (n / num) % MOD;
        int b = (a * (a+1) / 2) % MOD;
        int c = num % MOD;
        if(__builtin_popcount(bit) % 2) {
            int sub = (b * c) % MOD;
            ans = (ans - sub + MOD) % MOD;
        }
        else {
            (ans += b * c) %= MOD;
        }
    }
    return ans;
}

signed main() {
    cin >> N >> K;
    vector<int> div = divisor(K);
    int ans = 0;
    rep(i,0,div.size()) {
        (ans += solve(N/div[i], K/div[i]) * K) %= MOD;
    }
    cout << ans << endl;
    return 0;
}

解けるまでだいぶ時間かかったので、これは要復習かも。