Description：

求 $\begin{align}\sum_{i=1}^n\sum_{j=1}^n \max(i,j)\times \sigma(ij)\end{align}​$

$1\le n \le 10^7,1\le t \le 50000$

Solution：

原式 $\begin{align}=\sum_{i=1}^n\sum_{j=1}^ii\sigma_1(ij)-\sum_{i=1}^ni\sigma_1(i^2)\end{align}$ ，后面那部分可以线性筛。

$\sigma_1(i^2)=1+i+i^2(i\ is\ prime)\\\sigma_1((i\times prime[j])^2)=\sigma_1(i^2)\times \sigma_1(prime[j]^2)(i\%j\ne0)$

当 $i\%prime[j]=0$ 时，由于 $\sigma_1((p^k\times q\times p)^2)=\sigma_1((p^k\times p)^2)\times \sigma_1(q^2)=(\sigma_1((p^k)^2)+p^{2k+1}+p^{2k+2})\times \sigma_1(q^2)$ ，所以 $\sigma_1((i\times prime[j])^2)=(\sigma_1((prime[j]^{h[i]})^2)+prime[j]^{2h[i]+1}+prime[j]^{2h[i]+2})\times \sigma_1((\frac i {prime[j]^{h[i]}})^2)$

然后再考虑前半部分，有：

$\begin{align}原式=\sum_{i=1}^ni\sum_{j=1}^i\sigma_1(ij)=\sum_{i=1}^ni\sum_{j=1}^i\sum_{x|i}\sum_{y|j}x\frac j y[gcd(x,y)=1]=\sum_{i=1}^ni\sum_{j=1}^i\sum_{x|i}\sum_{y|j}x\frac j y\sum_{k|x,k|y}\mu(k)\end{align}$

把 $k$ 放到开始：

$\begin{align}原式&=\sum_{k=1}^n\mu(k)\sum_{i=1}^{n}i\sum_{j=1}^i\sum_{k|x,x|i}x\sum_{k|y,y|j}\frac j y\\&=\sum_{k=1}^n\mu(k)\sum_{i=1}^{\lfloor \frac n k\rfloor}ik\sum_{j=1}^i\sum_{x|i}xk\sum_{y|j}\frac j y\\&=\sum_{k=1}^n\mu(k)k^2\sum_{i=1}^{\lfloor \frac n k\rfloor}i\sum_{x|i}x\sum_{j=1}^i\sum_{y|j}\frac j y\end{align}$

设 $\begin{align}g(n)=n\sum_{k|n}k\sum_{i=1}^n\sum_{y|i}\frac i y=n\times\sigma_1(n)\times \sum_{i=1}^n\sigma_1(i)\end{align}$

$g$ 可以通过预处理 $\sigma_1$ 前缀和求。

于是：

$\begin{align}原式&=\sum_{k=1}^n\mu(k)k^2\sum_{i=1}^{\lfloor\frac n k\rfloor}g(i)\end{align}$

设 $d=i\times k$ ，于是有 $\begin{align}原式&=\sum_{d=1}^n\sum_{k|d}\mu(k)k^2g(\frac d k)\end{align}$

不妨设 $\begin{align}f(d)=\sum_{k|d}\mu(k)k^2g(\frac d k)\end{align}$ ，求出 $f$ 后就可以 $O(n\ln n)$ 求原式了。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

#define MAXN 1000001

typedef long long ll;

ll prime[MAXN],tot = 0;

bool isprime[MAXN];

ll sigma[MAXN],sigmas[MAXN],sigma2[MAXN];

ll divpow[MAXN];

ll mu[MAXN];

ll f[MAXN],g[MAXN],h[MAXN];

#define MOD 1000000007

inline void init()

{

	for(register ll i = 2;i < MAXN;++i)isprime[i] = true;

	divpow[1] = 1;sigma2[1] = 1;sigma[1] = 1;sigmas[1] = 1;mu[1] = 1;

	for(register ll i = 2;i < MAXN;++i)

	{

		if(isprime[i])

		{

			prime[++tot] = i;

			sigma2[i] = (1 + i + i * i) % MOD;

			divpow[i] = i;

			sigma[i] = 1 + i;

			mu[i] = -1;

		}

		for(register ll j = 1;j <= tot && i * prime[j] < MAXN;++j)

		{

			register ll k = i * prime[j];

			isprime[k] = false;

			if(i % prime[j] != 0)

			{

				divpow[k] = prime[j];

				sigma2[k] = sigma2[i] * sigma2[prime[j]] % MOD;

				sigma[k] = sigma[i] * sigma[prime[j]] % MOD;

				mu[k] = mu[i] * -1;

			}

			else

			{

				divpow[k] = divpow[i] * prime[j];

				sigma2[k] = (sigma2[divpow[i]] + divpow[i] * divpow[i] % MOD * prime[j] + divpow[i] * divpow[i] % MOD * prime[j] % MOD * prime[j]) * sigma2[i / divpow[i]] % MOD;

				sigma[k] = (sigma[divpow[i]] + divpow[i] * prime[j]) % MOD * sigma[i / divpow[i]] % MOD;

				mu[k] = 0;

				break;

			}

		}

		sigmas[i] = (sigmas[i - 1] + sigma[i]) % MOD;

	}

	for(register ll i = 1;i < MAXN;++i)g[i] = (i * sigma[i] % MOD * sigmas[i]) % MOD;

	for(register ll i = 1;i < MAXN;++i)

	{

		for(register ll j = i,k = 1;j < MAXN;j += i,++k)

			f[j] = (f[j] + mu[i] * ((i * i % MOD * g[k]) % MOD) + MOD) % MOD;

		f[i] = (f[i] + f[i - 1]) % MOD;

		h[i] = (h[i - 1] + i * sigma2[i]) % MOD;

	}

	return;

}

inline int read()

{

	register int 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;

}

int main()

{

	init();

	register int testcases = 0;

	scanf("%d",&testcases);

	register int k;

	for(register int i = 1;i <= testcases;++i)

	{

		k = read();

		printf("Case #%d: %lld\n",i,(2 * f[k] - h[k] + MOD) % MOD);

	}

	return 0;

}