luogu4240 毒瘤之神的考验(毒瘤乌斯反演)

link

题意：求出(sum_{i=1}^nsum_{j=1}^mvarphi(ij))，对998244353取模

多组数据，(Tle 10^4,n,mle 10^5)。

前置知识：(varphi(ij)=frac{varphi(i)varphi(j)gcd(i,j)}{varphi(gcd(i,j))})

证明：我是口胡呢还是好好证呢还是口胡吧

按照欧拉函数的计算式展开，会发现，左边是(ijprod_{p|i mathrm{color{green}{or}}p|j}frac{p-1}p)

右边是(frac{iprod_{p|i}frac{p-1}pjprod_{p|j}frac{p-1}pgcd(i,j)}{gcd(i,j)prod_{p|imathrm{color{green}{and}}p|j}frac{p-1}p})

显然，根据容斥原理，两边是相等的

然后推式子

(sum_{i=1}^nsum_{j=1}^mvarphi(ij))

(=sum_{i=1}^nsum_{j=1}^mfrac{varphi(i)varphi(j)gcd(i,j)}{varphi(gcd(ij))})

(=sum_{p=1}^nfrac p{varphi(p)}sum_{i=1}^nsum_{j=1}^mvarphi(i)varphi(j)[gcd(i,j)=p])

(=sum_{p=1}^nfrac p{varphi(p)}sum_{i=1}^{n/p}sum_{j=1}^{m/p}varphi(ip)varphi(jp)[gcd(i,j)=1])

(=sum_{p=1}^nfrac p{varphi(p)}sum_{i=1}^{n/p}sum_{j=1}^{m/p}varphi(ip)varphi(jp)sum_{d|i,d|j}mu(d))

(=sum_{p=1}^nfrac p{varphi(p)}sum_{d=1}^nmu(d)sum_{i=1}^{n/dp}sum_{j=1}^{m/dp}varphi(idp)varphi(jdp))

(=sum_{q=1}^nsum_{p|q}frac{pmu(frac qp)}{varphi(p)}sum_{i=1}^{n/q}sum_{j=1}^{m/q}varphi(iq)varphi(jq))

(=sum_{q=1}^nleft(sum_{p|q}frac{pmu(frac qp)}{varphi(p)} ight)left(sum_{i=1}^{n/q}varphi(iq) ight)left(sum_{i=1}^{m/q}varphi(iq) ight))

前面这一部分好处理--(O(nlog n))枚举倍数。后面？按照套路？数论分块？怎么分？？？？

观察了你谷题解后，终于懂了

设(sum(q)=sum_{p|q}frac{pmu(frac qp)}{varphi(p)})，显然可以在(O(nlog n))的时间复杂度内处理出来。

设(g(x,y)=sum_{i=1}^xvarphi(iy))，显然有递推式(g(x,y)=g(x-1,y)+varphi(xy))。

由于(xy<=n)，对于每个(x)，有(frac nx)的数值，我们可以通过动态申请内存，在(O(nlog n))的时间复杂度和空间复杂度内求出(g)数组。

设(T(n,a,b)=sum_{q=1}^nleft(sum_{p|q}frac{pmu(frac qp)}{varphi(p)} ight)left(sum_{i=1}^{a}varphi(iq) ight)left(sum_{i=1}^{b}varphi(iq) ight)=sum_{q=1}^nsum(q)g(a,q)g(b,q))

显然T的递推式为(T(n,a,b)=T(n-1,a,b)+sum(n)g(a,n)g(b,n))

根据数论分块那套理论，对于一个(n/q)和(m/q)相同的(q)的区间，当(n/q=a,m/q=b)时，这一区间的(ans=T(r,a,b)-T(l-1,a,b))，r和l是这一区间内的最大值和最小值

我们考虑预处理(n*B*B)范围的答案，B是我们钦定的一个数字，T数组开的空间复杂度为(O(nB^2))（实际上由于(a*n,b*nle 10^5)的限制，应该开不到(O(nB^2))。

对于每次询问，我们只能在(n/qle B)时候进行数论分块操作通过(T)数组计算答案，复杂度根据数论分块那套理论为(O(sqrt n))。

