[學習筆記]莫比烏斯反演
1、\(Crash\)的數字表格
\(assume\ n<m\)
\(\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)\)
\(\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]\frac{ij}{d}\)
\(\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]ij\)
\(\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{x|gcd(i,j)}\mu(x)ij\)
\(\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{x|i,x|j}\mu(x)ij\)
\(\sum_{x=1}^{n}\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\mu(x)ij[x|i,x|j]\)
\(assume\ i=ix,j=jx\)
\(\sum_{x=1}^{n}\mu(x)x^{2}\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{m}{dx}}ij\)
\(g(x)=\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{n}{dx}}ij=[\frac{\frac{n}{dx}(\frac{n}{dx}+1)}{2}]^{2}\)
\(\sum_{d=1}^{n}d\sum_{x=1}^{n/d}\mu(x)x^{2}g(x)\)
用數論分塊優化\(O(n\sqrt{n})\)
記得開\(O_{2}\)
\(Code\ Below:\)
#include <bits/stdc++.h> #define ll long long using namespace std; const int maxn=10000000+10; const int p=20101009; ll n,m,prim[maxn],vis[maxn],mu[maxn],sum[maxn],cnt,ans; ll a,b,d; void getmu(ll n){ mu[1]=1; register ll i,j; for(i=2;i<=n;i++){ if(!vis[i]){prim[++cnt]=i;mu[i]=-1;} for(j=1;i*prim[j]<=n&&j<=cnt;j++){ vis[i*prim[j]]=1; if(i%prim[j]==0) break; mu[i*prim[j]]=-mu[i]; } } for(i=1;i<=n;i++) sum[i]=(sum[i-1]+mu[i]*i*i%p)%p; } ll g(ll x){ return ((a/x)*(a/x+1)/2%p)*((b/x)*(b/x+1)/2%p)%p; } int main() { scanf("%lld%lld",&n,&m); if(n>m) swap(n,m); getmu(n); register ll l,r,num; for(d=1;d<=n;d++){ num=0;a=n/d;b=m/d; for(l=1;l<=a;l=r+1){ r=min(a/(a/l),b/(b/l)); num=(num+(sum[r]-sum[l-1])*g(l)%p)%p; } ans=(ans+d*num%p)%p; } printf("%lld\n",(ans+p)%p); return 0; }
2、簡單的數學題
\(\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)\)
\(\sum_{d=1}^{n}d\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==d]ij\)
\(\sum_{d=1}^{n}d^{3}\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[gcd(i,j)==1]ij\)
\(\sum_{d=1}^{n}d^{3}\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}\sum_{x|gcd(i,j)}\mu(x)ij\)
\(\sum_{d=1}^{n}d^{3}\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}\sum_{x=1}^{n/d}\mu(x)[x|i,x|j]ij\)
\(assume\ i=ix,j=jx\)
\(\sum_{d=1}^{n}d^{3}\sum_{x=1}^{n/d}\mu(x)x^{2}\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{n}{dx}}ij\)
\(g(x)=\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{n}{dx}}ij=[\frac{\frac{n}{dx}(\frac{n}{dx}+1)}{2}]^{2}\)
\(\sum_{d=1}^{n}d^{3}\sum_{x=1}^{n/d}\mu(x)x^{2}g(x)\)
突然發現時間複雜度\(O(n\sqrt{n})\)，解法是偽的
所以\(mobius\)卷積沒用了，要用\(phi\)卷積
\(\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)\)
\(\sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{k|i,k|j}\varphi(k)\)
\(\sum_{k=1}^{n}\varphi(k)\sum_{k|i}\sum_{k|j}ij\)
\(\sum_{k=1}^{n}\varphi(k)\sum_{i=1}^{n}\sum_{j=1}^{n}[k|i,k|j]ij\)
\(\sum_{k=1}^{n}\varphi(k)k^{2}\sum_{i=1}^{n/k}\sum_{j=1}^{n/k}ij\)
\(g(x)=(\sum_{i=1}^{x}i)^{2}=[\frac{x(x+1)}{2}]^{2}\)
\(\sum_{k=1}^{n}\varphi(k)k^{2}g(n/k)^{2}\)
#include <bits/stdc++.h>
#define ll long long
#define res register ll
using namespace std;
const int maxn=7500000+10;
ll p,n,prim[maxn],vis[maxn],phi[maxn],cnt,ans,inv2,inv6;
map<ll,ll> Phi;
ll fast_pow(ll a,ll b){
ll ret=1;
for(;b;b>>=1,a=a*a%p)
ret=ret*(b&1?a:1)%p;
return ret;
}
void getphi(ll n){
phi[1]=1;
res i,j;
for(i=2;i<=n;i++){
if(!vis[i]){prim[++cnt]=i;phi[i]=i-1;}
for(j=1;i*prim[j]<=n&&j<=cnt;j++){
vis[i*prim[j]]=1;
if(i%prim[j]==0){
phi[i*prim[j]]=phi[i]*prim[j]%p;
break;
}
phi[i*prim[j]]=phi[i]*(prim[j]-1)%p;
}
}
for(i=1;i<=n;i++) phi[i]=(i*i%p*phi[i]%p+phi[i-1])%p;
}
ll g(ll x){
x%=p;
return (x*(x+1)/2%p)*(x*(x+1)/2%p)%p;
}
ll f(ll x){
x%=p;
return x*(x+1)%p*(2*x+1)%p*inv6%p;
}
ll calc_phi(res n){
if(n<=7500000) return phi[n];
if(Phi[n]) return Phi[n];
res ans=g(n);
for(res l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans=(ans-(f(r)-f(l-1)+p)%p*calc_phi(n/l)%p)%p;
}
return Phi[n]=(ans+p)%p;
}
int main()
{
scanf("%lld%lld",&p,&n);
getphi(7500000);
inv2=fast_pow(2,p-2);
inv6=fast_pow(6,p-2);
for(res l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans=(ans+(calc_phi(r)-calc_phi(l-1)+p)%p*g(n/l)%p)%p;
}
printf("%lld\n",(ans+p)%p);
return 0;
}3、於神之怒加強版
\(assume\ n<m\)
\(\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)^{k}\)
\(\sum_{d=1}^{n}d^{k}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]\)
\(\sum_{d=1}^{n}d^{k}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]\)
\(\sum_{d=1}^{n}d^{k}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{x|gcd(i,j)}\mu(x)\)
\(\sum_{d=1}^{n}d^{k}\sum_{x=1}^{n}\mu(x)\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[x|i,x|j]\)
\(\sum_{d=1}^{n}d^{k}\sum_{x=1}^{n/d}\mu(x)\lfloor\frac{n}{dx}\rfloor\lfloor\frac{m}{dx}\rfloor\)
\(\sum_{t=1}^{n}\lfloor\frac{n}{t}\rfloor\lfloor\frac{m}{t}\rfloor\sum_{d|t}d^{k}\mu(\frac{t}{d})\)
\(f(x)=\sum_{d|t}d^{k}\mu(\frac{t}{d})\)
\(\sum_{t=1}^{n}f(t)\lfloor\frac{n}{t}\rfloor\lfloor\frac{m}{t}\rfloor\)
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=5000000+10;
const int p=1e9+7;
ll n,m,k,prim[maxn],vis[maxn],f[maxn],cnt,ans;
ll fast_pow(ll a,ll b){
ll ret=1;
for(;b;b>>=1,a=a*a%p)
ret=ret*(b&1?a:1)%p;
return ret;
}
void sieve(ll n){
f[1]=1;
ll i,j;
for(i=2;i<=n;i++){
if(!vis[i]){prim[++cnt]=i;f[i]=fast_pow(i,k)-1;}
for(j=1;i*prim[j]<=n&&j<=cnt;j++){
vis[i*prim[j]]=1;
if(i%prim[j]==0){
f[i*prim[j]]=f[i]*(f[prim[j]]+1)%p;
break;
}
f[i*prim[j]]=f[i]*f[prim[j]]%p;
}
}
for(i=1;i<=n;i++) f[i]=(f[i]+f[i-1])%p;
}
int main()
{
ll T,l,r;
scanf("%lld%lld",&T,&k);
sieve(5000000);
while(T--){
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
ans=0;
for(l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+(n/l)*(m/l)%p*(f[r]-f[l-1]+p)%p)%p;
}
printf("%lld\n",ans);
}
return 0;
}