Description
去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。
在整理以前的试题时，发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”，其中 表示i的约数个数。他现在长大了，题目也变难了。
求如下表达式的值：
一行一个整数ans，表示答案模1000000007的值。
Sample Input
2
Sample Output
8
HINT
对于100%的数据n <= 10^9。
题解：
解锁新技能：杜教筛。
再复习一下：
若$F(n)=sum_{i=1}^{n}f(i)，g(i)=sum_{j|i}f(i),G(n)=sum_{i=1}^{n}g(i)$，
则有：$G(n)=sum_{i=1}^{n}F(lfloorfrac{n}{i} floor)$。
即：$F(n)=G(n)-sum_{i=2}^{n}F(lfloorfrac{n}{i} floor)$
然后就可以上杜教筛了。
话归正题：
对于本题，其实和bzoj3994式子一样……
$sum_{i=1}^{n}sum_{j=1}^{n}d(ij)$
$=sum_{i=1}^{n}sum_{j=1}^{n}sum_{p=1}^{n^{2}}[p|ij]$
$=sum_{i=1}^{n}sum_{j=1}^{n}sum_{p=1}^{n^{2}}[frac{p}{(p,i)}|j]$
$=sum_{t=1}^{n}sum_{i=1}^{lfloorfrac{n}{t} floor}sum_{j=1}^{n}sum_{p=1}^{lfloorfrac{n^{2}}{t} floor}[p|j]*[(i,p)==1]$
$=sum_{t=1}^{n}sum_{i=1}^{lfloorfrac{n}{t} floor}sum_{p=1}^{lfloorfrac{n^{2}}{t} floor}sum_{d|(i,p)}mu(d)lfloor frac{n}{p} floor$
$=sum_{d=1}^{n}mu(d)sum_{t=1}^{lfloorfrac{n}{d} floor}sum_{i=1}^{lfloorfrac{n}{td} floor}sum_{p=1}^{lfloorfrac{n^{2}}{td} floor}lfloorfrac{n}{pd} floor$
$n>=pd,n^{2}>=tdn$
$ herefore Ans=sum_{d=1}^{n}mu(d)sum_{t=1}^{lfloorfrac{n}{d} floor}sum_{i=1}^{lfloorfrac{n}{td} floor}sum_{p=1}^{lfloorfrac{n}{td} floor}lfloorfrac{n}{pd} floor$
$Ans=sum_{d=1}^{n}mu(d)sum_{t=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{td} floorsum_{p=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{pd} floor$
$设f(i)=sum_{x=1}^{i}lfloorfrac{i}{x} floor$
$ herefore Ans=sum_{d=1}^{n}mu(d)f^2(lfloorfrac{n}{d} floor)$
然后分块+杜教筛即可。
代码：
1 #include<cstdio> 2 using namespace std; 3 typedef long long ll; 4 const ll N=(ll)3e6+1; 5 const ll mod=1000000007ll; 6 short miu[N]; 7 int prim[N/10],num; 8 bool vis[N]; 9 int sum[N]; 10 inline void init(){ 11 miu[1]=sum[1]=1; 12 for(int i=2;i<N;i++){ 13 if(!vis[i]){ 14 prim[++num]=i; 15 miu[i]=-1; 16 }for(int j=1;j<=num&&prim[j]*i<N;j++){ 17 vis[i*prim[j]]=1; 18 if(i%prim[j]==0){ 19 miu[i*prim[j]]=0; 20 break; 21 } 22 else 23 miu[i*prim[j]]=-miu[i]; 24 } 25 // printf("miu[%d]=%d ",i,miu[i]); 26 } 27 for(int i=2;i<N;i++) 28 sum[i]=sum[i-1]+miu[i]; 29 } 30 struct edges{ 31 ll v;int w;edges *last; 32 }edge[N/10],*head[76545];int cnt; 33 const int limit=76543; 34 inline void push(int u,ll v,int w){ 35 edge[++cnt]=(edges){v,w,head[u]};head[u]=edge+cnt; 36 } 37 inline ll Get_sum(ll x){ 38 if(x<N) return sum[x]; 39 int t=x%limit; 40 for(edges *i=head[t];i;i=i->last) 41 if(i->v==x) return i->w; 42 ll ans=1; 43 for(ll i=2,pos;i<=x;i=pos+1){ 44 pos=x/i; 45 pos=x/pos; 46 ans-=(pos-i+1)*Get_sum(x/pos); 47 } 48 push(t,x,ans); 49 return ans; 50 } 51 inline ll Get_F(ll x){ 52 ll ans=0; 53 for(ll i=1,pos;i<=x;i=pos+1){ 54 pos=x/i; 55 pos=x/pos; 56 ll t=x/i; 57 ans+=t*(pos-i+1); 58 ans%=mod; 59 }return ans; 60 } 61 inline ll solve(ll x){ 62 ll ans=0; 63 for(ll i=1,pos;i<=x;i=pos+1){ 64 pos=x/i; 65 pos=x/pos; 66 ll t=Get_F(x/i); 67 ans+=(Get_sum(pos)-Get_sum(i-1))%mod*t%mod*t%mod; 68 ans%=mod; 69 } 70 return (ans%mod+mod)%mod; 71 } 72 int main(){ 73 init(); 74 ll a; 75 scanf("%lld",&a); 76 printf("%lld ",solve(a)); 77 }