Description
去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了。
![](http://www.lydsy.com/JudgeOnline/upload/201507/abc.jpg)
Input
第一行一个整数n。
Output
一行一个整数ans,表示答案模1000000007的值。
Sample Input
Sample Output
HINT
对于100%的数据n <= 10^9。
数学问题 莫比乌斯反演 脑洞题
$ans(n)=sum_{i=1}^{n} sum_{j=1}^{n} f(i,j)$
$=sum_{i=1}^{n} sum_{j=1}^{n} sum_{k=1}^{n^2} [frac{k}{gcd(i,j)}|j]$
(枚举gcd)
$=sum_{d=1}^{n} d sum_{i=1}^{lfloor frac{n}{d}
floor}sum_{j=1}^{lfloor frac{n}{d}
floor}sum_{k=1}^{lfloor frac{n^2}{d}
floor} [k|j][gcd(k,i)=1] $
(消掉j)
$=sum_{d=1}^{n} d sum_{i=1}^{lfloor frac{n}{d}
floor}sum_{k=1}^{lfloor frac{n^2}{d}
floor} d lfloor frac{n}{dk}
floor [gcd(k,i)=1] $
(约掉一个d)
$=sum_{d=1}^{n} sum_{i=1}^{lfloor frac{n}{d}
floor}sum_{k=1}^{lfloor frac{n^2}{d}
floor} lfloor frac{n}{k}
floor [gcd(k,i)=1] $
$=sum_{d=1}^{n} sum_{i=1}^{lfloor frac{n}{d} floor}sum_{k=1}^{lfloor frac{n^2}{d} floor} lfloor frac{n}{k} floor sum_{t|(ik)} mu(t)$
$= sum_{t=1}^{1} mu(t) sum_{d=1}^{n} sum_{i=1}^{lfloor frac{n}{dt}
floor}sum_{k=1}^{lfloor frac{n^2}{dt}
floor} lfloor frac{n}{kt}
floor$
(注意到i对后面没有影响,把求和变为乘)
$= sum_{t=1}^{1} mu(t) sum_{d=1}^{n} lfloor frac{n}{dt}
floorsum_{k=1}^{lfloor frac{n^2}{dt}
floor} lfloor frac{n}{kt}
floor$
(注意到d在大于$ lfloor frac{n}{t}
floor$ 时,后面全是0,所以缩小求和上界,后面的k同理)
$= sum_{t=1}^{1} mu(t) sum_{d=1}^{lfloor frac{n}{t}
floor} lfloor frac{n}{dt}
floorsum_{k=1}^{lfloor frac{n}{t}
floor} lfloor frac{n}{kt}
floor$
$= sum_{t=1}^{1} mu(t) (sum_{d=1}^{lfloor frac{n}{t}
floor} lfloor frac{n}{dt}
floor)^2 $
$mu$的前缀和可以用bzoj3944的方法筛出来,后面的平方项可以分块计算。
也就是说分块套分块就可以愉快地出解了。
隔壁Sfailsth告诉我可以直接用Bzoj3994的结论来做
↓也就是从这个式子开始:
$sum_{i=1}^{N}sum_{j=1}^{M}d(ij)=sum_{i=1}^{N}sum_{j=1}^{M}lfloorfrac{N}{i}
floorlfloorfrac{M}{j}
floorlbrack gcd(i,j)=1
brack $
确实很方便
(做完这道题之后,陪伴了我们一年的权限号过期了)
1 #include<iostream> 2 #include<algorithm> 3 #include<cstring> 4 #include<cstdio> 5 #include<cmath> 6 #include<map> 7 #define LL long long 8 using namespace std; 9 const int mod=1e9+7; 10 const int mxn=5000010; 11 int read(){ 12 int x=0,f=1;char ch=getchar(); 13 while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();} 14 while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();} 15 return x*f; 16 } 17 int pri[mxn],mu[mxn],cnt=0; 18 bool vis[mxn]; 19 void init(){ 20 mu[1]=1; 21 for(int i=2;i<mxn;i++){ 22 if(!vis[i]){ 23 pri[++cnt]=i; 24 mu[i]=-1; 25 } 26 for(int j=1;j<=cnt && (LL)pri[j]*i<mxn;j++){ 27 int tmp=pri[j]*i; 28 vis[tmp]=1; 29 if(i%pri[j]==0){vis[tmp]=1;mu[tmp]=0;break;} 30 mu[tmp]=-mu[i]; 31 } 32 } 33 for(int i=2;i<mxn;i++){mu[i]=mu[i-1]+mu[i];} 34 return; 35 } 36 map<int,LL>mpmu; 37 int F(int x){ 38 if(x<mxn)return mu[x]; 39 if(mpmu.count(x))return mpmu[x]; 40 int res=1; 41 for(int i=2,pos;i<=x;i=pos+1){ 42 int y=x/i; 43 pos=x/y; 44 res=((LL)res-(LL)(pos-i+1)*F(y)%mod+mod)%mod; 45 } 46 mpmu[x]=res; 47 return res; 48 } 49 int calc(int x){ 50 int res=0; 51 for(int i=1,pos;i<=x;i=pos+1){ 52 int y=x/i; 53 pos=x/y; 54 res=(res+(pos-i+1)*(LL)y%mod)%mod; 55 } 56 return (LL)res*res%mod; 57 } 58 int main(){ 59 init(); 60 int n=read(); 61 LL ans=0; 62 for(int i=1,pos;i<=n;i=pos+1){ 63 int y=n/i;pos=n/y; 64 ans=((LL)ans+(((LL)F(pos)-F(i-1)%mod))*calc(y)%mod)%mod; 65 } 66 ans=(ans%mod+mod)%mod; 67 printf("%lld ",ans); 68 return 0; 69 }