问题为：$$ans=\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1} \frac{\phi(i)\phi(j^2)\phi(k^3)}{\phi(i)\phi(j)\phi(k)}\phi(gcd(i,j,k))\ \ \ 1\leq n\leq10^7$$
对上式作以下推导

\begin{gather*} 设x=p_1^{k_1}*……*p_n^{k_n}\\ \therefore \frac{\phi(x^2)}{\phi(x)}=\frac{(p_1-1)*p_1^{2k_1-1}*……*(p_n-1)*p_n^{2k_n-1}}{(p_1-1)*p_1^{k_1-1}*……*(p_n-1)*p_n^{k_n-1}}=x\\ ans=\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1}jk^2\phi(gcd(i,j,k))\\ ans=\sum^{n}_{d=1}\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1}jk^2\phi(d)[gcd(i,j,k)==d]\\ ans=\sum^{n}_{d=1}\phi(d)\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1}jk^2[gcd(i,j,k)==d]\\ 设f(n,d)=\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1}jk^2[gcd(i,j,k)==d]\\ 设F(n,d)=\sum^{n}_{i=1}\sum^{n}_{j=1}\sum^{n}_{k=1}jk^2[d|gcd(i,j,k)]\\ F(n,d)=d^3\sum^{\lfloor\frac{n}{d}\rfloor}_{i=1}1\sum^{\lfloor\frac{n}{d}\rfloor}_{j=1}j\sum^{\lfloor\frac{n}{d}\rfloor}_{k=1}k^2\\ F(n,d)=d^3\frac{\lfloor\frac{n}{d}\rfloor^3(\lfloor\frac{n}{d}\rfloor+1)^2(2\lfloor\frac{n}{d}\rfloor+1)}{12}\\ 由莫比乌斯反演得f(n,d)=\sum^{\lfloor\frac{n}{d}\rfloor}_{x=1}\mu(x)F(n,dx)\\ f(n,d)=d^3\sum^{\lfloor\frac{n}{d}\rfloor}_{x=1}\mu(x)\frac{\lfloor\frac{n}{xd}\rfloor^3(\lfloor\frac{n}{xd}\rfloor+1)^2(2\lfloor\frac{n}{xd}\rfloor+1)}{12}\\ ans=\sum^{n}_{d=1}\phi(d)*d^3\sum^{\lfloor\frac{n}{d}\rfloor}_{x=1}\mu(x)\frac{\lfloor\frac{n}{xd}\rfloor^3(\lfloor\frac{n}{xd}\rfloor+1)^2(2\lfloor\frac{n}{xd}\rfloor+1)}{12}\\ 设t=dx并且带入以上表达式得\\ ans=\sum^{n}_{t=1}\frac{\lfloor\frac{n}{t}\rfloor^3(\lfloor\frac{n}{t}\rfloor+1)^2(2\lfloor\frac{n}{t}\rfloor+1)}{12}d^3\sum_{d|t}\mu(\frac{t}{d})\phi(d)\\ \end{gather*}

到现在为止关于表达式的推导已经结束了，接下来就是怎么快速的计算出答案。

显然对于表达式的前半截$\sum^{n}_{t=1}\frac{\lfloor\frac{n}{t}\rfloor^3(\lfloor\frac{n}{t}\rfloor+1)^2(2\lfloor\frac{n}{t}\rfloor+1)}{12}$可以在每次询问的时候使用整数分块来回答。这儿还有一个不好处理的点，我们没办法求在模$2^{30}$的意义下12的逆元。于是我们观察上式$\lfloor\frac{n}{t}\rfloor和\lfloor\frac{n}{t}\rfloor+1$中必定有一个能够被2整除，于是我们在整数分块的时候先判断上式哪个能够被2整除并且除于两个2。然后再取模并对于最后的答案乘上3的逆元。计算3在模$2^{30}$的意义下的逆元不能用费马小定理，只能使用拓展欧几里得算法或者是欧拉定理。$\phi(2^{30})=\phi(1073741824)=536870912$最后求出来3的逆元为715827883。

