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}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ij)
(sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}sum_{x|gcd(i,j)}mu(x)ij)
(sum_{d=1}^{n}dsum_{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}dsum_{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}dsum_{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}dsum_{x=1}^{n/d}mu(x)x^{2}g(x))
用数论分块优化(O(nsqrt{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
",(ans+p)%p);
return 0;
}
2、简单的数学题
(sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j))
(sum_{d=1}^{n}dsum_{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(nsqrt{n}))，解法是伪的
所以(mobius)函数没用了，要用(phi)函数
(sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j))
(sum_{i=1}^{n}sum_{j=1}^{n}ijsum_{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
",(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)lfloorfrac{n}{dx} floorlfloorfrac{m}{dx} floor)
(sum_{t=1}^{n}lfloorfrac{n}{t} floorlfloorfrac{m}{t} floorsum_{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)lfloorfrac{n}{t} floorlfloorfrac{m}{t} floor)
#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
",ans);
}
return 0;
}