题目地址
真是道好题...
题目很直白了，推式子题，开推。
\(\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)·gcd(a[i],a[j])\)
后面的gcd(a[i],a[j])好像不是很好处理，但前面的gcd(i,j)是经典套路，所先把他们分开。
\(\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)==d]·gcd(a[i],a[j])\)
\(\sum_{d=1}^nd\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[gcd(i,j)==1]·gcd(a[i*d],a[j*d])\)
\(\sum_{d=1}^nd\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}\sum_{d'|gcd(i,j)}mu[d']·gcd(a[i*d],a[j*d])\)
\(\sum_{d=1}^nd\sum_{d'=1}^nmu[d']\sum_{i=1}^{n/dd'}\sum_{j=1}^{n/dd'}gcd(a[i*dd'],a[j*dd'])\)
令T=dd'
\(\sum_{T=1}^n\sum_{d'|T}T/d'·mu[d']*\sum_{i=1}^{n/T}\sum_{j=1}^{n/T}gcd(a[i*T],a[j*T])\)
前面是狄利克雷卷积，直接转成欧拉函数。
\(\sum_{T=1}^nphi[T]*\sum_{i=1}^{n/T}\sum_{j=1}^{n/T}gcd(a[i*T],a[j*T])\)
然后推后半部分，首先考虑\(\sum_{i=1}^n\sum_{j=1}^ngcd[a[i],a[j]]怎么推\)
其实差不多的，主要得敢推下去。
\(\sum_{d=1}d\sum_{i=1}^n\sum_{j=1}^n[gcd(a[i],a[j])==d]\)
\(\sum_{d=1}d\sum_{i=1}^n\sum_{j=1}^n[gcd(a[i]/d,a[j]/d)==1]\)
\(\sum_{d=1}d\sum_{i=1}^n\sum_{j=1}^n\sum_{d'|gcd(a[i]/d,a[j]/d)}mu[d']\)
\(\sum_{d=1}\sum_{d'=1}^nmu[d']·d\sum_{i=1}^n\sum_{j=1}^n[dd'|a[i] \& d'd|a[j]]\)
令T=dd'；
\(\sum_{T=1}^n\sum_{d'|T}T/d'·mu[d']\sum_{i=1}^n\sum_{j=1}^n[T|a[i] \& T|a[j]]\)
前后部分合一起。
\(\sum_{T=1}^nphi[T]\sum_{t=1}phi[t]\sum_{i=1}^n\sum_{j=1}^n[T|a[i*T] \& T|a[j*T]]\)
\(\sum_{T=1}^nphi[T]\sum_{t=1}phi[t](\sum_{i=1}^{n/T}[t|a[i*T])^2\)
推完了，然后就是各种预处理欧拉函数，预处理因数，具体看代码。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <queue>
#include <set>
#include <stack>
#include <time.h>
#include <map>
#include <algorithm>
#include <fstream>
//#include <unordered_map>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int maxn = 100000 + 100;
const int INF = 0x7fffffff;
const ll mod = 1e9+7;
const ll mod1 = 998244353;
const ll base = 137;
const double Pi = acos(-1.0);
const int G = 3;
bool vis[maxn];
int prime[maxn];
int phi[maxn];
int ccnt=0;
void euler()
{   phi[1]=1;
	vis[1]=true;
	for(int i=2;i<maxn;i++)
	{
		if(!vis[i])
		{
			prime[++ccnt]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=ccnt&&prime[j]*i<maxn;j++)
		{
			vis[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			phi[i*prime[j]]=phi[i]*phi[prime[j]];
		}
	}
}
int a[maxn];
vector<int>fac[maxn];
vector<int>v;
int cnt[maxn];
int cal(vector<int>v)
{
    int s=0;
    for(auto i:v)
    {
        for(auto j:fac[i])
        {
            s=(1ll*s+1ll*phi[j]*(2ll*cnt[j]+1)%mod)%mod;
            cnt[j]++;
        }
    }
    for(auto i:v)
    {
        for(auto j:fac[i])
        {
            cnt[j]=0;
        }
    }
    return s;
}
int main()
{
    euler();
    for(int i=1;i<=100000;i++)
    {
        for(int j=i;j<=100000;j+=i)
        {
            fac[j].push_back(i);
        }
    }
    int n;
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%d",&a[i]);
    }
    int ans=0;
    for(int i=1;i<=n;i++)
    {
        v.clear();
        for(int j=i;j<=n;j+=i)
        {
            v.push_back(a[j]);
        }
        ans=(1ll*ans+1ll*cal(v)*phi[i]%mod)%mod;
    }
    printf("%d\n",ans);
  //  system("pause");
}