对于(n/q>B)的部分，有(q<n/B)，暴力枚举(q)，通过(g)数组计算答案，这一部分单次计算的复杂度为(O(n/B))。

总复杂度为(O(nlog n+nB^2+T(sqrt n+n/B)))。实测，B开到50左右跑的快一点，且内存占用超小。

下面是乱七八糟的代码= =

注意讲文明，new来的内存要主动回收垃圾

注意取模（这题如果写的复杂度没错的话不卡常，开#define int long long也是没问题的

#include <cstdio>
#include <functional>
using namespace std;

const int p = 998244353;
const int b = 50;
bool vis[100010];
int prime[100010], tot, fuck = 100000;
int mu[100010], phi[100010], invphi[100010];
int sum[100010];
int *g[100010], *t[100][100]; //注意这里t数组下标是[2][3][1]

int qpow(int x, int y)
{
	int res = 1;
	for (x %= p; y > 0; x = x * (long long)x % p, y >>= 1) if (y & 1) res = res * (long long)x % p;
	return res;
}

int main()
{
	//线性筛phi,mu，预处理前面的部分
	phi[1] = mu[1] = invphi[1] = 1;
	for (int i = 2; i <= fuck; i++)
	{
		if (vis[i] == false) prime[++tot] = i, phi[i] = i + (mu[i] = -1);
		for (int j = 1; j <= tot && i * prime[j] <= fuck; j++)
		{
			vis[i * prime[j]] = true;
			if (i % prime[j] == 0) { phi[i * prime[j]] = phi[i] * prime[j]; break; }
			phi[i * prime[j]] = phi[i] * (prime[j] - 1);
			mu[i * prime[j]] = -mu[i];
		}
		invphi[i] = qpow(phi[i], p - 2);
		if (phi[i] * (long long)invphi[i] % p != 1) { return -233; printf("cnm
"); }
	}
	for (int pp = 1; pp <= fuck; pp++)
		for (int q = pp, d = 1; q <= fuck; q += pp, d++)
			sum[q] = (sum[q] + pp * (long long)invphi[pp] % p * mu[d]) % p, sum[q] += (sum[q] < 0 ? p : 0);

	//处理g数组
	for (int i = 1; i <= fuck; i++)
	{
		g[i] = new int[(fuck / i) + 1], g[i][0] = 0;
		for (int j = 1, sb = fuck / i; j <= sb; j++)
			g[i][j] = (g[i][j - 1] + phi[i * j]) % p;
	}
	
	//处理t数组 注意有第一维<=第二维，因为下面我们强制n<=m了
	for (int j = 1; j <= b; j++)
		for (int k = j; k <= b; k++)
		{
			int len = fuck / max(j, k);
			t[j][k] = new int[len + 1], t[j][k][0] = 0;
			for (int i = 1; i <= len; i++)
				t[j][k][i] = (t[j][k][i - 1] + sum[i] * (long long)g[i][j] % p * g[i][k] % p) % p;
		}
	
	//处理询问
	int tat;
	scanf("%d", &tat);
	while (tat --> 0)
	{
		int n, m, res = 0;
		scanf("%d%d", &n, &m);
		if (n > m) swap(n, m);
		//对于n/q>b的部分，暴力，通过g数组和sum数组计算计算
		for (int i = 1, sb = m / b; i <= sb; i++)
			res = (res + sum[i] * (long long)g[i][n / i] % p * g[i][m / i] % p) % p;
		//对于n/q<b的部分，数论分块，通过b数组计算
		for (int i = m / b + 1, j; i <= n; i = j + 1)
		{
			j = min(n / (n / i), m / (m / i));
			res = (res + t[n / i][m / i][j] - t[n / i][m / i][i - 1]) % p, res += (res < 0 ? p : 0);
		}
		printf("%d
", res);
	}
	
	//垃圾回收
	for (int i = 1; i <= fuck; i++)
		delete []g[i], g[i] = 0;
	for (int i = 1; i <= b; i++)
		for (int j = i; j <= b; j++)
			delete[] t[i][j], t[i][j] = 0;
	return 0;
}