Luogu P3768 简单的数学题

非常恶心的一道数学题，推式子推到吐血。

光是(gcd)求和我还是会的，但是多了个(ij)是什么鬼东西。

[sum_{i=1}^nsum_{j=1}^nijgcd(i,j)=sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^nij[gcd(i,j)=d] ]

很套路的把后面的(d)提出来：

[sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^nij[gcd(i,j)=d]=sum_{d=1}^nd^3sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor} ij[gcd(i,j)=1] ]

利用(sum_{d|n}mu(d)=[n=1])替换掉([gcd(i,j)=1])就有：

[sum_{d=1}^nd^3sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor} ij[gcd(i,j)=1]=sum_{d=1}^nd^3sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor} ijsum_{x|gcd(i,j)}mu(x) ]

然后第一个难点来了，我们把(x)弄到前面去枚举，然后利用和式的变换，枚举(i,j)是(x)的多少倍，然后用分配率凑在一起就有：

[sum_{d=1}^nd^3sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{n}{d} floor} ijsum_{x|gcd(i,j)}mu(x)=sum_{d=1}^nd^3sum_{x=1}^{lfloorfrac{n}{d} floor}mu(x) (sum_{i=1}^{lfloorfrac{n}{dx} floor})^2x^2 ]

PS:如果上面理解了那么下面就不难了，我们设(T=dx)带进去就有：

[sum_{d=1}^nd^3sum_{x=1}^{lfloorfrac{n}{d} floor}mu(x) (sum_{i=1}^{lfloorfrac{n}{dx} floor})^2x^2=sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2sum_{d|T}d^3mu(frac{T}{d})(frac{T}{d})^2 ]

后面那个东西可以直接提出来就得到：

[sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2sum_{d|T}d^3mu(frac{T}{d})(frac{T}{d})^2=sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2sum_{d|T}dT^2mu(frac{T}{d}) ]

把(T^2)提到前面去：

[sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2sum_{d|T}dT^2mu(frac{T}{d})=sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2T^2sum_{d|T}dmu(frac{T}{d}) ]

后面那个(sum_{d|T}dmu(frac{T}{d}))显然就是狄利克雷卷积的形式啊，这就等于((idastmu)(T))

然后回忆一下，((idastmu)(T))不就是(phi)么，因此带进去就有：

[sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2T^2sum_{d|T}dmu(frac{T}{d})=sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2T^2phi(T) ]

把(sum_{i=1}^{lfloorfrac{n}{T} floor}i)用等差数列求和的公式展开就变成：

[sum_{T=1}^n(sum_{i=1}^{lfloorfrac{n}{T} floor}i)^2T^2phi(T)=sum_{T=1}^n(frac{(1+lfloorfrac{n}{T} floor)cdotlfloorfrac{n}{T} floor}{2})^2T^2phi(T) ]

如果我们对(lfloorfrac{n}{T} floor)做除法分块，那么现在要求的就是(T^2phi(T))的前缀和了

考虑用杜教筛的方法化式子，我们设(f(i)=i^2phi(i))，这个时候我们要找一个(g)和(f)卷起来可以方便求。

由于(i^2)很烦，我们考虑先把它消去。令(g(i)=i^2)，则有：

(h(n)=(fast g)(n)=sum_{d|n} d^2phi(d)(frac{n}{d})^2=n^2ast(phiast epsilon)(n)=n^3)

所以用杜教筛的套路式就是(g(1)f(n)=h(n)-sum_{d=2}^ng(d)f(lfloorfrac{n}{d} floor))

即(f(n)=sum_{i=1}^ni^3-sum_{d=2}^ng(d)f(lfloorfrac{n}{d} floor))，我们知道(sum_{i=1}^ni^3=(frac{n(n+1)}{2})^2,sum_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6})，所以都可以(O(1))求啦。

CODE

#include<cstdio>
#include<map>
#define RI register int
#define RL register LL
using namespace std;
typedef long long LL;
const int P=10000000;
map <LL,int> sum_f; LL n; int prime[P+5],phi[P+5],cnt,f[P+5],mod,inv2,inv6,ans; bool vis[P+5];
inline void inc(int &x,int y)
{
    if ((x+=y)>=mod) x-=mod;
}
inline void dec(int &x,int y)
{
    if ((x-=y)<0) x+=mod;
}
inline int quick_pow(int x,int p,int mul=1)
{
    for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline int sum(int x)
{
    return 1LL*x*(x+1)%mod*inv2%mod;
}
inline int sqr(int l,int r)
{
    int t=1LL*r*(r+1)%mod*((r<<1LL)+1)%mod*inv6%mod;
    dec(t,1LL*l*(l-1)%mod*((l<<1LL)-1)%mod*inv6%mod); return t;
}
#define Pi prime[j]
inline void Euler(void)
{
    vis[1]=phi[1]=1; RI i,j; for (i=2;i<=P;++i)
    {
        if (!vis[i]) prime[++cnt]=i,phi[i]=i-1;
        for (j=1;j<=cnt&&i*Pi<=P;++j)
        {
            vis[i*Pi]=1; if (i%Pi) phi[i*Pi]=phi[i]*(Pi-1);
            else { phi[i*Pi]=phi[i]*Pi; break; }
        }
    }
    for (i=1;i<=P;++i) f[i]=f[i-1],inc(f[i],1LL*phi[i]*i%mod*i%mod);
}
#undef Pi
inline int F(LL x)
{
    if (x<=P) return f[x]; if (sum_f.count(x)) return sum_f[x];
    int ans=sum(x%mod); ans=1LL*ans*ans%mod; for (RL l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        int t=sqr(l%mod,r%mod);
        dec(ans,1LL*t*F(x/l)%mod);
    }
    return sum_f[x]=ans;
}
int main()
{
    //freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);
    scanf("%d%lld",&mod,&n); Euler(); inv2=quick_pow(2,mod-2);
    inv6=quick_pow(6,mod-2); for (RL l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l); int t=sum(n/l%mod); t=1LL*t*t%mod;
        int Fs=F(r); dec(Fs,F(l-1)); inc(ans,1LL*t*Fs%mod);
    }
    return printf("%d",ans),0;
}