矩阵求和

经过重重笔试面试的考验，小明成功进入 Macrohard 公司工作。
今天小明的任务是填满这么一张表：
表有 n 行 n 列，行和列的编号都从1算起。
其中第 i 行第 j 个元素的值是 gcd(i, j)的平方，
gcd 表示最大公约数，以下是这个表的前四行的前四列：
1 1 1 1
1 4 1 4
1 1 9 1
1 4 1 16

小明突然冒出一个奇怪的想法，他想知道这张表中所有元素的和。
由于表过于庞大，他希望借助计算机的力量。

「输入格式」
一行一个正整数 n 意义见题。

「输出格式」
一行一个数，表示所有元素的和。由于答案比较大，请输出模 (10^9 + 7)(即：十亿零七) 后的结果。

「样例输入」
4

「样例输出」
48

「数据范围」
对于 30% 的数据，n <= 1000
存在 10% 的数据，n = 10^5
对于 60% 的数据，n <= 10^6
对于 100% 的数据，n <= 10^7

资源约定：
峰值内存消耗（含虚拟机） < 256M
CPU消耗 < 2000ms

请严格按要求输出，不要画蛇添足地打印类似：“请您输入...” 的多余内容。

注意：
main函数需要返回0;
只使用ANSI C/ANSI C++ 标准;
不要调用依赖于编译环境或操作系统的特殊函数。
所有依赖的函数必须明确地在源文件中 #include <xxx>
不能通过工程设置而省略常用头文件。

提交程序时，注意选择所期望的语言类型和编译器类型。

莫比乌斯反演不会，但是可以明白下面这种解法。我们需要求n以内一共有几个k*k(1<=k<=n),也就是说gcd(i,j)=k的(i,j)共有多少对，可以转化为n/k以内gcd(i,j)=1的对数，这样就可以用欧拉筛法求一下欧拉函数，然后前缀和就是满足条件的对数，然后求总和。

代码：

#include <iostream>
#define mod 1000000007
#define MAX 10000005
using namespace std;
typedef long long ll;
ll phi[MAX],s[MAX],prime[MAX];
bool vis[MAX];
void init(int n) {
    s[1] = phi[1] = 1;
    for(int i = 2;i <= n;i ++){
        if(!vis[i]) {
            phi[i] = i - 1;
            prime[++ prime[0]] = i;
        }
        for(int j = 1;j <= prime[0] && i * prime[j] <= n;j ++) {
            vis[i * prime[j]] = true;
            if(i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
        s[i] = s[i - 1] + 2 * phi[i];
    }
}

int main(){
    int n;
    cin>>n;
    init(n);
    ll ans = 0;
    for(int i = 1;i <= n;i ++)
        ans = (ans + s[n / i] * i % mod * i) % mod;
    cout<<ans<<endl;
    return 0;
}