[bzoj3512]DZY Loves Math IV

给定n,m，求

模10^9+7的值。

听说是杜教筛裸题…

可能是我太鶸了.感觉好难啊..

Orz w_yqts

Orz xudyh

(i,j)=gcd(i,j)

s(i,j)=∑i=1mφ(ni)

ans=∑i=1ns(i,m)

对于n%p=0,有φ(n*p)=φ(n)*p

所以只需要求x=Πp1i的s(x,m).

令d=(n,i)(从某处开始d为变量)

于是有i%d=0,φ(d*i)=φ(i)*d

∵n中最高次数为一次

∴(n/d,i*d)=1,即n/d与i*d互质

∴φ(n∗i)=φ(nd∗i∗d)=φ(nd)∗φ(i)∗d

∵∑d|nφ(d)=n

∴φ(n∗i)=φ(nd)∗φ(i)∗∑e|dφ(e)=φ(nd)∗φ(i)∗∑e|dφ(de)=φ(i)∗∑e|dφ(de)∗φ(nd)=φ(i)∗∑e|dφ(ne)

s(n,m)=∑i=1n∑d|nφ(nd)∗∑i=1⌊md⌋φ(di)=∑d|nφ(nd)∗s(d,⌊md⌋)

s(1,m)用杜教筛求

然后记忆化一下就好了.

#include <bits/stdc++.h>
using namespace std;
#define p 1000000007
#define ll long long
#define N 2500005
map <ll,ll> f,s;
int pn,pr
,pri
,flag
;
ll phi
,Ans
,ph
,ttt;
inline void init()
{
pri[1]=flag[1]=phi[1]=ph[1]=1LL;
for (int i=2;i<N;++i)
{
if (!flag[i]) pr[++pn]=pri[i]=i,phi[i]=i-1;
ph[i]=phi[i];
for (int j=1;j<=pn && pr[j]*i<N;++j)
{
flag[i*pr[j]]=1;
if (i%pr[j]==0) {phi[i*pr[j]]=phi[i]*pr[j];pri[i*pr[j]]=pri[i];break;}
phi[i*pr[j]]=phi[i]*(pr[j]-1);
pri[i*pr[j]]=pri[i]*pr[j];
}
}
for (int i=1;i<N;++i) phi[i]+=phi[i-1];
}
inline ll calc(int n)
{
if (n<N) return phi
;
if (f
) return f
;
ll res=(ll)(n+1)*n>>1;
for (int i=2,pos;i<=n;i=pos+1)
{
pos=n/(n/i);
res-=calc(n/i)*(pos-i+1);
}
return f
=res;
}
inline int sum(int n,int m)
{
if (!m) return 0;
if (s[ttt*(n-1)+m]) return s[ttt*(n-1)+m];
if (n==1) return calc(m)%p;
ll ans=0LL;
for (int i=1;i*i<=n;++i)
if (n%i==0)
{
ans+=ph[n/i]*sum(i,m/i)%p;
if (i*i!=n) ans+=ph[i]*sum(n/i,m/(n/i))%p;
}
return s[ttt*(n-1)+m]=ans%p;
}
int main()
{
init();
ll n,m;
cin>>n>>m;ttt=(ll)m;
ll ans=0LL;
for (int i=1;i<=n;++i)
{
if (pri[i]!=i) Ans[i]=Ans[pri[i]]*(i/pri[i])%p;else
Ans[i]=sum(i,m);
ans=(ans+Ans[i])%p;
}
cout<<ans<<endl;
}