Luogu3768 简单的数学题
https://www.luogu.com.cn/problem/P3768
欧拉反演/杜教筛
欧拉反演公式:
\[\sum_{d|n} \varphi(d)=n
\]
本题:
\[\sum_{i=1}^n \sum_{j=1}^n ij\gcd(i,j) \\
对\gcd(i,j)进行欧拉反演 \\
\sum_{i=1}^n \sum_{j=1}^n ij\sum_{k|\gcd(i,j)}\varphi(k)\\
=\sum_{i=1}^n \sum_{j=1}^n ij\sum_{k|i,k|j}\varphi(k)\\
=\sum_{k=1}^n \varphi(k)k^2 (\sum_{i=1}^{\lfloor \frac{n}{k}\rfloor}i)^2\\
进行杜教筛\\
设f(n)=\varphi(i)i^2 ,g(n)=n^2\\
(f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})\\
=\sum_{d|n}\varphi(d)d^2(\frac{n}{d})^2\\
=\sum_{d|n}\varphi(d)n^2\\
=n^3\\
\sum_{i=1}^n (f*g)(i)\\
=\sum_{i=1}^n \sum_{d|i} f(d)*g(\frac{i}{d})\\
=\sum_{i=1}^n \sum_{d|i} f(\frac{i}{d})*g(d)\\
=\sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}f(i)\\
=\sum_{d=1}^n g(d)S(\lfloor \frac{n}{d}\rfloor)\\
根据公式\sum_{i=1}^n i^3=\frac{n^2(n+1)^2}{4}\\
\sum_{i=1}^n (f*g)(i)\\
=\sum_{i=1}^n i^3\\
=\frac{n^2(n+1)^2}{4}\\
列出杜教筛核心公式,当d=1时\\
g(1)S(n)=\sum_{i=1}^n (f*g)(i)-\sum_{i=2}^n g(i)*S(\lfloor \frac{n}{i}\rfloor)\\
代入g,(f*g)\\
S(n)=\frac{n^2(n+1)^2}{4}-\sum_{i=2}^n i^2*S(\lfloor \frac{n}{i}\rfloor)\\
再观察原式\\
\sum_{k=1}^n \varphi(k)k^2 (\sum_{i=1}^{\lfloor \frac{n}{k}\rfloor}i)^2\\
在数论分块的同时,设分块区间为[l,r],那么计算(S(r)-S(l-1))*(\sum_{i=1}^{\lfloor \frac{n}{\lfloor \frac{n}{l} \rfloor }\rfloor}i)^2即可
\]
\(map\)可以跑过去,\(unordered\_map\)略快
\(C++ Code:\)
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#define ll long long
#define N 5000000
using namespace std;
ll inv6,cnt,prime[N+10],p,n,phi[N+10],s[N+10];
bool pri[N+10];
map<ll,ll>a;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if (!b)
{
x=1;
y=0;
return a;
}
ll z=exgcd(b,a%b,x,y);
ll t=x;
x=y;
y=t-a/b*y;
return z;
}
ll inv(ll w)
{
ll x,y;
ll c=exgcd(w,p,x,y);
ll v=p/c;
x=(x%v+v)%v;
return x;
}
void Pre()
{
phi[1]=1;
for (ll i=2;i<=N;i++)
{
if (!pri[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for (ll j=1;j<=cnt;j++)
{
ll g=i*prime[j];
if (g>N)
break;
pri[g]=true;
if (i%prime[j]==0)
{
phi[g]=phi[i]*prime[j];
break;
}
phi[g]=phi[i]*phi[prime[j]];
}
}
for (ll i=1;i<=N;i++)
phi[i]%=p;
for (ll i=1;i<=N;i++)
s[i]=(s[i-1]+phi[i]*i%p*i%p)%p;
}
ll _s(ll x)
{
x%=p;
return x*(x+1)%p*(x+x+1)%p*inv6%p;
}
ll S(ll n)
{
if (n<=N)
return s[n];
if (a[n])
return a[n];
ll z=n%p;
ll ans=z*(z+1)/2;
ans%=p;
ans=ans*ans%p;
ll l,r;
for (ll l=2;l<=n;l=r+1)
{
r=n/(n/l);
ll o=S(n/l);
ll e=_s(r)-_s(l-1);
e=(e%p+p)%p;
ans=(ans-o*e%p)%p;
ans=(ans%p+p)%p;
}
ans=(ans%p+p)%p;
return a[n]=ans;
}
int main()
{
scanf("%lld%lld",&p,&n);
Pre();
inv6=inv(6);
ll ans=0;
ll l,r;
for (l=1;l<=n;l=r+1)
{
r=n/(n/l);
ll w=n/l;
w%=p;
w=w*(w+1)/2;
w%=p;
w=w*w%p;
ll e=S(r),o=S(l-1);
e=(e-o)%p;
e=(e%p+p)%p;
ans=(ans+e*w%p)%p;
}
ans=(ans%p+p)%p;
cout << ans << endl;
return 0;
}