对于后半截$设g(n)=d^3\sum_{d|n}\mu(\frac{n}{d})\phi(d)$，本来对于一般的题目这儿用$O(NlogN)$的方法预处理是完全足够的，但是$O(NlogN)$对于$10^7$还是太慢了啊。于是我们我们考虑到$\mu(x)和\phi(x)$为两个积性函数。
根据积性函数的性质，我们知道积性函数的狄利克雷卷积还是积性函数。所以$g(n)$为一个积性函数。我们可以使用线性筛来推得$g(n)$，要线性筛推$g(i*j)$对于i与j互质的情况$g(i*j)=g*(i)*g(j)$，我们可以把i与j互质的情况简化为需要推导出$g(p^k)$

\begin{align*} g(p^k) &=\sum^{k}_{i=0}\phi(p^k)\mu(p^{k-i})\\ g(p^k) &=\phi(p^k)\mu(1)+\phi(p^{k-1})\mu(p)\\ g(p^k) &=\phi(p^k)-\phi(p^{k-1})\\ \end{align*} $$若k == 1 g(p) = p - 2$$ \begin{align*} 若k != 1\ g(p) &= p^k-p^{k-1}-(p^{k-1}-p^{k-2})\\ g(p) &= p^k-2p^{k-1}+p^{k-2}\\ g(p) &= (p-1)^2p^{k-2}\\ \end{align*}

这样的话就可以用线性筛递推求得g(n)

#include <bits/stdc++.h>

const int MAXN = 10000000;
const long long MOD = 1073741824;

using namespace std;

int T,n,cnt,inv;
int prime[MAXN+5],num[MAXN+5];
long long g[MAXN+5],ans;
bool vis[MAXN+5];

int pow_mod(long long a,int b){
    long long ans = 1;
    while(b){
        if(b & 1) ans = ans * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    ans %= MOD;
    return ans;
}

inline void init()
{
    g[1] = 1;
    for(int i = 2;i <= MAXN;i++){
        if(!vis[i]){
            prime[++cnt] = i;
            g[i] = i - 2;
            num[i] = 1;
        }
        for(int j = 1;j <= cnt && prime[j] * i <= MAXN;j++){
            vis[prime[j]*i] = true;
            if(i%prime[j] == 0){
                num[i*prime[j]] = num[i] + 1;
                if(num[i] == 1){
                    if(prime[j] == 2) g[i*prime[j]] = g[i/prime[j]];
                    else g[i*prime[j]] = g[i]/(prime[j]-2)*(prime[j]-1)*(prime[j]-1);
                }
                else g[i*prime[j]] = g[i]*prime[j];
                break;
            }
            else {
                g[prime[j]*i] = g[prime[j]] * g[i];
                num[prime[j]*i] = 1;
            }
        }
    }
    
    for(int i = 1;i <= MAXN;i++) g[i] = (g[i-1]+g[i]*i%MOD*i%MOD*i%MOD)%MOD; 
        
}

int main()
{
    init();
    
    inv = 715827883;
    scanf("%d",&T);
    while(T--){
        scanf("%d",&n);
        ans = 0;
        int R;
        long long m = 0,ret = 0;
        for(int i = 1;i <= n;i = R + 1)
        {
            R = n/(n/i);
            m = n/i;ret = (g[R]-g[i-1]+MOD)%MOD;
            if(m%2 == 0)
            {
                ans += ((m/2)*(m/2)%MOD*m%MOD*(m+1)%MOD*(m+1)%MOD*(m*2+1)*inv%MOD*ret%MOD)%MOD;
            }
            else {
                ans += (m*m%MOD*m%MOD*((m+1)/2)%MOD*((m+1)/2)%MOD*(m*2+1)*inv%MOD*ret%MOD)%MOD;
            }
            ans %= MOD;
            //printf("i=%d R=%d m=%lld ret=%lld ans=%lld\n",i,R,m,ret,ans);
        }
    
        printf("%lld\n",ans);
    }
    return 0;
}