PASSWORD4.5--Calc

题目链接
题意简述：
求\(\sum_{s\in[1,n]\bigwedge s\notin p_i}\sum_{j=1}^{m}sj(s,j)\;mod\;10^{9} + 7\)（\((s,j)\)表示\(s,j\)的最大公约数）
标签：欧拉反演？杜教筛，数论分块
题解：
考虑\(k=0\)，要求式子为\(\sum_{i=1}^n\sum_{j=1}^{m}ij(i,j)\;mod\;10^{9} + 7\)
显然莫比乌斯反演（欧拉反演），于是开始推式子：（不妨设\(n<m\)）

\[\sum_{i=1}^n\sum_{j=1}^{m}ij(i,j) \\ = \sum_{i=1}^n\sum_{j=1}^{m}ij\sum_{d|(i,j)}\varphi(d) \\ = \sum_{d=1}^n\varphi(d)\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}ijd^{2} \\ = \sum_{d=1}^nd^{2}\varphi(d)\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}i\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}j \]

可以先线性筛出\(f(d) = d^{2}\varphi(d)\)前缀和，再数论分块求解，复杂度:\(O(\sqrt{n})\)
若\(k>0\)，可将\(p_i\)看成断点，将中间的每一个区间用前缀和求解
由于出题人较为duliu，\(n<=10^9\)，无法线性筛，于是考虑配合线性筛\(O(n^{\frac{2}{3}})\)复杂度的杜教筛
目标：\(\sum_{i=1}^nf(i)\)
过程：

\[h=g*f \\ let\;g(n)=n^{2},h(n) = \sum_{d|n}f(d)g(\frac{n}{d}) = \sum_{d|n}n^{2}\varphi(d) = n^3,S(n) = \sum_{i=1}^nf(i) \\ \sum_{i=1}^n{h_i} = \sum_{i=1}^n\sum_{d|i}g(d)f(\frac{i}{d}) = \sum_{d=1}^ng(d)\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}f(i) = \sum_{d=1}^ng(d)S(\left \lfloor \frac{n}{d} \right \rfloor) \\ \sum_{i=1}^n{h_i}=g(1)S(n) + \sum_{d=2}^ng(d)S(\left \lfloor \frac{n}{d} \right \rfloor) \\ \sum_{i=1}^n{i^3} = S(n) + \sum_{d=2}^nd^{2}S(\left \lfloor \frac{n}{d} \right \rfloor) \]

令\(t=10^9\)，线性筛前\(t^{\frac{2}{3}}\thickapprox 10^6\)个，后用杜教筛，复杂度：\(O(n^{\frac{2}{3}})\)
代码：

#include <iostream>
#include <cstdio>
#include <unordered_map>

using namespace std;

typedef long long ll;
const int N = 1e6 + 5;
const ll MOD = 1e9 + 7;
const ll inv6 = 166666668;
const ll inv2 = 500000004;
const int K = 20;

unordered_map<ll, int> f;
ll n, m, sf[N + 5], Ans;
int tot, prime[N + 5], phi[N + 5], d[K], k;
bool v[N + 5];

inline ll mul(ll a, ll b)
{
	return a * b % MOD;
}

inline void Prime(int n)
{
	v[1] = phi[1] = 1;
	for (int i = 2; i <= n; ++i)
	{
		if (!v[i])
		{
			prime[++tot] = i;
			phi[i] = i - 1;
		}
		for (int j = 1; j <= tot; ++j)
		{
			if (prime[j] > n / i)
				break ;
			v[i * prime[j]] = 1;
			if (i % prime[j])
				phi[i * prime[j]] = phi[i] * phi[prime[j]];
			else
				phi[i * prime[j]] = phi[i] * prime[j];
		}
	}
	for (int i = 1; i <= n; ++i)
	{
		sf[i] = mul(mul(i, i), phi[i]);
		sf[i] = (sf[i] + sf[i - 1]) % MOD;
	}
}

inline ll squaresum(ll l, ll r)
{
	ll s = (mul(r, mul(r + 1, 2 * r + 1)) - mul(l - 1, mul(l, 2 * l - 1)) + MOD) % MOD;
	return (s * inv6) % MOD;
}

inline ll powsum(ll x)
{
	ll s = mul(x, mul(x + 1, inv2));
	return (s * s) % MOD;
}

inline ll s_f(int n)
{
	if (n <= N)
		return sf[n];
	if (f[n])
		return f[n];
	ll ret = 0;
	for (int l = 2, r; l <= n; l = r + 1)
	{
		r = min(n / (n / l), n);
		ret -= squaresum(l, r) * s_f(n / l);
		ret = (ret + MOD) % MOD;
	}
	ret += powsum(n);
	ret = (ret + MOD) % MOD;
	return f[n] = ret;
}

inline ll calc(ll n, ll m)
{
	ll s = min(n, m), ans = 0;
	for (ll l = 1, r; l <= s; l = r + 1)
	{
		r = min(n / (n / l), m / (m / l));
		ll s1 = mul(mul(1 + (n / l), n / l), inv2);
		ll s2 = mul(mul(1 + (m / l), m / l), inv2);
		s2 = mul(s1, s2);
		ans += mul((s_f(r) - s_f(l - 1) + MOD) % MOD, s2);
		ans = (ans + MOD) % MOD;
	}
	return ans;
}

int main()
{
	ios::sync_with_stdio(false);
//	freopen("calc.in", "r", stdin);
//	freopen("calc.out", "w", stdout);
	cin >> n >> m >> k;
	for (int i = 1; i <= k; ++i)
		cin >> d[i];
	Prime(min(max(n, m), N * 1ll));
	if (!k)
	{
		cout << calc(n, m) << endl;
		return 0;
	}
	d[0] = 0;
	for (int i = 1; i <= k; ++i)
	{
		Ans += (calc(d[i] - 1, m) - calc(d[i - 1], m) + MOD) % MOD;
		Ans = (Ans + MOD) % MOD;
	}
	if (d[k] != n)
	{
		Ans += (calc(n, m) - calc(d[k], m) + MOD) % MOD;
		Ans = (Ans + MOD) % MOD;
	}
	cout << Ans << endl;
	return 0;
}