最小公倍数计数

51nod 1222

求满足$a \leq lcm(i, j) \leq b$ 的二元组$[i, j]$的个数

$a, b \leq 10^{11}$

---

设

$$S(n) = \sum_{i=1}^n \sum_{j=1}^i [\dfrac{ij}{gcd(i, j)} \leq n ]$$

则

$$ans = S(b) - S(a-1)$$

$$S(n) = \sum_{d=1}^n \sum_{i=1}^n \sum_{j=1}^n [ij \leq dn] [gcd(i, j) = d]$$

$i$/=$d$, $j$/=$d$

$$S(n) = \sum_{d=1}^n \sum_{i=1}^{n/d} \sum_{j=1}^{n/d} [ idjd \leq dn][gcd(id, jd) = d]$$

$$S(n) = \sum_{d=1}^n \sum_{i=1}^{n/d} \sum_{j=1}^{n/d} [ ij \leq \dfrac{n}{d}][gcd(i, j) = 1]$$

$$S(n) = \sum_{d=1}^n \sum_{i=1}^{n/d} \sum_{j=1}^{n/d} [ ij \leq \dfrac{n}{d}] \sum_{k|i 且 k|j} \mu(k)$$

枚举一下$k$

$$S(n) = \sum_{d=1}^n \sum_{k=1}^{n/d} \mu(k) \sum_{i=1}^{n/dk} \sum_{j=1}^{n/dk} [ikjk \leq \dfrac{n}{d}]$$

$$S(n) = \sum_{k=1}^n \mu(k) \sum_{d=1}^{n/k} \sum_{i=1}^{n/dk} \sum_{j=1}^{n/dk} [ij \leq \dfrac{n}{dk^2}]$$

我们发现如果$k > \sqrt{n}, ij$就必须小于$0$了

$$S(n) = \sum_{k=1}^{\sqrt{n}} \mu(k) \sum_{d=1}^{n/k} \sum_{i=1}^{n/dk} \sum_{j=1}^{n/dk} [ijd \leq \dfrac{n}{k^2}]$$

此时的$i, j$代表$idk, jdk$
我们发现后面的问题就是求出有多少满足$abc \leq n$的三元组$[a, b, c]$
我们设 $a \leq b \leq c$，这样就可以枚举出有序的所有三元组。
$if( a \ne b \ne c), ans += 3$
$if( a = b \ne c ), ans += 2$
$if( a = b = c ), ans += 1$

code：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2000000;

bool is_prime[N]; 
int tot = 0, prime[N], mobius[N]; 
ll l, r, ans = 0;

inline void calc_mobius( int size )
{
    memset( is_prime, true, sizeof( is_prime ) );
    is_prime[1] = false; mobius[1] = 1;

    for( int i = 2; i <= size; i ++ )
    {
        if( is_prime[i] )
            prime[++tot] = i, mobius[i] = -1;

        for( int j = 1; j <= tot; j ++ )
        {
            int p = prime[j];
            if( i * p > size ) break;
            is_prime[i*p] = false;

            if( i % p == 0 )
            {
                mobius[i*p] = 0;
                break;
            }

            mobius[i*p] = -mobius[i];
        }
    }
}

ll calc( ll n )
{
    ll ret = 0;

    for( ll k = 1; k*k <= n; k ++ )
    {
        if( !mobius[k] ) continue;

        ll lim = n / ( k * k ), sum = 0;
        for( ll d = 1; d*d*d <= lim; d ++ )
        {
            for( ll i = d+1; d*i*i <= lim; i ++ )
            {
                ll j = lim / ( i * d ); 
                if( j < i ) continue;

                sum += 2; // i == j
                if( j < i + 1 ) continue;
                sum += 3 * ( j - i ); 
            }

            ll tmp = lim / ( d * d );
            if( tmp < d ) continue;

            sum ++; // d == i == j
            if( tmp - d < 1 ) continue;
            sum += ( tmp - d ) * 2; // d == i
        }

        ret += mobius[k] * sum;
    }

    return ret;
}

int main()
{
    calc_mobius( 1200000 );

    scanf( "%lld%lld", &l, &r );

    ans = calc( r ); 
    ans -= calc( l - 1 );

    printf( "%lld\n", ans );

    return 0;
}