您的位置:首页 > 其它

BZOJ3512 DZY Loves Math IV(杜教筛+线性筛)

2018-09-24 10:41 351 查看

  注意到n很小,考虑枚举i。现在要求的是f(n,m)=Σφ(in) (i=1~m)。显然当n没有平方因子时,φ(in)=φ(i)·φ(n/gcd(i,n))·gcd(i,n)。利用φ*1=id又可得φ(i,n)=φ(i)·Σφ(n/d) (d|gcd(i,n))。改为枚举d就可以得到f(n,m)=Σφ(n/d)*f(d,m/d) (d|n),记忆化搜索求解。n有平方因子时可以发现只要把平方因子提出来最后再乘上就行了,除去平方因子的数可以线性筛得到。

  当n=1时无法继续递归,答案即为φ的前缀和,杜教筛即可。复杂度应该是O(n√m+m2/3)左右,不是很会证。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;
int read()
{
int x=0,f=1;char c=getchar();
while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
return x*f;
}
#define N 100010
#define P 1000000007
int n,m,prime[N<<5],phi[N<<5],p[N<<5],ans=0,cnt=0;
bool flag[N<<5];
map<int,int> f,g
;
int getphi(int n)
{
if (n<(N<<5)) return phi
;
if (f
) return f
;
int s=1ll*n*(n+1)/2%P;
for (int i=2;i<=n;i++)
{
int t=n/(n/i);
s=(s-1ll*(t-i+1)*getphi(n/i)%P+P)%P;
i=t;
}
return f
=s;
}
int calc(int n,int m)
{
if (!m) return 0;
if (n==1) return getphi(m);
if (g
[m]) return g
[m];
int x=n,s=0;n=p
;
for (int i=1;i*i<=n;i++)
if (n%i==0)
{
s=(s+1ll*(getphi(n/i)-getphi(n/i-1)+P)*calc(i,m/i)%P)%P;
if (i*i<n) s=(s+1ll*(getphi(i)-getphi(i-1)+P)*calc(n/i,m/(n/i))%P)%P;
}
s=1ll*s*(x/n)%P;
return g
[m]=s;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("bzoj3512.in","r",stdin);
freopen("bzoj3512.out","w",stdout);
const char LL[]="%I64d\n";
#else
const char LL[]="%lld\n";
#endif
n=read(),m=read();
flag[1]=1,phi[1]=1,p[1]=1;
for (int i=2;i<(N<<5);i++)
{
if (!flag[i]) prime[++cnt]=i,phi[i]=i-1,p[i]=i;
for (int j=1;j<=cnt&&prime[j]*i<(N<<5);j++)
{
flag[prime[j]*i]=1;
if (i%prime[j]==0) {phi[prime[j]*i]=phi[i]*prime[j];p[prime[j]*i]=p[i];break;}
else phi[prime[j]*i]=phi[i]*(prime[j]-1),p[prime[j]*i]=p[i]*prime[j];
}
}
for (int i=1;i<(N<<5);i++) phi[i]=(phi[i-1]+phi[i])%P;
for (int i=1;i<=n;i++)
ans=(ans+calc(i,m))%P;
cout<<ans;
return 0;
}

 

内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: