bzoj3512 DZY Loves Math IV

题目链接：bzoj3512

题目大意：

给定n,m,求∑ni=1∑mj=1ϕ(ij)模109+7。

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

话说，样例可真良心。。直接给你个极限的

Input

100000 1000000000

Output

857275582

题解：

杜教筛、数论

设S(n,m)=∑mi=1ϕ(ni)。

那么有Ans=∑ni=1S(i,m)。

膜了一发题解，大概有两种吧。有种用到了一个显然的但是我不会的一个结论。

所以我去学了另一种虽然有点慢但是是我不会上面那个还能做的(期间还膜了膜奥爷爷和男神orz)。

主要是考虑如何弄开S(n,m)=∑mi=1ϕ(ni)。

那么重点就是这个n，将其分解质因数，分类：

- 若n有质因子p的指数大于1，设为pr,r>1，那么ϕ(ni)=pr−1∗ϕ(npr−1i)。设P为n的所有质因子指数超过1的那部分，即nP=p1p2…pn，就是剩下的质因子指数都为1了。所以就有S(n,m)=P∗S(nP,m)。

- 对于满足|μ(n)|=1的n，设其有一质因子p，则有S(n,m)=ϕ(p)∗S(np,m)+S(n,⌊mp⌋)为什么呢？

- 从这个式子S(n,m)=∑mi=1ϕ(ni)出发

- 对于p不能整除i的，ϕ(ni)=ϕ(p)∗ϕ(npi)，即∑i=1nϕ(ni)=ϕ(p)∗S(np,m)=(p−1)∗S(np,m)

- 那么对于p|i，ϕ(ni)=p∗ϕ(npi)，即∑i=1nϕ(ni)=p∗S(np,m)

- 即对所有p|i，都少了S(np,m)=∑mi=1ϕ(npi)，换元把np−>i，即∑⌊mp⌋i=1ϕ(ni)=S(n,⌊mp⌋)

- 综上，就有S(n,m)=ϕ(p)∗S(np,m)+S(n,⌊mp⌋)

然后呢，就是杜教筛的帕了，用其处理求欧拉函数的前缀和。

同样记忆化搜索。

对于S函数，要考虑边界判断，n=1时m=1时等等等等。

也是记忆化递归。

不同的是，这个不能用map来记忆化了，会MLE，状态数太多了。

所以就开个f[i]只存S(i,m)的。

怎么说，这种做法真的贼慢。。十几秒。。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<map>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
#define N 1000100

const LL ny=500000004;
const LL mod=1000000007;
map<int,LL> Phi;bool ispri
;
int cnt,M,pri[N/10],p
,minp
;LL phi
,lim,sum
;
LL mymax(LL x,LL y){return (x>y)?x:y;}
int mymin(int x,int y){return (x<y)?x:y;}
void pre()
{
cnt=0;phi[1]=sum[1]=1;p[1]=1;minp[1]=1;
for (int i=2;i<=lim;i++)
{
if (!ispri[i]) {pri[++cnt]=i;phi[i]=(i-1)%mod;minp[i]=i;p[i]=1;}
for (LL j=1;j<=cnt && i*pri[j]<=lim;j++)
{
ispri[i*pri[j]]=true;
minp[i*pri[j]]=pri[j];
if (i%pri[j]==0)
{
p[i*pri[j]]=p[i]*pri[j]%mod;
phi[i*pri[j]]=phi[i]*pri[j]%mod;
break;
}
p[i*pri[j]]=p[i];
phi[i*pri[j]]=phi[i]*(pri[j]-1)%mod;
}
}
for (LL i=2;i<=lim;i++) sum[i]=(sum[i-1]+phi[i])%mod;
}
LL solve(int n)
{
if (n<=lim) return sum
;
if (Phi.count(n)) return Phi
;
LL ans=(LL)n*(n+1)%mod*ny%mod;int i,r;
for (int i=2;i<=n;i=r+1)
{
r=mymin(n,n/(n/i));
ans-=(LL)(r-i+1)*solve(n/i)%mod;
ans%=mod;if (ans<0) ans+=mod;
}
Phi
=ans;
return ans;
}
LL f
;
LL s(int n,int m)
{
if (m==0) return 0;
if (m==M && f
) return f
;
if (n==1) return solve(m);
if (m==1) return phi
;
LL k=minp
;
LL ans=((LL)(k-1)*s(n/k,m)%mod+s(n,m/k)%mod)%mod;
if (m==M) f
=ans;
return ans;
}
int main()
{
//freopen("a.in","r",stdin);
//freopen("a.out","w",stdout);
int n,m,i;LL ans=0;
scanf("%d%d",&n,&m);M=m;
lim=(LL)trunc(pow(m,2.0/3.0));
lim=mymax(lim,(LL)n);pre();
for (i=1;i<=n;i++)
ans=(ans+(LL)p[i]*s(i/p[i],m)%mod)%mod;
if (ans<0) ans+=mod;
printf("%lld\n",ans);
return 0;
}