这题公式真tm难推……为了这题费了我一个草稿本……
woc……在51Nod上码LaTeX码了两个多小时……
一开始码完了前半段，刚码完后半段突然被51Nod吃了，重新码完后半段之后前半段又被吃了，吓得我赶紧换Notepad++接着写……
有的细节懒得再码了，这么一坨LaTeX估计也够你们看了……
egin{equation}
ans=sum_{i=1}^nsum_{j=1}^n [i,j]\
=2sum_{i=1}^nsum_{j=1}^i [i,j]-frac{n(n+1)}2\
Letspace s(n)=sum_{i=1}^nsum_{j=1}^i [i,j],f(n)=sum_{i=1}^n [i,n]\
f(n)=sum_{i=1}^n [i,n]\
=sum_{i=1}^nfrac{in}{(i,n)}\
=nsum_{i=1}^nfrac i{(i,n)}\
=nsum_{d|n}sum_{i=1}^n[(i,n)=d]frac i d\
=nsum_{d|n}sum_{i=1}^{frac n d}[(i,frac n d)=1]i\
=nsum_{d|n}sum_{i=1}^{d}[(i,d)=1]i\
=nsum_{d|n}frac{phi(d)d+[d=1]}2\
=nfrac{1+sum_{d|n}phi(d)d}2\
s(n)=sum_{i=1}^n f(i)\
=frac{sum_{i=1}^n i(1+sum_{d|i}phi(d)d)}2\
=frac{sum_{i=1}^n i+sum_{i=1}^n isum_{d|i}phi(d)d}2\
=frac{frac{n(n+1)}2+sum_{i=1}^n isum_{d|i}phi(d)d}2\
=frac{frac{n(n+1)}2+sum_{d=1}^nphi(d)dsum_{d|i}i}2\
=frac{frac{n(n+1)}2+sum_{d=1}^nphi(d)d^2sum_{i=1}^{lfloorfrac n d
floor}i}2\
=frac{frac{n(n+1)}2+sum_{i=1}^n isum_{d=1}^{lfloorfrac n i
floor}phi(d)d^2}2\
ans=2s(n)-frac{n(n+1)}2\
=sum_{i=1}^n isum_{d=1}^{lfloorfrac n i
floor}phi(d)d^2\
Let space h(d)=phi(d)d^2,g(n)=sum_{d=1}^nh(d)\
n=sum_{d|n}phi(d)\
n^3=sum_{d|n}phi(d)n^2\
=sum_{d|n}phi(d)d^2(frac n d)^2\
=sum_{d|n}h(d)(frac n d)^2\
sum_{i=1}^n i^3=sum_{i=1}^nsum_{d|i}h(d)(frac i d)^2\
=sum_{d=1}^n h(d)sum_{d|i}(frac i d)^2\
=sum_{d=1}^n h(d)sum_{i=1}^{lfloorfrac n d
floor}i^2\
=sum_{i=1}^n i^2sum_{d=1}^{lfloorfrac n i
floor}h(d)\
=sum_{i=1}^n i^2 g(lfloorfrac n i
floor)\
g(n)=sum_{i=1}^n i^3-sum_{i=2}^ni^2 g(lfloorfrac n i
floor)
end{equation}
然后就是杜教筛的形式了，上杜教筛即可
egin{equation}
sum_{i=1}^n i^3=(frac{n(n+1)}2)^2\
sum_{i=1}^n i^2=frac{n(n+1)(2n+1)}6\
ans=sum_{i=1}^n i g(lfloorfrac n i
floor)
end{equation}
在外面套上一层分块不会影响复杂度，利用g的定义，筛出$phi$之后即可算出较小的g，较大的g直接杜教筛算即可，总复杂度$O(n^{frac 2 3})$
贴两份代码(虽然Python2的代码用Python2和Pypy2交都过不去......)：
1 ''' 2 h(i)=phi(d)*d^2 3 g(i)=sum{h(j)|1<=j<=i} 4 g(n)=sum{i^3|1<=i<=n}-sum{i^2*g(n/i)|2<=i<=n} 5 线筛预处理一部分g，大一些的部分直接上杜教筛即可 6 s_3(n)=s_1(n)^2,s_2(n)=n(n+1)(2n+1)/6 7 ''' 8 p=1000000007 9 table_size=8000000 10 def get_table(n): 11 global phi 12 notp=[False for i in xrange(n+1)] 13 prime=[] 14 cnt=0 15 phi[1]=1 16 for i in xrange(2,n+1): 17 if not notp[i]: 18 prime.append(i) 19 cnt+=1 20 phi[i]=i-1 21 for j in xrange(cnt): 22 if i*prime[j]>n: 23 break 24 notp[i*prime[j]]=True 25 if i%prime[j]: 26 phi[i*prime[j]]=phi[i]*(prime[j]-1) 27 else: 28 phi[i*prime[j]]=phi[i]*prime[j] 29 break 30 for i in xrange(2,n+1): 31 phi[i]=phi[i]*i*i%p 32 phi[i]=(phi[i]+phi[i-1])%p 33 def s1(n): 34 return (n*(n+1)>>1)%p 35 def s2(n): 36 return (n*(n+1)*((n<<1)+1)>>1)/3%p 37 def S(n): 38 if n<table_size: 39 return phi[n] 40 elif hashmap.has_key(n): 41 return hashmap[n] 42 ans=n*(n+1)/2 43 ans*=ans 44 ans%=p 45 i=2 46 while i<=n: 47 last=n/(n/i) 48 #print 'last=%d'%last 49 ans-=(s2(last)-s2(i-1))*S(n/i)%p 50 ans%=p 51 i=last+1 52 if ans<0: 53 ans+=p 54 hashmap[n]=ans 55 return ans 56 n=input() 57 hashmap=dict() 58 table_size=min(table_size,n) 59 phi=[0 for i in xrange(table_size+1)] 60 get_table(table_size) 61 #print 'table OK' 62 ans=0 63 i=1 64 while i<=n: 65 last=n/(n/i) 66 ans+=S(n/i)*(s1(last)-s1(i-1))%p 67 ans%=p 68 i=last+1 69 print ans
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<ext/pb_ds/assoc_container.hpp> 5 #include<ext/pb_ds/hash_policy.hpp> 6 #define s1(n) ((long long)(n)%p*(((n)+1)%p)%p*inv_2%p) 7 #define s2(n) ((long long)(n)%p*(((n)+1)%p)%p*((((long long)(n)%p)<<1)%p+1)%p*inv_6%p) 8 using namespace std; 9 using namespace __gnu_pbds; 10 const int table_size=10000010,maxn=table_size+10,p=1000000007,inv_2=500000004,inv_6=166666668; 11 void get_table(int); 12 int S(long long); 13 bool notp[maxn]={false}; 14 int prime[maxn]={0},phi[maxn]={0}; 15 gp_hash_table<long long,int>hashmap; 16 long long n; 17 int main(){ 18 scanf("%lld",&n); 19 get_table(min((long long)table_size,n)); 20 int ans=0; 21 for(long long i=1,last;i<=n;i=last+1){ 22 last=n/(n/i); 23 ans+=S(n/i)*((s1(last)-s1(i-1))%p)%p; 24 ans%=p; 25 } 26 if(ans<0)ans+=p; 27 printf("%d",ans); 28 return 0; 29 } 30 void get_table(int n){ 31 phi[1]=1; 32 for(int i=2;i<=n;i++){ 33 if(!notp[i]){ 34 prime[++prime[0]]=i; 35 phi[i]=i-1; 36 } 37 for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++){ 38 notp[i*prime[j]]=true; 39 if(i%prime[j])phi[i*prime[j]]=phi[i]*(prime[j]-1); 40 else{ 41 phi[i*prime[j]]=phi[i]*prime[j]; 42 break; 43 } 44 } 45 } 46 for(int i=2;i<=n;i++){ 47 phi[i]=(long long)phi[i]*i%p*i%p; 48 phi[i]=(phi[i]+phi[i-1])%p; 49 } 50 } 51 int S(long long n){ 52 if(n<=table_size)return phi[n]; 53 else if(hashmap.find(n)!=hashmap.end())return hashmap[n]; 54 int ans=s1(n)*s1(n)%p; 55 for(long long i=2,last;i<=n;i=last+1){ 56 last=n/(n/i); 57 ans-=S(n/i)*((s2(last)-s2(i-1))%p)%p; 58 ans%=p; 59 } 60 if(ans<0)ans+=p; 61 return hashmap[n]=ans; 62 } 63 /* 64 h(i)=phi(d)*d^2 65 g(i)=sum{h(j)|1<=j<=i} 66 g(n)=sum{i^3|1<=i<=n}-sum{i^2*g(n/i)|2<=i<=n} 67 ans=sum{i*g(n/i)|1<=i<=n} 68 线筛预处理一部分g，大一些的部分直接上杜教筛即可 69 s_3(n)=s_1(n)^2,s_2(n)=n(n+1)(2n+1)/6 70 */