mayoko’s diary

プロコンとかいろいろ。

ABC#20D問題:LCM Rush

とりあえず本番部分点取れたのは良かった。満点解法は答えみてなんとなく理解してからもめっちゃ混乱して3時間くらいかかった。

問題:http://abc020.contest.atcoder.jp/tasks/abc020_d

解法:Kの各gcdの値gごとに(GCD(n, K)=g)となるnの和を計算する。これを満たす各nについて、LCM(n, K) = n * (K/g)となるので、この和sumを計算できれば、各gごとにsum * (K/g)を計算することで答えが得られる。
じゃあ各gごとにどうやって和を計算するのってことになるが、GCD(n,K) = gのとき、少なくともnはgの倍数であるので,GCD(n, K) = GCD(g*m, K) = g*GCD(m, K/g)となる(mは1〜N/g)。よって、GCD(n,K)=gとなるnの和は、mとK/gで互いに素になるものの和を計算すれば良い。
じゃあ今度はmとK/gで互いに素になるものの和はどう計算できるんですかってことになる。これは包除原理的にやっていけば良い。まず単純和を計算して、次に2で割り切れるものの和を計算して引いて、3で割り切れるものの和を計算して引いて、6で割り切れるものの和を計算して足して、...というような感じ。これらはすべて等差数列になるのでO(1)で計算できる。
以下ソースコード

#include <bits/stdc++.h>

#define rep(i, n) for (int i = 0; i < (int)n; i++)

using namespace std;
typedef long long ll;
typedef pair<ll, int> pli;

const int MAX = 100010;
bool isPrime[MAX];

void primeInit() {
    for (int i = 2; i < MAX; i++) isPrime[i] = true;
    for (int i = 2; i < MAX; i++) {
        if (isPrime[i]) {
            for (int j = 2; j * i < MAX; j++) {
                isPrime[i*j] = false;
            }
        }
    }
}

vector<pli> bunkai(ll n) {
    vector<ll> div;
    for (ll i = 2; i * i <= n; i++) {
        if (n % i == 0 && isPrime[i]) div.push_back(i);
    }
    vector<pli> ret;
    for (auto el : div) {
        int cnt = 0;
        while (n % el == 0) {
            cnt++;
            n /= el;
        }
        ret.push_back(pli(el, cnt));
    }
    if (n != 1) ret.push_back(pli(n, 1));
    return ret;
}

const ll MOD = 1000000007;
ll N, K;
ll ans = 0;

void dfs(int cur, ll g, vector<pli> pf) {
    if (cur == (int)pf.size()) {
        ll M = N / g;
        int size = pf.size();
        int S = 1<<(size);
        ll sum = 0;
        for (int s = 0; s < S; s++) {
            int cnt = 0;
            ll start = 1;
            bool ng = false;
            for (int i = 0; i < size; i++) {
                if (((s>>i) & 1)) {
                    if (pf[i].second > 0) {
                        cnt++;
                        start *= pf[i].first;
                    } else {
                        ng = true;
                        break;
                    }
                }
            }
            if (ng) continue;
            ll num = M / start;
            ll p = num * (num+1) / 2 * start;
            p %= MOD;
            if (cnt % 2 == 0) {
                sum += p;
            } else {
                sum -= p;
            }
            sum %= MOD;
            if (sum < 0) sum += MOD;
        }
//        printf("g is %lld sum is %lld\n", g, sum);
        ans += sum*K;
        ans %= MOD;
    } else {
        ll now = 1;
        for (int i = 0; i <= pf[cur].second; i++) {
            pf[cur].second -= i; // 残る素因数
            dfs(cur+1, g*now, pf);
            pf[cur].second += i;
            now *= pf[cur].first;
        }
    }
}

int main(void) {
    primeInit();
    cin >> N >> K;
    vector<pli> pf = bunkai(K); // prime factor
//    cout << "pf size is " << pf.size() << endl;
    dfs(0, 1, pf);
    cout << ans << endl;
    return 0;
}