AtCoder Beginner Contest 020 D - LCM Rush

問題記録 AtCoder ABC ABC #020

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

問題概要

整数N, Kが与えられる。 \sum_{i=1}^{N}{LCM(i, K)}を求めよ。

1 <= N, K <= 109

解法

（要点：「約数にxをもつ集合」から「約数に2xを持つ集合」, 「約数に3xを持つ集合」, …を引いていけば「GCDがxの集合」を出せる）

LCM(n, K) = nK / GCD(n, K) である。ここでGCDの定義より、GCD(n, K)が取る値は必ずKの約数のうちのどれかになる。Kの約数ごとに場合分けして求める値を計算することにすると、K / GCD(n, K)が定数になるためΣの外にくくりだすことができる。以上よりこの問題は結局、Kのすべての約数xについて「KとのGCDがxであるような整数i (1 <= i <= N)の和 × K × x」を求め、それらを足し合わせる問題と言い換えることができる。なお109以下の数であれば約数の個数は最大でも1000ちょっとくらいしかないので、これは十分に列挙可能である。

残る問題は「KとのGCDがxであるような整数i (1 <= i <= N)の和」をどう計算するかということである。ここで一旦GCDという条件を緩めて「約数にxを持つような整数i (1 <= i <= Nの集合」がどうなるかを考えてみると、単純な等差数列になる。例えばx = 1のときは1, 2, …, Nになるし、x = 2のときは2, 4, …となる（形式的には初項x, 公差x, 項数[N/x]の等差数列）。ここで「約数にxをもつ集合」から「約数に2xを持つ集合」「約数に3xを持つ集合」…を引いていくと、残ったものが「GCDがxである集合」になる（ただし2x, 3x等がKの約数でない場合は飛ばす）。であるから、Kの約数xの大きい方から見ていって「約数にxを持つ集合」（＝等差数列）の和を求め、Kの約数であるような2x, 3x, …, の場合の値を引いていけばKのすべての約数xについてGCD(i, K) = xの場合の値が出せる。最後にそれらを足し合わせたものが答え。

等差数列の和は公式に当てはめればO(1)で、Kの約数の個数は高々1000個くらいになるので十分に間に合う。

感想

前のこどふぉで解いた問題が「LCM Rushっぽい」と言われていたのでそろそろ解けるかなと思ってやってみたら、解けた。これでABC-Dが全埋めということになって嬉しい。

コード (D言語)

import std.stdio, std.array, std.string, std.conv, std.algorithm;
import std.typecons, std.range, std.random, std.math, std.container;
import std.numeric, std.bigint, core.bitop;
 
immutable long MOD = 10^^9 + 7;
 
void main() {
    auto s = readln.split.map!(to!long);
    auto N = s[0];
    auto K = s[1];
 
    long[] F;
    for (long i = 1; i * i <= K; i++) {
        if (i * i == K) F ~= i;
        else if (K % i == 0) F ~= i, F ~= K / i;
    }
    F.sort();
 
 
    auto ans = new long[](F.length);
    auto sums = new long[](F.length);
    
    for (int i = F.length.to!int - 1; i >= 0; i--) {
        long n = N / F[i];
        sums[i] = n * (2 * F[i] % MOD + (n - 1) * F[i] % MOD) % MOD * powmod(2, MOD-2, MOD) % MOD;
 
        for (int j = i + 1; j < F.length; j++) {
            if (F[j] % F[i] == 0) {
                sums[i] -= sums[j] % MOD;
                sums[i] = (sums[i] + MOD) % MOD;
            }
        }
        
        ans[i] = K * sums[i] % MOD * powmod(F[i], MOD-2, MOD) % MOD;
    }
 
    long ansans = 0;
    foreach (a; ans) ansans = (ansans + a) % MOD;
    ansans = (ansans + MOD) % MOD;
    ansans.writeln;
}
 
 
long powmod(long a, long x, long m) {
    long ret = 1;
    while (x) {
        if (x % 2) ret = ret * a % m;
        a = a * a % m;
        x /= 2;
    }
    return ret;
}