Luogu P3768 簡單的數學題
非常噁心的一道數學題，推式子推到吐血。
光是\(\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\)
\[\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;
}