【CF1174E】Ehab and the Expected GCD Problem

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9+7;
const int N = 1e6+5;
int n;
ll siz[N];
ll d[2][25][25],pow2[25],pow3[25];
void ycl(){
    int o=0;
    pow2[0] = pow3[0] = 1;
    for(int i=1;i<=n;i++){
        siz[i] = n/i;
    }
    for(int i=1;i<=20;i++){
        pow2[i] = pow2[i-1] * 2 % mod;
    }
    for(int i=1;i<=13;i++){
        pow3[i] = pow3[i-1] * 3 % mod;

    }
}
bool flag = 0;
int main() {
    cin>>n;
    if(n <= 2) {
        cout<<1<<endl;
        return 0;
    }
    ycl();
    int nn = 0;
    while((1<<(nn+1))<=n) nn++;
    d[1][nn][0] = 1;
    if(pow2[nn-1]*3<=n) d[1][nn-1][1] = 1;
    int o = 1;
    ll ans = 0;
    for(int len=2;len<=n;len++){
        bool flag = 0;
        o = o^1;
        memset(d[o],0,sizeof(d[o]));
        for(int i=0;pow2[i]<=n;i++){
            for(int j=0;j<=1&&pow2[i]*pow3[j]<=n;j++) {
                int now = pow2[i]*pow3[j];
                d[o][i][j] += d[o^1][i][j]*(siz[now]-(len-1))%mod;
                d[o][i][j] %= mod;
                if(i>=1){
                    d[o][i-1][j] += d[o^1][i][j]*(siz[pow2[i-1]*pow3[j]]-siz[now])%mod;
                }
                if(j>=1){
                    d[o][i][j-1] += d[o^1][i][j]*(siz[pow2[i]*pow3[j-1]]-siz[now])%mod;
                }
            }
        }
    }
    cout<<d[o][0][0]<<endl;
}