Description：

给出六个数组 $A,B,C,D,E,F$ ，求： $$ \sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^nA_iB_jC_kD_{\gcd(i,j)}E_{\gcd(j,k)}F_{\gcd(k,i)} $$ $1\leqslant n\leqslant 10^5$

Solution：

首先发现最不好处理的是 $\gcd$ ，于是考虑把 $\gcd$ 拆开，由莫比乌斯反演有一个定理： $$ f(n)=\sum_{d|n}\sum_{d'|d}\mu(d')f(\frac d{d'}) $$ 我们用 $\times$ 来表示迪利克雷卷积，那么有： $$ f(n)=\sum_{d|n}(\mu\times f)(d) $$ 那么我们就可以把 $\gcd$ 反演掉： $$ \begin{align} &\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^nA_iB_jC_kD_{\gcd(i,j)}E_{\gcd(j,k)}F_{\gcd(k,i)}\\ =&\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^nA_iB_jC_kE_{\gcd(k,j)}\sum_{p|\gcd(i,j)}(\mu\times D)(p)\sum_{q|\gcd(i,k)}(\mu\times F)(q)\\ =&\sum_{p=1}^n\sum_{q=1}^n(\mu\times D)(p)(\mu\times F)(q)\sum_{p|i,q|i}^nA_i\sum_{p|j}^n\sum_{q|k}^nB_jC_kE_{\gcd(k,j)}\\ d&=\gcd(p,q)\\ =&\sum_{d=1}^n\sum_{p=1}^{\lfloor\frac nd\rfloor}(\mu\times D)(pd)\sum_{q=1}^{\lfloor\frac nd\rfloor}(\mu\times F)(qd)[\gcd(p,q)=1]\sum_{pqd|i}A_i\sum_{pd|j}^n\sum_{qd|k}^nB_jC_kE_{\gcd(k,j)}\\ \end{align} $$ 考虑最后面那一个求和式： $$ \begin{align} &\sum_{pq|i}^{\lfloor\frac nd\rfloor}A(id)\sum_{p|j}^{\lfloor\frac nd\rfloor}\sum_{q|k}^{\lfloor\frac nd\rfloor}B(jd)C(kd)E(\gcd(k,j)d)\\ =&\sum_{pq|i}^{\lfloor\frac nd\rfloor}A(id)\sum_{p|j}^{\lfloor\frac nd\rfloor}\sum_{q|k}^{\lfloor\frac nd\rfloor}B(jd)C(kd)E_d(\gcd(k,j))\\ =&\sum_{pqd|i}^nA(i)\sum_{p|j}^{\lfloor\frac nd\rfloor}\sum_{q|k}^{\lfloor\frac nd\rfloor}B(jd)C(kd)\sum_{d'|k,d'|j}(\mu\times E_d)(d')\\ &A'(d)=\sum_{d|i}^nA(i)\\ =&A'(pqd)\sum_{p|j}^{\lfloor\frac nd\rfloor}B_d(j)\sum_{q|k}^{\lfloor\frac nd\rfloor}C_d(k)\sum_{d'|k,d'|j}(\mu\times E_d)(d')\\ \end{align} $$ 带回原式： $$ \begin{align} =&\sum_{d=1}^n\sum_{p=1}^{\lfloor\frac nd\rfloor}(\mu\times D)(pd)\sum_{q=1}^{\lfloor\frac nd\rfloor}(\mu\times F)(qd)[\gcd(p,q)=1]A'(pqd)\sum_{p|j}^{\lfloor\frac nd\rfloor}B_d(j)\sum_{q|k}^{\lfloor\frac nd\rfloor}C_d(k)\sum_{d'|k,d'|j}(\mu\times E_d)(d')\\ \end{align} $$

因为要求 $p\times q\times d\leqslant n$ ，也就是说 $p\times q\leqslant \lfloor\frac nd\rfloor$ ，所以 $q$ 那一个循环是跑不满的。

然后这里用了一个很神奇的操作，就是观察上面那个式子，因为 $p\times q\leqslant \lfloor\frac nd\rfloor$ ，所以我们可以分 $p\geqslant q$ 和 $p<q$ 两种情况分别计算，这样做的好处是总有一个变量 $\leqslant \sqrt{\lfloor\frac nd\rfloor}$ 因此可以暴力枚举它，里面还能再套一个 $\lfloor\frac nd\rfloor$ 。

以 $p<q$ 为例，考虑枚举 $p$ ： $$ \begin{align} =&\sum_{d=1}^n\sum_{p=1}^{\lfloor\frac nd\rfloor}(\mu\times D)(pd)\sum_{q=1}^{\lfloor\frac nd\rfloor}(\mu\times F)(qd)[\gcd(p,q)=1]A'(pqd)\sum_{p|j}^{\lfloor\frac nd\rfloor}B_d(j)\sum_{q|k}^{\lfloor\frac nd\rfloor}C_d(d)\sum_{d'|k,d'|j}(\mu\times E_d)(d')\\ \end{align} $$ 最后面那个 $\sum$ 非常不好处理因为和两个循环变量都有关，我们想办法搞掉它。

也就是说我们希望求出来一个 $S(p,q,d)$ 表示： $$ S(p,q,d)=\sum_{p|j}^{\lfloor\frac nd\rfloor}B_d(j)\sum_{q|k}^{\lfloor\frac nd\rfloor}C_d(d)\sum_{d'|k,d'|j}(\mu\times E_d)(d')\\ $$ $S(p,q,d)$ 需要 $O(\lfloor\frac nd\rfloor)$ 或类似的预处理 $q\in[1,\lfloor\frac nd\rfloor]$ 时的值，并能够 $O(1)$ 查询。 $$ S(p,q,d)=\sum_{d'=1}^{\lfloor\frac nd\rfloor}(\mu\times E_d)(d')\sum_{p|j,d'|j}^{\lfloor\frac nd\rfloor}B_d(j)\sum_{q|k,d'|k}^{\lfloor\frac nd\rfloor}C_d(k) $$ 但是这样还是很难处理，因为后两重和式每个都有两个限制，考虑简化。

因为我们暴力枚举了 $p​$ ，因此可以把与 $q​$ 有关的放到最外层循环，这样内层与 $q​$ 无关，可以在 $p​$ 的每次循环预处理，即： $$ \begin{align} S(p,q,d)&=\sum_{d'=1}^{\lfloor\frac nd\rfloor}(\mu\times E_d)(d')\sum_{p|j,d'|j}^{\lfloor\frac nd\rfloor}B_d(j)\sum_{q|k,d'|k}^{\lfloor\frac nd\rfloor}C_d(k)\\ &=\sum_{q|k}^{\lfloor\frac nd\rfloor}C_d(k)\sum_{d'|k}^{\lfloor\frac nd\rfloor}(\mu\times E_d)(d')\sum_{p|j,d'|j}^{\lfloor\frac nd\rfloor}B_d(j)\\ \end{align} $$ 然后就可以做了。

复杂度嘛，因为里面实际上枚举了两层 $\lfloor\frac nd\rfloor$ ，但是 $p\times q\leqslant \lfloor\frac nd\rfloor$ 所以跑不满，听说复杂度是 $O(n^{1.5}\log\log n)$ 的。

从这道题学到的东西：

1、一定不要枚举上界宽松的变量，而是枚举上界紧的变量，并用他们去限制上界松的变量。

2、 $f(n)=\sum_{d|n}(\mu\times f)(d)$ ，常用来拆开函数中嵌套的 $\gcd$ 。

3、如果出现 $a\times b\leqslant n$ 的情况要想到分开处理。

4、一个神奇的计算迪利克雷卷积的方式。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n;

#define MAXN 100010

typedef unsigned long long ull;

ull A[MAXN],B[MAXN],C[MAXN],D[MAXN],E[MAXN],F[MAXN];

inline ull rd()

{

	register ull res = 0;register char c = getchar();

	while(!isdigit(c))c = getchar();

	while(isdigit(c))res = (res << 1) + (res << 3) + c - '0',c = getchar();

	return res;

}

ull prime[MAXN],tot = 0;

bool isprime[MAXN];

ull mu[MAXN];

inline int gcd(int a,int b){return (b == 0 ? a : gcd(b,a % b));}

ull Ed[MAXN];

ull s1[MAXN],s2[MAXN];

#define R register

int main()

{freopen("sum25.in","r",stdin);

	scanf("%d",&n);

	for(R int i = 1;i <= n;++i)A[i] = rd();

	for(R int i = 1;i <= n;++i)B[i] = rd();

	for(R int i = 1;i <= n;++i)C[i] = rd();

	for(R int i = 1;i <= n;++i)D[i] = rd();

	for(R int i = 1;i <= n;++i)F[i] = rd();

	for(R int i = 1;i <= n;++i)E[i] = rd();

	for(R int i = 2;i <= n;++i)isprime[i] = true;

	mu[1] = 1;

	for(R int i = 2;i <= n;++i)

	{

		if(isprime[i])prime[++tot] = i,mu[i] = -1;

		for(R int j = 1;j <= tot && i * prime[j] <= n;++j)

		{

			R int k = i * prime[j];

			isprime[k] = false;

			if(i % prime[j] == 0){mu[k] = 0;break;}

			else mu[k] = -mu[i];

		}

	}

	for(R int i = 1;i <= tot && prime[i] <= n;++i)for(R int j = n / prime[i],k = j * prime[i];j >= 1;--j,k -= prime[i])A[j] += A[k];

	R ull ans = 0;

	for(R int i = 1;i <= tot && prime[i] <= n;++i)for(R int j = n / prime[i],k = j * prime[i];j >= 1;--j,k -= prime[i])F[k] -= F[j],D[k] -= D[j];

	for(R int d = 1;d <= n;++d)

	{

		for(R int i = 1;i <= n / d;++i)Ed[i] = E[i * d];

		for(R int i = 1;i <= tot && prime[i] <= n / d;++i)for(R int j = n / d / prime[i],k = j * prime[i];j >= 1;--j,k -= prime[i])Ed[k] -= Ed[j];

		for(R int p = 1;p * p <= n / d;++p)

		{

			for(R int i = 1;i <= n / d;++i)s1[i] = (i % p == 0) * B[i * d],s2[i] = (i % p == 0) * C[i * d];

			for(R int i = 1;i <= tot && prime[i] <= n / d;++i)for(R int j = n / d / prime[i],k = j * prime[i];j >= 1;--j,k -= prime[i])s1[j] += s1[k],s2[j] += s2[k];

			for(R int i = 1;i <= n / d;++i)s1[i] *= Ed[i],s2[i] *= Ed[i];

			for(R int i = 1;i <= tot && prime[i] <= n / d;++i)for(R int j = prime[i],k = 1;j <= n / d;++k,j += prime[i])s1[j] += s1[k],s2[j] += s2[k];

			for(R int i = 1;i <= n / d;++i)s1[i] *= C[i * d],s2[i] *= B[i * d];

			for(R int i = 1;i <= tot && prime[i] <= n / d;++i)for(R int j = n / d / prime[i],k = j * prime[i];j >= 1;--j,k -= prime[i])s1[j] += s1[k],s2[j] += s2[k];

			for(R int q = p + 1;p * q <= n / d;++q)if(gcd(p,q) == 1)ans += A[d * p * q] * s1[q] * F[q * d] * D[p * d] + A[d * p * q] * s2[q] * D[q * d] * F[p * d];

			if(p == 1)ans += A[d] * s2[1] * D[d] * F[d];

		}

	}

	cout << ans << endl;

	return 0;

}