2016广东省省赛H题

题意：计算：

$$\sum_{i=1}^n\sum_{j=1}^mi^2j^2gcd(i,j)$$ （时限1000ms）

input:

T(t<=1000),接下来t行，m,n(n,m<1e6).

题目难点在于多了i,j的平方，由对称性，设n < m，可以这样化简：

$$\begin{align}&\sum_{i=1}^n\sum_{j=1}^mi^2j^2gcd(i,j)\=& \sum_{i=1}^n\sum_{j=1}^m\sum_{d|i,d|j}i^2j^2d[(\frac id,\frac jd)=1]\=& \sum_{d=1}^n\sum_{i=1}^{\lfloor \frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor}i^2j^2d^5[(i,j)=1]\tag 1\=& \sum_{d=1}^n\sum_{i=1}^{\lfloor \frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor}i^2j^2d^5\sum_{d‘|i,d‘|j}\mu(d‘)\=& \sum_{d=1}^n\sum_{d‘=1}^{\lfloor \frac nd\rfloor}\sum_{i=1}^{\lfloor \frac n{dd‘} \rfloor}\sum_{j=1}^{\lfloor \frac m{dd‘} \rfloor}i^2j^2d‘^4d^5\mu(d‘)\\stackrel{d‘=\frac kd}=&\sum_{k=1}^n\sum_{d|k}\sum_{i=1}^{\lfloor \frac nk\rfloor}\sum_{j=1}^{\lfloor \frac mk\rfloor}i^2j^2k^4d\mu(\frac k{d})\=& \sum_{k=1}^n\sum_{i=1}^{\lfloor \frac nk\rfloor}\sum_{j=1}^{\lfloor \frac mk\rfloor}i^2j^2k^4(\mu*Id) \tag 2 \=& \sum_{k=1}^n\sum_{i=1}^{\lfloor \frac nk\rfloor}\sum_{j=1}^{\lfloor \frac mk\rfloor}i^2j^2k^4\phi(k)\=& \sum_{k=1}^nk^4\phi(k)(\sum_{i=1}^{\lfloor \frac nk\rfloor}i^2)(\sum_{j=1}^{\lfloor \frac mk\rfloor}j^2)\tag 3\end{align}$$

推导到这里就可以解决了。

$k^4\phi(k)$这部分可以打表求前缀和，$\sum_{i=1}^{\lfloor \frac nk\rfloor}i^2，\sum_{j=1}^{\lfloor \frac mk\rfloor}j^2$这部分打表或者用平方和公式。对于固定的$n, \lfloor\frac nk\rfloor$共有$O(\sqrt{n})$个取值，所以总复杂度是$O(\sqrt n + \sqrt m)$

(2)式中符号*代表卷积，$Id$函数代表$Id(n)=n$;他们的卷积是$\phi(n)$,下面给出证明：
由于$\mu,Id$都是积性函数，所以他们卷积也是积性函数.
令$n=\prod_{i=1}^s p_i^{k_i}$那么：

$$\sum_{d|n}d\mu(\frac nd) =\prod_{i=1}^s(\sum_{j=0}^{k_i}p_i^j\mu(\frac {p_i^{k_i}}{p_i^j})) =\prod_{i=1}^s(p_i^{k_i}-p_i^{k_i-1})=\phi(n)$$

对比(1),(3)式子，感觉很相似，可是在推导过程从莫比乌斯反演然后卷积才得出欧拉函数，感觉是绕了一下，不知有没有更好的推导方法。

后来发现原来第一步可以这样化简：

这里简单补一个证明：
$\sum_{k|n}\phi(k)=n$
$\phi(n)$是积性函数，同样设$n=\prod_{i=1}^s p_i^{k_i}$
那么有：
$\sum_{k|n}\phi(k)=\prod_{i=1}^s\sum_{j=0}^{k_i}\phi(p_i^j)=\prod_{i=1}^s(1+\sum_{j=1}^{k_i}(p_i^j-p_i^{j-1}))=\prod_{i=1}^sp_i^{k_i}=n$
所以其实对于gcd的某些问题，用欧拉函数处理要比反演要优。mark。
相关代码：

const int maxn = 1000006;
int prime[maxn];
int phi[maxn];
void getP();//欧拉函数打表,之后要和i^4相乘再求前缀和
int sum2[maxn];//∑i^2

int main()
{
    //...
    int ans=0;
    if(n>m) swap(n,m);
    for(int i=1,j;i<=n;i=j+1)
    {
        j=min(m/(m/i),n/(n/i));
        ans+=(phi[j]-phi[i-1])*sum2[m/i]*sum2[n/i];//这里还没考虑溢出
        ans%=mod;
    }
    //...
}