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

题目描述

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

输入

第一行一个正整数$n$。

输出

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

样例

样例输入

3

样例输出

2

数据范围

$n\leq 10^{10}$

题解

神题，难度非常$NOIp(rofessional)$。
首先考虑到大于不好求，显然小于等于更好求，我们用补集法求出$lcm$小于等于$n$的，再与总方案数作差即可。

$$\sum_{i=1}^n\sum_{j=1}^n lcm(i,j)\leq n$$

$$\sum_{i=1}^n\sum_{j=1}^n \frac{i\cdot j}{gcd(i,j)}\leq n$$

现在按套路，换枚举$d=gcd(i,j)$试试

$$\sum_{d=1}^n\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=d]\&[\frac{i\cdot j}{d}\leq n]$$

$$\sum_{d=1}^n\sum_{i'=1}^{\frac{n}{d}}\sum_{j'=1}^{\frac{n}{i'\cdot d}}[gcd(i',j')=1]$$

莫比乌斯函数可以出场了！

$$\sum_{d=1}^n\sum_{i'=1}^{\frac{n}{d}}\sum_{j'=1}^{\frac{n}{i'\cdot d}}\mu(gcd(i',j'))$$

再换枚举$d'=gcd(i',j')$。

$$\sum_{d=1}^n\sum_{d'=1}^{\sqrt{\frac{n}{d}}}\mu(d')\times[i''\cdot j''\leq \frac{n}{d'^2\cdot d}的组数]$$

嗯……我们把$d'$放在最外边。因为$d\times d'^2\leq n$，所以可以发现$d'$不超过$\sqrt{n}$。

$$\sum_{d'=1}^{\sqrt{n}}\mu(d')\times[d\cdot 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\n", ans);
}