「模拟赛20181010」太阳神 莫比乌斯反演

Posted 2021-01-14 modeststarlight

题目描述

太阳神拉很喜欢最小公倍数，有一天他想到了一个关于最小公倍数的题目。
求满足如下条件的数对((a,b))对数:(a,b)均为正整数且(a,bleq n) 而(lcm(a,b)>n)。其中的(lcm)当然表示最小公倍数。答案对(1000000007)取模。

输入

第一行一个正整数(n)。

输出

一行一个整数表示答案，对(1000000007)取模。

样例

样例输入

3

样例输出

2

数据范围

(nleq 10^{10})

题解

神题，难度非常(NOIp(rofessional))。
首先考虑到大于不好求，显然小于等于更好求，我们用补集法求出(lcm)小于等于(n)的，再与总方案数作差即可。
[sum_{i=1}^nsum_{j=1}^n lcm(i,j)leq n]
[sum_{i=1}^nsum_{j=1}^n frac{icdot j}{gcd(i,j)}leq n]
现在按套路，换枚举(d=gcd(i,j))试试
[sum_{d=1}^nsum_{i=1}^{n}sum_{j=1}^{n}[gcd(i,j)=d]&[frac{icdot j}{d}leq n]]
[sum_{d=1}^nsum_{i'=1}^{frac{n}{d}}sum_{j'=1}^{frac{n}{i'cdot d}}[gcd(i',j')=1]]
莫比乌斯函数可以出场了！
[sum_{d=1}^nsum_{i'=1}^{frac{n}{d}}sum_{j'=1}^{frac{n}{i'cdot d}}mu(gcd(i',j'))]
再换枚举(d'=gcd(i',j'))。
[sum_{d=1}^nsum_{d'=1}^{sqrt{frac{n}{d}}}mu(d') imes[i''cdot j''leq frac{n}{d'^2cdot d}的组数]]
嗯……我们把(d')放在最外边。因为(d imes d'^2leq n)，所以可以发现(d')不超过(sqrt{n})。
[sum_{d'=1}^{sqrt{n}}mu(d') imes[dcdot i''cdot j''leq frac{n}{d'^2}的组数]]
好了，我们只要求出三元组数量即可！可以令(d'leq i'' leq j'')，然后排列组合计算……

而如果规定了大小顺序，可以节省时间。比如将第一个枚数举到三次根即停止，第二个数枚举到剩余的平方根即停止，最后一个数可以直接计算。

这样的复杂度可以证明是(O(n^{frac{2}{3}}))的，同时，杜教筛也可以达到同样的复杂度，但是会有更大的常数。
(Code:)

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define N 100005
#define ll long long
#define mod 1000000007
template<typename Mytype>void Read(Mytype &p)
{
    p = 0;
    char c = getchar();
    for (; c < '0' || c > '9'; c = getchar());
    for (; c >= '0' && c <= '9'; c = getchar())p = p * 10 + c - '0';
}
ll n;
int mu[N], h[N], p[N], cnt;
int main()
{
    Read(n);
    mu[1] = 1;
    for (int i = 2; i <= N - 5; i++)
    {
        if (!h[i])
        {
            p[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt; j++)
        {
            if (i * p[j] > N - 5)
                break;
            h[i * p[j]] = 1;
            if (i % p[j] == 0)
            {
                mu[i * p[j]] = 0;
                break;
            }
            mu[i * p[j]] = -mu[i];
        }
    }
    int m = int(sqrt(n) + 0.5), ans = 0;
    for (int i = 1; i <= m; i++)
    {
        ll w = n / i / i, sum = 0;
        for (ll a = 1; a * a * a <= w; a++)
        {
            for (ll b = a; b * b <= w / a; b++)
            {
                ll c = w / a / b - b + 1;
                if (a == b)
                    sum = (sum + 1 + (c - 1) * 3) % mod;
                else
                    sum = (sum + 3 + (c - 1) * 6) % mod;
            }
        }
        sum %= mod;
        ans = (ans + mu[i] * sum) % mod;
    }
    ans = -ans;
    ans = (ans + (n % mod) * (n % mod)) % mod;
    if (ans < 0)ans += mod;
    printf("%d
", ans);
}