BZOJ3512: DZY Loves Math IV

题意求∑ni=1∑mj=1φ(ij)，(n≤105,m≤109)∑i=1n∑j=1mφ(ij)，(n≤105,m≤109)

直接尝试化简式子好像化不出来，因为n的取值范围不大，考虑枚举i，对于每个i，计算∑mj=1φ(ij)∑j=1mφ(ij)即∑mi=1φ(ni)∑i=1mφ(ni)（下文都用这个柿子叙述）

然后这个式子也不能画的样子，但因为这是φφ，所以当i有因子pikipiki时，可以把piki−1piki−1提到式子前面，于是当n=∏Ki=1pikin=∏i=1Kpiki时，定义D=∏Ki=1piki−1D=∏i=1Kpiki−1，

柿子可以变成D×∑mi=1φ(nDi)D×∑i=1mφ(nDi)，然后对于nDnD，它没有平方因子，考虑一个质因子pipi，如果ii不含pipi，可以将它提到前面变成pi−1pi−1，如果i含有这个因子，将它提到前面应该是pipi，这导致我们少算了∑m/pii=1φ(nDpiipi)=∑m/pii=1φ(nDi)∑i=1m/piφ(nDpiipi)=∑i=1m/piφ(nDi)

然后递归求f(n,m)f(n,m)就好了，要记录一下f(i,m)f(i,m)的情况

code：

#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;

const int maxn = 9e6;
const ll Mod = 1000000007;
int p[700000],pr,phi[maxn],sp[maxn],mp[maxn];
map<int,int>h;
bool v[maxn];
void get_prime()
{
phi[1]=sp[1]=1;
for(int i=2;i<maxn;i++)
{
if(!v[i])
{
p[++pr]=i;
phi[i]=i-1;
mp[i]=i;
}
for(int j=1;j<=pr&&i*p[j]<maxn;j++)
{
int k=i*p[j];
v[k]=true;
mp[k]=p[j];
if(i%p[j]==0)
{
phi[k]=phi[i]*p[j];
break;
}
phi[k]=phi[i]*phi[p[j]];
}
sp[i]=(sp[i-1]+phi[i])%Mod;
}
}
ll N2=(Mod+1)/2,N3=(Mod+1)/3;
int get_s(int x)
{
if(x<maxn) return sp[x];
if(h.count(x)>0) return h[x];
ll ret=(ll)x*(x+1)%Mod*N2%Mod; int nex;
for(int l=2,r=x;l<=r;l=nex+1)
{
nex=r/(r/l);
(ret-=(ll)(nex-l+1)*get_s(r/l)%Mod)%=Mod;
}
if(ret<0)ret+=Mod;
h[x]=ret;
return ret;
}
ll msum;
ll sf[210000]; int M;
ll f(int n,int m)
{
if(m==M&&sf
) return sf
;
if(n==1) return get_s(m);
if(m==1) return phi
;
if(!n||!m) return 0;
ll ps=mp
;
return ((ps-1)%Mod*f(n/mp
,m)%Mod+f(n,m/mp
))%Mod;
}

int main()
{
get_prime();
int n; scanf("%d%d",&n,&M);
msum=get_s(M);
ll ret=msum; sf[1]=msum;
for(int i=2;i<=n;i++)
{
int x=i,y=1;
while(x!=1)
{
int k=mp[x]; y*=k;
while(x%k==0) x/=k;
}
x=i/y;
sf[i]=(ll)x*f(y,M)%Mod;
(ret+=sf[i])%=Mod;
}
printf("%lld\n",ret);

return 0;
}