G. GCD Festival（莫比乌斯反演）

题目地址
真是道好题…
题目很直白了，推式子题，开推。
${\sum }_{i=1}^n{\sum }_{j=1}^{n}gcd(i,j)\cdot gcd(a[i],a[j])$
后面的gcd(a[i],a[j])好像不是很好处理，但前面的gcd(i,j)是经典套路，所先把他们分开。
${\sum }_{d=1}^{n}d{\sum }_{i=1}^n{\sum }_{j=1}^n[gcd(i,j)==d]\cdot gcd(a[i],a[j])$
${\sum }_{d=1}^{n}d{\sum }_{i=1}^{n/d}{\sum }_{j=1}^{n/d}[gcd(i,j)==1]\cdot gcd(a[i\ast d],a[j\ast d])$
${\sum }_{d=1}^{n}d{\sum }_{i=1}^{n/d}{\sum }_{j=1}^{n/d}{\sum }_{d^{\prime }|gcd(i,j)}mu[d^{\prime }]\cdot gcd(a[i\ast d],a[j\ast d])$
${\sum }_{d=1}^{n}d{\sum }_{d^{\prime }=1}^{n}mu[d^{\prime }]{\sum }_{i=1}^{n/d{d}^{\prime }}{\sum }_{j=1}^{n/d{d}^{\prime }}gcd(a[i\ast d{d}^{\prime }],a[j\ast d{d}^{\prime }])$
令T=dd’
${\sum }_{T=1}^n{\sum }_{d^{\prime }|T}T/d^{\prime }\cdot mu[d^{\prime }]\ast {\sum }_{i=1}^{n/T}{\sum }_{j=1}^{n/T}gcd(a[i\ast T],a[j\ast T])$
前面是狄利克雷卷积，直接转成欧拉函数。
${\sum }_{T=1}^{n}phi[T]\ast {\sum }_{i=1}^{n/T}{\sum }_{j=1}^{n/T}gcd(a[i\ast T],a[j\ast T])$
然后推后半部分，首先考虑${\sum }_{i=1}^n{\sum }_{j=1}^{n}gcd[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^{\prime }|gcd(a[i]/d,a[j]/d)}mu[d^{\prime }]$
${\sum }_{d=1}{\sum }_{d^{\prime }=1}^{n}mu[d^{\prime }]\cdot d{\sum }_{i=1}^n{\sum }_{j=1}^n[d{d}^{\prime }|a[i]\&d^{\prime }d|a[j]]$
令T=dd’；
${\sum }_{T=1}^n{\sum }_{d^{\prime }|T}T/d^{\prime }\cdot mu[d^{\prime }]{\sum }_{i=1}^n{\sum }_{j=1}^n[T|a[i]\&T|a[j]]$
前后部分合一起。
${\sum }_{T=1}^{n}phi[T]{\sum }_{t=1}phi[t]{\sum }_{i=1}^n{\sum }_{j=1}^n[T|a[i\ast T]\&T|a[j\ast T]]$
${\sum }_{T=1}^{n}phi[T]{\sum }_{t=1}phi[t]({\sum }_{i=1}^{n/T}[t|a[i\ast 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");
}