[HAOI2011] Problem b

复习一下莫比乌斯反演

首先很显然用一下容斥把它转化成求 (ans=sum_{i=1}^a sum_{j=1}^b [{gcd(i,j)=d}])

我们可以定义 f(d) 和 F(d) 如下：

(f(d)=sum_{i=1}^Nsum_{j=1}^M[gcd(i,j)=d])

(F(d)=sum_{i=1}^Nsum_{j=1}^M[d|gcd(i,j)])

发现

(sum_{n|d}f(d)=F(n)=lfloorfrac Nn floorlfloorfrac Mn floor)

莫比乌斯反演，得到：

(f(n)=sum_{n|d}mu(lfloorfrac dn floor)F(d))

于是
(ans=f(d)=sum_{d|p}mu(lfloorfrac pd floor)F(p))

换个元
(ans=sum_{p'}mu(p')F( p' d ) =sum_{p'=1}^{min(lfloorfrac Nd floor,lfloorfrac Md floor)}mu(p')lfloorfrac N{p'd} floor lfloor frac M{p'd} floor)

把 (p') 写作 (p) ，得到

(ans=sum_{p}mu(p)F( p d ) =sum_{p=1}^{min(lfloorfrac Nd floor,lfloorfrac Md floor)}mu(p)lfloorfrac N{pd} floor lfloor frac M{pd} floor)

#include <bits/stdc++.h>
using namespace std;

#define int long long
const int N = 1000005;
const int MAXN = 1000005;

bool isNotPrime[MAXN + 1];
int mu[MAXN + 1], phi[MAXN + 1], primes[MAXN + 1], cnt;
inline void euler() {
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    phi[1] = 1;
    for (int i = 2; i <= MAXN; i++) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= cnt; j++) {
            int t = i * primes[j];
            if (t > MAXN) break;
            isNotPrime[t] = true;
            if (i % primes[j] == 0) {
                mu[t] = 0;
                phi[t] = phi[i] * primes[j];
                break;
            } else {
                mu[t] = -mu[i];
                phi[t] = phi[i] * (primes[j] - 1);
            }
        }
    }
    for(int i=1;i<=MAXN;i++) mu[i]+=mu[i-1];
}

int solve(int n,int m,int k) {
    n/=k; m/=k;
    if(n==0 || m==0) return 0;
    int ans=0,l=1,r=0;
    if(n>m) swap(n,m);
    while(l<=n) {
        r=min(n/(n/l),m/(m/l));
        ans+=(mu[r]-mu[l-1])*(n/l)*(m/l);
        l=r+1;
    }
    return ans;
}

signed main() {
    euler();
    int t,a,b,c,d,k;
    ios::sync_with_stdio(false);
    cin>>t;
    while(t--) {
        cin>>a>>b>>c>>d>>k; --a; --c;
        cout<<solve(b,d,k)-solve(a,d,k)-solve(b,c,k)+solve(a,c,k)<<endl;
    }
}