CF 900 D - Unusual Sequences

問題

Problem - D - Codeforces

問題概要

数列について、かつを満たすような数列が何種類存在するか、で割った余りを求めよ。

• 

アイデア

まず、数列の全てのgcdがであるということは、数列の全ての要素はの倍数であるということになる。よってその和もの倍数でないと数列は存在し得ない。

つまり、がで割り切れない時、答えは0であり、割り切れる時、条件はかつと言い換えられる。

では、数列の全要素の和がで全要素のgcdが1であるような数列の個数をで表す。つまり、この問題で求めたいものはである。

また、数列の全要素の和がである(gcdは問わない)ような数列の個数をとすると、とには以下のような関係がある：

上の式からだけをΣの外に出して移項すると、

また、は直接計算することが出来る。個の1が並んでいる状況を想定して、個の間にそれぞれカンマを入れるか否かによって全ての数列を対応付ける事できる。よって、である。以上のことから、

となる。の約数の個数はであるから、後はこれを再帰的にメモ化しながら計算すれば十分間に合う(本当に十分速かった(15ms)(submission))。

メビウスの反転公式による高速化

この式に注目しよう：

の約数を列挙して、の和を考えているが、がの約数全てを動く時、によって列挙されるものもまたの約数全てということになるので、以下の用に書いても差し支えないことに気づく：

上の式にはメビウスの反転公式が適用できる：

これによって、の約数を列挙し、直接を計算できる(solve2で実装)。

実装(C++)

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define all(x) (x).begin(),(x).end()
#define pb push_back
#define fi first
#define se second

const ll mod = 1e9+7;

ll mod_pow(ll x, ll n){
    ll ret = 1;
    while(n){
        if(n&1) (ret*=x)%=mod;
        (x*=x)%=mod;
        n>>=1;
    }
    return ret;
}

// nの約数におけるメビウス関数の値 O(sqrt(n))
map<int,int> moebius(int n){
    vector<int> primes;
    // nの素因数を列挙
    for(int i=2; i*i<=n; ++i){
        if(n%i==0) primes.push_back(i);
        while(n%i==0) n/=i;
    }
    if(n>1) primes.push_back(n);

    map<int,int> ret;
    int SZ = primes.size();
    // 2^SZ通りを試す(nの約数の個数よりは少ない)
    rep(i,1<<SZ){
        int mu = 1, d = 1;
        rep(j,SZ){
            if(i>>j&1){
                mu *= -1;
                d *= primes[j];
            }
        }
        ret[d] = mu;
    }
    return ret;
}

map<ll,ll> m1;
ll solve1(ll n){
    if(m1.count(n)) return m1[n];

    ll ret = mod_pow(2,n-1);
    for(int i=1; i*i<=n; ++i){
        if(n%i==0){
            if(i!=n) ret = (ret - solve1(i) + mod)%mod;
            if(n/i!=n && n/i!=i) ret = (ret - solve1(n/i) + mod)%mod;
        }
    }
    return m1[n] = ret;
}

ll solve2(ll n){
    ll ret = 0;
    map<int,int> m = moebius(n);
    for(int i=1; i*i<=n; ++i){
        if(n%i==0){
            (ret += m[n/i]*mod_pow(2,i-1)+mod) %= mod;
            if(n/i!=i) (ret += m[i]*mod_pow(2,n/i-1)+mod) %= mod;
        }
    }
    return ret;
}

int main(){
    ll x,y;
    cin >>x >>y;

    ll ans = 0;
    if(y%x==0){
        // ans = solve1(y/x);
        ans = solve2(y/x);
    }
    cout << ans << endl;
    return 0;
}