求 S(n) = ∑_i∈[1,n]∑_j∈[1,n] F(i,j),n<=10⁹

    F(i,j) = ∑_d∈[1,i*j] [d | i*j]

1. 化简S(n)

step 1.

    F(i,j) = ∑_d∈[1,i*j] [d/gcd(i,d) | i*j/gcd(i,d)]

             = ∑_d∈[1,i*j] [d/gcd(i,d) | j]

step 2.

    S(n) = ∑_i∈[1,n]∑_j∈[1,n]∑_d∈[1,i*j] [d/gcd(i,d) | j]

           = ∑_i∈[1,n]∑_d∈[1,i*n] [n/(d/gcd(i,d)]

           = ∑_i∈[1,n]∑_d∈[1,i*n] [n*gcd(i,d)/d]

step 3.

    S(n) = ∑_i∈[1,n]∑_d∈[1,i*n]∑_k [gcd(i,d)==k]*[n*k/d]

           = ∑_k∑_i∈[1,n/k]∑_d∈[1,n] [gcd(i,d)==1]*[n/d]

    S(n) = ∑_i∈[1,n]∑_d∈[1,n]∑_k∈[1,n/i] [gcd(i,d)==1]*[n/d]

            = ∑_i∈[1,n]∑_d∈[1,n] [gcd(i,d)==1]*[n/d]*[n/i]

    S(n) = ∑_i∈[1,n]∑_j∈[1,n] [gcd(i,j)==1]*[n/i]*[n/j]

step 4.

    S(n) = ∑_d μ(d)*∑_i∈[1,n/d] ∑_j∈[1,n/d] [n/(id)]*[n/(jd)]

           = ∑_d μ(d)*∑_i∈[1,n/d] [n/(id)]*∑_j∈[1,n/d] [n/(jd)]

    令g(n) = ∑_i∈[1,n] [n/i],则S(n) = ∑_d μ(d)*g([n/d])²

    只计算所有不同的[n/d],可以做到O(n^2/3)

2.0 计算M(n) = ∑_i∈[1,n] μ(i)

    M(n) = ∑_i∈[1,n] μ(i)

            = ∑_i∈[1,n] (∑_d|i μ(d) - ∑_d|i,d<i μ(d))

            = 1 - ∑_i∈[1,n]∑_d|i,d<i μ(d)

            = 1 - ∑_d|i,d<i∑_i∈[1,n] μ(d)

            = 1 - ∑_d∈[1,n] ∑_i∈[2,n/d] μ(d)

            = 1 - ∑_i∈[2,n] ∑_d∈[1,n/i] μ(d)

            = 1 - ∑_i∈[2,n] M([n/i])

    筛出n^2/3以内的M，查询的时候查表O(1)回答.

    对于大于n^2/3的情况，我们用hash表来记忆化.

    这样能做到O(n^2/3)的复杂度.

2.1 计算F(n) = ∑_i∈[1,n] phi(i)

     由n = ∑_d|n phi(d) ,我们可以类似地得到F(n)的计算方法

     F(n) = ∑_i∈[1,n] phi(i)

            = ∑_i∈[1,n] (∑_d|i phi(d) - ∑_d|i,d<i phi(d))

            = n(n+1)/2 - ∑_i∈[1,n]∑_d|i,d<i phi(d)

            = n(n+1)/2 - ∑_i∈[2,n] F([n/i])

3. 码码码

#include <vector>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#define MAXN 1000005
#define SIZE ((1<<18)-1)
#define Mod 1000000007
using namespace std;

typedef long long ll;

namespace Mobius{
	bool ban[MAXN];
	vector<int>p;
	int n,mu[MAXN],ptot;
	inline void linear_shake(int n){
		mu[1]=1;
		for(int i=2;i<=n;++i){
			if(!ban[i]){
				++ptot;
				p.push_back(i);
				mu[i]=Mod-1;
			}
			for(int j=0,k;j<ptot;++j){
				if((ll)i*p[j]>n)break;
				ban[k=i*p[j]]=1;
				if(i%p[j]){
					mu[k]=Mod-mu[i];
				}
				else{
					mu[k]=0;
					break;
				}
			}
		}
		for(int i=2;i<=n;++i){
			mu[i]=(mu[i]+mu[i-1])%Mod;
		}
	}
	int set[MAXN],f[MAXN];
	int head[SIZE+1],tot,next[MAXN];
	inline int M(int x){
		if(!ptot)linear_shake(n=1000000);
		if(x<=n)return mu[x];
		int p=x&SIZE;
		for(int i=head[p];i;i=next[i])
			if(set[i]==x){
				return f[i];
			}
		int res=1;
		for(int i=2,j;i<=x;i=j+1){
			j=x/(x/i);
			res=(res+Mod-ll(j-i+1)*M(x/i)%Mod)%Mod;
		}
		int o=++tot;
		set[o]=x,f[o]=res;
		next[o]=head[p],head[p]=o;
		return res;
	}
};

inline void set_IO(){
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
}

inline int g(int n){
	int res=0;
	for(int i=1,j;i<=n;i=j+1){
		j=n/(n/i);
		res=(res+ll(j-i+1)*(n/i)%Mod)%Mod;
	}
	return res;
}

int n,ans;

int main(){
	set_IO();
	scanf("%d",&n);
	for(int i=1,j,k;i<=n;i=j+1){
		j=n/(n/i);
		k=g(n/i),k=(ll)k*k%Mod;
		k=(ll)k*((Mobius::M(j)+Mod-Mobius::M(i-1))%Mod)%Mod;
		ans=(ans+k)%Mod;
	}
	printf("%d\n",ans);
	return 0;
}