BZOJ3512：DZY Loves Math IV （杜教筛）

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=3512

题目分析：很早就听说这是道神仙题，所以一直没有做。后来洛谷的某场比赛，T1和这题很像，然而我推了1h也化不开ϕ(ij)ϕ(ij)（也许我智商太低），最后在该题取得了0分的好成绩QAQ。于是我决定过来学一学关于ϕϕ的一些新姿势。

首先有以下四条式子成立：

令d=gcd(i,j)d=gcd(i,j)：

d=∑x|i,x|jϕ(x)d=∑x|i,x|jϕ(x)

ϕ(ij)=ϕ(id)ϕ(j)dϕ(ij)=ϕ(id)ϕ(j)d

ϕ(ij)=ϕ(i)ϕ(j)dϕ(d)ϕ(ij)=ϕ(i)ϕ(j)dϕ(d)

若x无平方因子，d为x的约数，e为d的约数，则有：

ϕ(xd)ϕ(de)=ϕ(xe)ϕ(xd)ϕ(de)=ϕ(xe)

第一条式子其实就是n=∑d|nϕ(d)n=∑d|nϕ(d)的变形。第二条式子中由于idid与jj互质，所以ϕ(id)ϕ(j)=ϕ(ijd)ϕ(id)ϕ(j)=ϕ(ijd)，又因为d是ijdijd的因数，所以等式成立。第三条式子要从ϕϕ的计算方法考虑：ϕ(n)=∏ipki−1i(pi−1)ϕ(n)=∏ipiki−1(pi−1)，意味着某个质因数pipi在n中第一次出现，贡献是pi−1pi−1，再出现时贡献是pipi。而i和j共有的质因子，在计算ϕ(i)ϕ(j)ϕ(i)ϕ(j)时乘了两次pi−1pi−1，所以要除回ϕ(d)ϕ(d)，再乘以d。第四条式子显然正确。

然而知道了上面4条式子，我们远远还不足以做出这题。由于这题解法十分脑洞，我直接开始推倒吧：

设S(n,m)=∑mi=1ϕ(in)S(n,m)=∑i=1mϕ(in)，则ans=∑ni=1S(i,m)ans=∑i=1nS(i,m)。由于n比较小，可以直接枚举。考虑怎么计算S(n,m)S(n,m)。设n=xyn=xy，其中x是n的最大无平方因子约数，这个可以通过线性筛直接预处理。通过ϕϕ的计算方式可知ϕ(n)=ϕ(x)yϕ(n)=ϕ(x)y，所以我们把y提出来：

S(n,m)=y∑i=1mϕ(ix)S(n,m)=y∑i=1mϕ(ix)

用第二条式子代换可得：

S(n,m)=y∑i=1mϕ(i)ϕ(xd)dS(n,m)=y∑i=1mϕ(i)ϕ(xd)d

其中d=gcd(x,i)d=gcd(x,i)。然后用第一条式子将后面的d换掉：

S(n,m)=y∑i=1mϕ(i)ϕ(xd)∑e|x,e|iϕ(de)S(n,m)=y∑i=1mϕ(i)ϕ(xd)∑e|x,e|iϕ(de)

将e的枚举条件写成e|x,e|ie|x,e|i是因为d=gcd(x,i)d=gcd(x,i)，合并了d之后方便些。

然后用第四条式子将后面的两个ϕϕ合并，并交换i和e的枚举顺序，可得：

S(n,m)=y∑e|xϕ(xe)∑i=1⌊me⌋ϕ(ie)S(n,m)=y∑e|xϕ(xe)∑i=1⌊me⌋ϕ(ie)

所以：

S(n,m)=y∑e|xϕ(xe)S(e,⌊me⌋)S(n,m)=y∑e|xϕ(xe)S(e,⌊me⌋)

现在这变成了递归的形式，当n=1时便无法再递归，此时要用杜教筛求出ϕϕ的前缀和。由于⌊⌊mx⌋y⌋=⌊mxy⌋⌊⌊mx⌋y⌋=⌊mxy⌋，于是递归时第二维的值始终是⌊mk⌋⌊mk⌋的形式，它只有2m−−√2m个取值。而且这些值和杜教筛递归时求的值是一样的，所以这个式子复杂度的上界为O(nm−−√+n23)O(nm+n23)，对S(n,m)引入记忆化可以在4s内跑过极限数据。从这里也可以看出，杜教筛时记忆化的个数远远小于2m−−√2m。

这题的关键在于将y提出来，因为S(n,m)中n已经是一个定值。如果不提出y，便无法合并两个ϕϕ，就很难计算前缀和。另外，递归时如果用map进行记忆化，不要用下标符号，因为它会返回0，而原先的答案模109+7109+7可能就是0，会增加时间，所以要调用count函数。

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
#include<map>
using namespace std;

const int maxn=1000100;
const long long Mod=1000000007;
typedef long long LL;

LL phi[maxn];
LL sum[maxn];
int low[maxn];

bool vis[maxn];
int prime[maxn];
int cur=0;

#define P pair<int,int>
#define MP(x,y) make_pair(x,(long long)y)

map <P,long long> a;
map <int,long long> b;
int n,m;

void Linear_shaker()
{
phi[1]=1;
low[1]=1;
for (int i=2; i<maxn; i++)
{
if (!vis[i]) phi[i]=i-1,low[i]=i,prime[++cur]=i;
for (int j=1; j<=cur && i*prime[j]<maxn; j++)
{
int k=i*prime[j];
vis[k]=true;
if (i%prime[j])
{
phi[k]=phi[i]*(prime[j]-1);
low[k]=low[i]*prime[j];
}
else
{
phi[k]=phi[i]*prime[j];
low[k]=low[i];
break;
}
}
}
for (int i=1; i<maxn; i++) sum[i]=(sum[i-1]+phi[i])%Mod;
}

LL Sum(LL M)
{
return M*(M+1LL)/2LL%Mod;
}

LL Dfs(int M)
{
if (M<maxn) return sum[M];
if (b.count(M)) return b[M];
LL temp=Sum(M),last;
for (LL i=2; i<=M; i=last+1)
{
last=M/(M/i);
temp=(temp- (last-i+1)*Dfs(M/i)%Mod +Mod)%Mod;
}
b.insert( MP(M,temp) );
return temp;
}

LL Calc(int N,int M)
{
if ( !N || !M ) return 0;
if (N==1) return Dfs(M);
if (M==1) return phi
;
if ( a.count( MP(N,M) ) ) return a[ MP(N,M) ];
LL x=low
,y=N/x,ue=(long long)floor( sqrt( (double)x )+0.1 ),temp=0;
for (LL e=1; e<=ue; e++) if (x%e==0)
{
temp=(temp+ phi[x/e]*Calc(e,M/e)%Mod )%Mod;
if (e!=x/e)
{
e=x/e;
temp=(temp+ phi[x/e]*Calc(e,M/e)%Mod )%Mod;
e=x/e;
}
}
temp=temp*y%Mod;
a.insert( MP( MP(N,M) ,temp) );
return temp;
}

int main()
{
freopen("3512.in","r",stdin);
freopen("3512.out","w",stdout);

Linear_shaker();
scanf("%d%d",&n,&m);
if (n>m) swap(n,m);
LL ans=0;
for (int i=n; i>=1; i--) ans=(ans+ Calc(i,m) )%Mod;
printf("%I64d\n",ans);

return 0;
}