はまやんはまやんはまやん

hamayanhamayan's blog

Divide Both [AtCoder Beginner Contest 206(Sponsored by Panasonic) E]

https://atcoder.jp/contests/abc206/tasks/abc206_e

前提知識

解説

https://atcoder.jp/contests/abc206/submissions/23623547

しょっぱなから発想の転換が必要な問題。
類題を多く解いているので、すぐ思いついたが、そうでないと難しい。

発想の転換

xとyを列挙して、条件を満たすかどうかを確認していくような流れでは今回は解けない。
そうではなく、最大公約数がgとなるようなx,yの組がどれだけあるかという目線で解いていく。
最大公約数gを全探索して、条件を満たす(x,y)の組がどれだけあるかを計算することで問題を考えていこう。

少しだけ簡単な問題を解く

条件の少し緩くして考えてみよう。
条件を「g != 1」だけにした状態で問題が解けるかを考えてみる。
この問題が解ければ、大体解法は理解したようなもの。

数値の上限が106なので、同様にgの上限も106である。
これは全探索できそうである。
gを全探索して、最小公倍数がgとなる個数を数え上げよう。

最小公倍数がgとなるためには少なくとも、x,yはgの倍数である必要がある。
よって、[L,R]の中のgの倍数の個数、R/g-(L-1)/gがx,yの候補数であり、
(R/g-(L-1)/g)2が答えの組み合わせのように見える。
だが、これではgの倍数同士の組であることは言えるが、x=2g, y=2gのような最小公倍数がgの倍数となってしまう組も含まれる。
(R/g-(L-1)/g)2の中で適さないのは最小公倍数が2g, 3g, 4g, ...のものであり、これを取り除くいい方法がある。
約数系包除原理である。

約数系包除原理

エラトステネスの篩というか、調和級数的計算量というか、
g=1~Nの全てに対して2g, 3g, ...をN以下まで探索するようなループを書くと
全体の計算量はO(NlogN)になるというものがある。
何を言っているのか分からないという場合は、エラトステネスの篩のアルゴリズムを学ぶとどういうことか分かると思う。
コードベースで説明すると「rep(i,1,N) for(j=i;j<=N;j+=i) というループ構造はO(NlogN)で行える」ということ。

ここで言っていることは丁度、さっきやりたかったことである。
以下のような配列を用意して、適さない2g, 3g, 4g, ...を取り除くようなコードを書こう。

cnt[g] := 最小公倍数がgであるときの(x,y)の組合せ

実装の注意点としては、cnt[g]を計算するときにはcnt[2g], cnt[3g], ...を引く必要があるので、gは降順に計算していくこと。
これで、基本的にはcnt[g] = (R/g-(L-1)/g)2なのだが、これでは最小公倍数が2g, 3g, とかが混ざっているので、
gの倍数がR以下になるまでcnt[2g], cnt[3g], ...を引けば、ちょうど最小公倍数がgになるときの組み合わせが求められる。

元々の問題は?

後はcntの総和を取れば、少しだけ簡単にした問題は解ける。
元々の問題では、x != g, y != gという条件がついているので、別途、x=gまたはy=gのパターンを引いてやろう。
こうなるのは[L,R]にgが含まれているときだけで、

  • x=g、かつ、yはgの倍数(R/g-(L-1)/g個)を引く
  • y=g、かつ、xはgの倍数(R/g-(L-1)/g個)を引く
  • これだとx=g, y=gが2回引かれているので1回引きにするために+1する

のように微調整してやることで元々の問題も解ける。

int L, R;
ll cnt[1010101];
//---------------------------------------------------------------------------------------------------
void _main() {
    cin >> L >> R;

    rrep(g, 1010100, 2) {
        cnt[g] = 1LL * (R / g - (L - 1) / g) * (R / g - (L - 1) / g);
        int x = 2 * g;
        while (x <= R) {
            cnt[g] -= cnt[x];
            x += g;
        }
    }

    ll ans = 0;
    rep(g, 2, 1010101) {
        ans += cnt[g]; 
        if (L <= g && g <= R) {
            ans -= R / g - (L - 1) / g;
            ans -= R / g - (L - 1) / g;
            ans++;
        }
    }
    cout << ans << endl;
}