Bzoj3512 DZY Loves Math IV

Time Limit: 15 Sec  Memory Limit: 128 MB
Submit: 447  Solved: 223

Description

给定n,m，求  模10^9+7的值。

Input

仅一行，两个整数n,m。

Output

仅一行答案。

Sample Input

100000 1000000000

Sample Output

857275582
数据规模：
1<=n<=10^5，1<=m<=10^9，本题共4组数据。

HINT

 

Source

By Jc 鸣谢杜教

 

数学问题 莫比乌斯反演 杜教筛

 

数学题的麻烦在于，如果你拖久了不更博，就不得不再推一遍式子了。

 

注意到n比较小，只有1e5，这是一个可以枚举的范围，我们来看看如果枚举n会怎么样：
设$S(n,m)=\sum_{i=1}^{m} \varphi(i*n)$
有一个结论：
若$|\mu(n)|=1$，即n是无平方数，则：
$$S(n,m)=\sum_{i=1}^{m} \varphi(i) \sum_{d|gcd(i,n)} \varphi(\frac{n}{d})$$
然而博主并不知道如何能想到这个结论。手推一下发现它是对的。
http://www.cnblogs.com/ziliuziliu/p/6508562.html
↑这里有一个关于它的神奇证明

于是当n是无平方数的时候可以这么求：
$$S(n,m)=\sum_{d|n} \varphi(\frac{n}{d}) \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor} \varphi(i*d)$$
$$S(n,m)=\sum_{d|n} \varphi(\frac{n}{d}) S(d,\lfloor \frac{m}{d} \rfloor)$$
记忆化一下即可。

而当$\mu(n)=0$的时候就很方便了：
我们知道

$\varphi(n)= n*\Pi_{i=1}^{k} (1-(1/p_k)) $

其中p数组是n的全部质因数。我们已经可以算出每个质因数都只出现一次时的答案，那么对于剩下的质因数，直接乘进去就行。
也就是说，找到最大的x使得$|\mu(\frac{n}{x})|=1$,
$$S(n,m)=x*S(\frac{n}{x},m)$$
后面这个S用最上面的方法求。
递归边界是$S(1,m)$，也就是 $\mu $ 的前缀和，这个可以用杜教筛算。
最后，$ans=\sum_{i=1}^{n}S(n,m)$

 

#include<iostream>
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<map>
#define LL long long
using namespace std;
const int mod=1e9+7;
const int mxn=5000010;
const int inv2=500000004;
int read(){
int x=0,f=1;char ch=getchar();
while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();}
return x*f;
}
int pri[mxn>>1],cnt=0;
int phi[mxn],mx[mxn];
int smm[mxn];
void init(){
phi[1]=1;mx[1]=1;
for(int i=2;i<mxn;i++){
if(!phi[i]){
pri[++cnt]=i;
phi[i]=i-1;
mx[i]=1;
}
for(int j=1;j<=cnt && pri[j]*i<mxn;j++){
int tmp=pri[j]*i;
if(i%pri[j]==0){
mx[tmp]=mx[i]*pri[j];
phi[tmp]=phi[i]*pri[j];
break;
}
phi[tmp]=phi[i]*(pri[j]-1);
mx[tmp]=mx[i];
}
}
for(int i=1;i<mxn;i++)
smm[i]=((LL)smm[i-1]+phi[i])%mod;
return;
}
map<int,int>mphi;
int Calc_phi(int x){
if(x<mxn)return smm[x];
if(mphi.count(x))return mphi[x];
int res=(LL)x*(x+1)%mod*inv2%mod;
for(int i=2,pos;i<=x;i=pos+1){
int y=x/i; pos=x/y;
res=((LL)res-(pos-i+1)*Calc_phi(y))%mod;
}
res=(res+mod)%mod;
mphi[x]=res;
return res;
}
map<LL,int>mpS;
int solve(int n,int m){
if(!m)return 0;
if(n==1)return Calc_phi(m);
LL tar=(LL)n*mod+m;
if(mpS.count(tar))return mpS[tar];
int res=0;
for(int i=1;i*i<=n;i++){
if(n%i==0){
res=((LL)res+phi[n/i]*(LL)solve(i,m/i))%mod;
if(i*i!=n)res=((LL)res+phi[i]*(LL)solve(n/i,m/(n/i)))%mod;
}
}
mpS[tar]=res;
return res;
}
int n,m;
int main(){
int i,j;
init();
n=read();m=read();
int ans=0;
for(i=1;i<=n;i++)
ans=(ans+(LL)mx[i]*solve(i/mx[i],m))%mod;
ans=(ans+mod)%mod;
printf("%d\n",ans);
return 0;
}