bzoj 3512: DZY Loves Math IV 杜教筛

题意

给出n,m，求$\sum_{i=1}^n\sum_{j=1}^n\varphi(ij)$。答案模1e9+7.
$n<=10^5,m<=10^9$

分析

看到n辣么小，考虑枚举n。
设

$$S(n,m)=\sum\limits_{i=1}^m\varphi(ni)$$
那么我们要求的就是 $$\sum\limits_{i=1}^nS(i,m)$$
首先有一个结论就是说若 $|\mu(n)|=1$那么有 $$\varphi(nk)=\sum_{d|n}\varphi(\frac{n}{d})\varphi(k)$$
证明：设 $d=gcd(n,k)$，那么有 $$\varphi(nk)=d*\varphi(\frac{n}{d}k)=d*\varphi(\frac{n}{d})\varphi(k)$$
把 $\sum\limits_{i|d}\varphi(i)=d$带进去得 $$\varphi(nk)=\varphi(k)\sum_{i|d}\varphi(i)\varphi(\frac{n}{d})=\varphi(k)\sum_{d|n}\varphi(\frac{n}{d})$$
证毕。
把结论代进 $S(n,m)$的式子里可以得到 $$S(n,m)=\sum_{i=1}^m\varphi(ni)$$
$$=\sum_{i=1}^m\sum_{d|n,d|i}\varphi(i)\varphi(\frac{n}{d})$$
$$=\sum_{d|n}\varphi(\frac{n}{d})\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\varphi(di)$$
$$=\sum_{d|n}\varphi(\frac{n}{d})S(d,\lfloor\frac{m}{d}\rfloor)$$
当 $\mu(n)=0$的时候，我们可以找到最大的 $k$满足$k|n,|\mu(k)|=1$，那么 $S(n,m)=\frac{kS(k,m)}{n}$。
当 $n>1$时可以根据上式求 $S(n,m)$，当 $n=1$时就可以直接杜教筛了。
记得用map记忆化一下。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;

typedef long long LL;

const int N=1000005;
const int MOD=1000000007;

int n,m,phi[N],low[N],tot,prime[N],s[N];
bool not_prime[N];
map<LL,int> ma;

void get_prime(int n)
{
    phi[1]=low[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i,phi[i]=i-1,low[i]=i;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                low[i*prime[j]]=low[i];
                break;
            }
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
            low[i*prime[j]]=low[i]*prime[j];
        }
    }
    for (int i=1;i<=n;i++) s[i]=s[i-1]+phi[i],s[i]-=s[i]>=MOD?MOD:0;
}

LL point(int x,int y)
{
    return (LL)(x-1)*m+(LL)y;
}

int solve(int n,int m)
{
    if (m==0) return 0;
    if (m==1) return phi[n];
    if (n==1)
    {
        if (m<=1000000) return s[m];
        if (ma[point(n,m)]) return ma[point(n,m)];
        int ans=(LL)m*(m+1)/2%MOD;
        for (int i=2,last;i<=m;i=last+1)
        {
            last=m/(m/i);
            ans+=MOD-(LL)(last-i+1)*solve(1,m/i)%MOD;
            ans-=ans>=MOD?MOD:0;
        }
        return ma[point(n,m)]=ans;
    }
    else
    {
        int k=low[n],ans=0;
        if (ma[point(k,m)]) return (LL)ma[point(k,m)]*(n/k)%MOD;
        for (int i=1;i*i<=k;i++)
            if (k%i==0)
            {
                ans+=(LL)phi[k/i]*solve(i,m/i)%MOD;
                ans-=ans>=MOD?MOD:0;
                if (i*i==k) break;
                ans+=(LL)phi[i]*solve(k/i,m/(k/i))%MOD;
                ans-=ans>=MOD?MOD:0;
            }
        ma[point(k,m)]=ans;
        ans=(LL)ans*(n/k)%MOD;
        return ans;
    }
}

int main()
{
    get_prime(1000000);
    scanf("%d%d",&n,&m);
    if (n>m) swap(n,m);
    int ans=0;
    for (int i=1;i<=n;i++) ans+=solve(i,m),ans-=ans>=MOD?MOD:0;
    printf("%d",ans);
    return 0;
}