题目描述
太阳神拉很喜欢最小公倍数，有一天他想到了一个关于最小公倍数的题目。
求满足如下条件的数对((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);
}