非常恶心的一道数学题，推式子推到吐血。
光是\(\gcd\)求和我还是会的，但是多了个\(ij\)是什么鬼东西。
\[\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=d]\]
很套路的把后面的\(d\)提出来：
\[\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[\gcd(i,j)=d]=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor} ij[\gcd(i,j)=1]\]
利用\(\sum_{d|n}\mu(d)=[n=1]\)替换掉\([\gcd(i,j)=1]\)就有：
\[\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor} ij[\gcd(i,j)=1]=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor} ij\sum_{x|\gcd(i,j)}\mu(x)\]
然后第一个难点来了，我们把\(x\)弄到前面去枚举，然后利用和式的变换，枚举\(i,j\)是\(x\)的多少倍，然后用分配率凑在一起就有：
\[\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor} ij\sum_{x|\gcd(i,j)}\mu(x)=\sum_{d=1}^nd^3\sum_{x=1}^{\lfloor\frac{n}{d} \rfloor}\mu(x) (\sum_{i=1}^{\lfloor\frac{n}{dx} \rfloor})^2x^2\]
PS:如果上面理解了那么下面就不难了，我们设\(T=dx\)带进去就有：
\[\sum_{d=1}^nd^3\sum_{x=1}^{\lfloor\frac{n}{d} \rfloor}\mu(x) (\sum_{i=1}^{\lfloor\frac{n}{dx} \rfloor})^2x^2=\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2\sum_{d|T}d^3\mu(\frac{T}{d})(\frac{T}{d})^2\]
后面那个东西可以直接提出来就得到：
\[\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2\sum_{d|T}d^3\mu(\frac{T}{d})(\frac{T}{d})^2=\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2\sum_{d|T}dT^2\mu(\frac{T}{d})\]
把\(T^2\)提到前面去：
\[\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2\sum_{d|T}dT^2\mu(\frac{T}{d})=\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2T^2\sum_{d|T}d\mu(\frac{T}{d})\]
后面那个\(\sum_{d|T}d\mu(\frac{T}{d})\)显然就是狄利克雷卷积的形式啊，这就等于\((id\ast\mu)(T)\)
然后回忆一下，\((id\ast\mu)(T)\)不就是\(\phi\)么，因此带进去就有：
\[\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2T^2\sum_{d|T}d\mu(\frac{T}{d})=\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2T^2\phi(T)\]
把\(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i\)用等差数列求和的公式展开就变成：
\[\sum_{T=1}^n(\sum_{i=1}^{\lfloor\frac{n}{T} \rfloor}i)^2T^2\phi(T)=\sum_{T=1}^n(\frac{(1+\lfloor\frac{n}{T} \rfloor)\cdot\lfloor\frac{n}{T} \rfloor}{2})^2T^2\phi(T)\]
如果我们队\(\lfloor\frac{n}{T}\rfloor\)做除法分块，那么现在要求的就是\(T^2\phi(T)\)的前缀和了
考虑用杜教筛的方法化式子，我们社\(f(i)=i^2\phi(i)\)，这个时候我们要找一个\(g\)和\(f\)卷起来可以方便求。
由于\(i^2\)很烦，我们考虑先把它消去。令\(g(i)=i^2\)，则有：
\(h(n)=(f\ast g)(n)=\sum_{d|n} d^2\phi(d)(\frac{n}{d})^2=n^2\ast(\phi\ast \epsilon)(n)=n^3\)
所以用杜教筛的套路式就是\(g(1)f(n)=h(n)-\sum_{d=2}^ng(d)f(\lfloor\frac{n}{d}\rfloor)\)
即\(f(n)=\sum_{i=1}^ni^3-\sum_{d=2}^ng(d)f(\lfloor\frac{n}{d}\rfloor)\)，我们知道\(\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;
}原文地址：https://www.cnblogs.com/cjjsb/p/9965908.html