BZOJ3512：DZY Loves Math IV （杜教筛）

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

---

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

首先有以下四条式子成立：
令$d=gcd(i,j)$：

$$d=\sum_{x|i,x|j}\phi(x)$$

$$\phi(ij)=\phi(\frac{i}{d})\phi(j)d$$

$$\phi(ij)=\frac{\phi(i)\phi(j)d}{\phi(d)}$$

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

$$\phi(\frac{x}{d})\phi(\frac{d}{e})=\phi(\frac{x}{e})$$

第一条式子其实就是$n=\sum_{d|n}\phi(d)$的变形。第二条式子中由于$\frac{i}{d}$与$j$互质，所以$\phi(\frac{i}{d})\phi(j)=\phi(\frac{ij}{d})$，又因为d是$\frac{ij}{d}$的因数，所以等式成立。第三条式子要从$\phi$的计算方法考虑：$\phi(n)=\prod_{i}p_i^{k_i-1}(p_i-1)$，意味着某个质因数$p_i$在n中第一次出现，贡献是$p_i-1$，再出现时贡献是$p_i$。而i和j共有的质因子，在计算$\phi(i)\phi(j)$时乘了两次$p_i-1$，所以要除回$\phi(d)$，再乘以d。第四条式子显然正确。

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

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

$$S(n,m)=y\sum_{i=1}^{m}\phi(ix)$$

用第二条式子代换可得：

$$S(n,m)=y\sum_{i=1}^{m}\phi(i)\phi(\frac{x}{d})d$$

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

$$S(n,m)=y\sum_{i=1}^{m}\phi(i)\phi(\frac{x}{d})\sum_{e|x,e|i}\phi(\frac{d}{e})$$

将e的枚举条件写成$e|x,e|i$是因为$d=gcd(x,i)$，合并了d之后方便些。
然后用第四条式子将后面的两个$\phi$合并，并交换i和e的枚举顺序，可得：

$$S(n,m)=y\sum_{e|x}\phi(\frac{x}{e})\sum_{i=1}^{\lfloor\frac{m}{e}\rfloor}\phi(ie)$$

所以：

$$S(n,m)=y\sum_{e|x}\phi(\frac{x}{e})S(e,\lfloor\frac{m}{e}\rfloor)$$

现在这变成了递归的形式，当n=1时便无法再递归，此时要用杜教筛求出$\phi$的前缀和。由于$\lfloor\frac{\lfloor\frac{m}{x}\rfloor}{y}\rfloor=\lfloor\frac{m}{xy}\rfloor$，于是递归时第二维的值始终是$\lfloor\frac{m}{k}\rfloor$的形式，它只有$2\sqrt{m}$个取值。而且这些值和杜教筛递归时求的值是一样的，所以这个式子复杂度的上界为$O(n\sqrt{m}+n^{\frac{2}{3}})$，对S(n,m)引入记忆化可以在4s内跑过极限数据。从这里也可以看出，杜教筛时记忆化的个数远远小于$2\sqrt{m}$。

这题的关键在于将y提出来，因为S(n,m)中n已经是一个定值。如果不提出y，便无法合并两个$\phi$，就很难计算前缀和。另外，递归时如果用map进行记忆化，不要用下标符号，因为它会返回0，而原先的答案模$10^9+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[N];
    if ( a.count( MP(N,M) ) ) return a[ MP(N,M) ];
    LL x=low[N],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;
}