bzoj 4916: 神犇和蒟蒻【欧拉函数+莫比乌斯函数+杜教筛】

居然扒到了学长出的题
和3944差不多（？），虽然一眼看上去很可怕但是仔细观察发现，对于mu来讲，答案永远是1（对于带平方的，mu值为0，1除外），然后根据欧拉筛的原理，\( \sum_{i=1}^{n}\phi(i^2)=\sum_{i=1}^{n}\phi(i)*i \)，然后就可以正常推了：
设
\[ g(n)=\sum_{i=1}^{n}i\sum_{d=1}^{i}[d|i]\phi(d)=\sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}{6} \]
\[ s(n)=\sum_{i=1}^{n}i\phi(i) \]那么把g展开：
\[ g(n)=\sum_{i=2}^{n}i\sum_{d=1}^{i-1}[d|i]\phi(d)+s(n) \]
\[ s(n)=g(n)-\sum_{i=2}^{n}i\sum_{d=1}^{i-1}[d|i]\phi(d) \]
\[ =g(n)-\sum_{k=2}^{n}k\sum_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}d\phi(d) \]
\[ =g(n)-\sum_{k=2}^{n}k*s(\left \lfloor \frac{n}{k} \right \rfloor) \]
\[ =\frac{n(n+1)(2n+1)}{6}-\sum_{k=2}^{n}k*s(\left \lfloor \frac{n}{k} \right \rfloor) \]
这就是标准的杜教筛递归子问题形式了，直接求解即可。

#include<iostream>
#include<cstdio>
#include<map>
#include<cstring>
using namespace std;
const int N=1000005,m=1000000,mod=1e9+7;
int phi[N],mi[N],q[N],tot,n,k,s[N],ans[N];
bool v[N];
map<long long,int>mp;
int S(int n,int l)
{
    if(l<=1)
        return phi[n*l];
    if(n==1)
    {
        if(l<=m)
            return s[l];
        if(ans[k/l]!=-1)
            return ans[k/l];
        long long re=(long long)l*(l+1)/2%mod;
        for(int i=2,la;i<=l;i=la+1)
        {
            la=l/(l/i);
            if(l/i<=m)
                re=(re-(long long)s[l/i]*(la-i+1)%mod)%mod;
            else
                re=(re-(long long)S(1,l/i)*(la-i+1)%mod)%mod;
        }
        return ans[k/l]=(re%mod+mod)%mod;
    }
    if(mp[(long long)n*mod+l])
        return mp[(long long)n*mod+l];
    long long re=0ll;
    for(int i=1;i*i<=n;i++)
        if(n%i==0)
        {
            re=(re+(long long)phi[n/i]*S(i,l/i)%mod)%mod;
            if(i*i!=n)
                re=re+(long long)phi[i]*S(n/i,l/(n/i))%mod;
        }
    return mp[(long long)n*mod+l]=(re%mod+mod)%mod;
}
// int S(int n,int l)
// {
    // if (l<=1) return phi[n*l];
    // if (n==1)
    // {
        // if (l<=m) return s[l];
        // if (ans[k/l]!=-1) return ans[k/l];
        // int re=(int)l*(l+1)/2%mod;
        // for (int i=2,la;i<=l;i=la+1)
        // {
            // la=l/(l/i);
            // if (l/i<=m) re=re-(int)(la-i+1)*s[l/i]%mod+mod;
            // else re=re-(int)(la-i+1)*S(1,l/i)%mod+mod;
        // }
        // return ans[k/l]=re%mod;
    // }
    // else
    // {
        // if (mp[(int)n*mod+l]) return mp[(int)n*mod+l];
        // int re=0;
        // for (int i=1;i*i<=n;i++)
            // if (n%i==0)
            // {
                // re=re+(int)phi[n/i]*S(i,l/i)%mod;
                // if (i*i!=n) re=re+(int)phi[i]*S(n/i,l/(n/i))%mod;
            // }
        // return mp[(int)n*mod+l]=re%mod;
    // }
// }
int main()
{
    memset(ans,-1,sizeof(ans));
    mi[1]=phi[1]=1;
    for(int i=2;i<=m;i++)
    {
        if(!v[i])
        {
            q[++tot]=i;
            phi[i]=i-1;
            mi[i]=i;
        }
        for(int j=1;j<=tot&&i*q[j]<=m;j++)
        {
            int k=i*q[j];
            v[k]=1;
            if(i%q[j]==0)
            {
                phi[k]=phi[i]*q[j];
                mi[k]=mi[i];
                break;
            }
            phi[k]=phi[i]*(q[j]-1);
            mi[k]=mi[i]*q[j];
        }
    }
    for(int i=1;i<=m;i++)
        s[i]=(s[i-1]+phi[i])%mod;
    scanf("%lld%lld",&n,&k);
    if(n>k)
        swap(n,k);
    long long ans=0ll;
    for(int i=1;i<=n;i++)
        ans=(ans+((long long)i/mi[i]*S(mi[i],k)%mod))%mod;
    printf("%lld\n",(ans%mod+mod)%mod);
    return 0;
}

转载于:https://www.cnblogs.com/lokiii/p/8330788.html