「模拟赛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);
}