喪心病狂的杜教篩......
這題公式真tm難推……為了這題費了我一個草稿本……
woc……在51Nod上碼LaTeX碼了兩個多小時……
一開始碼完了前半段，剛碼完後半段突然被51Nod吃了，重新碼完後半段之後前半段又被吃了，吓得我趕緊換Notepad++接着寫……
有的細節懶得再碼了，這麼一坨LaTeX估計也夠你們看了……
\begin{equation}
ans=\sum_{i=1}^n\sum_{j=1}^n [i,j]\\
=2\sum_{i=1}^n\sum_{j=1}^i [i,j]-\frac{n(n+1)}2\\
Let\space s(n)=\sum_{i=1}^n\sum_{j=1}^i [i,j],f(n)=\sum_{i=1}^n [i,n]\\
f(n)=\sum_{i=1}^n [i,n]\\
=\sum_{i=1}^n\frac{in}{(i,n)}\\
=n\sum_{i=1}^n\frac i{(i,n)}\\
=n\sum_{d|n}\sum_{i=1}^n[(i,n)=d]\frac i d\\
=n\sum_{d|n}\sum_{i=1}^{\frac n d}[(i,\frac n d)=1]i\\
=n\sum_{d|n}\sum_{i=1}^{d}[(i,d)=1]i\\
=n\sum_{d|n}\frac{\phi(d)d+[d=1]}2\\
=n\frac{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 i\sum_{d|i}\phi(d)d}2\\
=\frac{\frac{n(n+1)}2+\sum_{i=1}^n i\sum_{d|i}\phi(d)d}2\\
=\frac{\frac{n(n+1)}2+\sum_{d=1}^n\phi(d)d\sum_{d|i}i}2\\
=\frac{\frac{n(n+1)}2+\sum_{d=1}^n\phi(d)d^2\sum_{i=1}^{\lfloor\frac n d\rfloor}i}2\\
=\frac{\frac{n(n+1)}2+\sum_{i=1}^n i\sum_{d=1}^{\lfloor\frac n i\rfloor}\phi(d)d^2}2\\
ans=2s(n)-\frac{n(n+1)}2\\
=\sum_{i=1}^n i\sum_{d=1}^{\lfloor\frac n i\rfloor}\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}^n\sum_{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}^{\lfloor\frac n d \rfloor}i^2\\
=\sum_{i=1}^n i^2\sum_{d=1}^{\lfloor\frac n i\rfloor}h(d)\\
=\sum_{i=1}^n i^2 g(\lfloor\frac n i\rfloor)\\
g(n)=\sum_{i=1}^n i^3-\sum_{i=2}^ni^2 g(\lfloor\frac n i\rfloor)
\end{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(\lfloor\frac n i\rfloor)
在外面套上一層分塊不會影響複雜度，利用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 View Code
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 */ 233333333