math
求 S(n) = ∑i∈[1,n]∑j∈[1,n] F(i,j),n<=109
F(i,j) = ∑d∈[1,i*j] [d | i*j]
1. 化简S(n)
step 1.
F(i,j) = ∑d∈[1,i*j] [d/gcd(i,d) | i*j/gcd(i,d)]
= ∑d∈[1,i*j] [d/gcd(i,d) | j]
step 2.
S(n) = ∑i∈[1,n]∑j∈[1,n]∑d∈[1,i*j] [d/gcd(i,d) | j]
= ∑i∈[1,n]∑d∈[1,i*n] [n/(d/gcd(i,d)]
= ∑i∈[1,n]∑d∈[1,i*n] [n*gcd(i,d)/d]
step 3.
S(n) = ∑i∈[1,n]∑d∈[1,i*n]∑k [gcd(i,d)==k]*[n*k/d]
= ∑k∑i∈[1,n/k]∑d∈[1,n] [gcd(i,d)==1]*[n/d]
S(n) = ∑i∈[1,n]∑d∈[1,n]∑k∈[1,n/i] [gcd(i,d)==1]*[n/d]
= ∑i∈[1,n]∑d∈[1,n] [gcd(i,d)==1]*[n/d]*[n/i]
S(n) = ∑i∈[1,n]∑j∈[1,n] [gcd(i,j)==1]*[n/i]*[n/j]
step 4.
S(n) = ∑d μ(d)*∑i∈[1,n/d] ∑j∈[1,n/d] [n/(id)]*[n/(jd)]
= ∑d μ(d)*∑i∈[1,n/d] [n/(id)]*∑j∈[1,n/d] [n/(jd)]
令g(n) = ∑i∈[1,n] [n/i],则S(n) = ∑d μ(d)*g([n/d])2
只计算所有不同的[n/d],可以做到O(n2/3)
2.0 计算M(n) = ∑i∈[1,n] μ(i)
M(n) = ∑i∈[1,n] μ(i)
= ∑i∈[1,n] (∑d|i μ(d) - ∑d|i,d<i μ(d))
= 1 - ∑i∈[1,n]∑d|i,d<i μ(d)
= 1 - ∑d|i,d<i∑i∈[1,n] μ(d)
= 1 - ∑d∈[1,n] ∑i∈[2,n/d] μ(d)
= 1 - ∑i∈[2,n] ∑d∈[1,n/i] μ(d)
= 1 - ∑i∈[2,n] M([n/i])
筛出n2/3以内的M，查询的时候查表O(1)回答.
对于大于n2/3的情况，我们用hash表来记忆化.
这样能做到O(n2/3)的复杂度.
2.1 计算F(n) = ∑i∈[1,n] phi(i)
由n = ∑d|n phi(d) ,我们可以类似地得到F(n)的计算方法
F(n) = ∑i∈[1,n] phi(i)
= ∑i∈[1,n] (∑d|i phi(d) - ∑d|i,d<i phi(d))
= n(n+1)/2 - ∑i∈[1,n]∑d|i,d<i phi(d)
= n(n+1)/2 - ∑i∈[2,n] F([n/i])
3. 码码码
#include <vector>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#define MAXN 1000005
#define SIZE ((1<<18)-1)
#define Mod 1000000007
using namespace std;
typedef long long ll;
namespace Mobius{
bool ban[MAXN];
vector<int>p;
int n,mu[MAXN],ptot;
inline void linear_shake(int n){
mu[1]=1;
for(int i=2;i<=n;++i){
if(!ban[i]){
++ptot;
p.push_back(i);
mu[i]=Mod-1;
}
for(int j=0,k;j<ptot;++j){
if((ll)i*p[j]>n)break;
ban[k=i*p[j]]=1;
if(i%p[j]){
mu[k]=Mod-mu[i];
}
else{
mu[k]=0;
break;
}
}
}
for(int i=2;i<=n;++i){
mu[i]=(mu[i]+mu[i-1])%Mod;
}
}
int set[MAXN],f[MAXN];
int head[SIZE+1],tot,next[MAXN];
inline int M(int x){
if(!ptot)linear_shake(n=1000000);
if(x<=n)return mu[x];
int p=x&SIZE;
for(int i=head[p];i;i=next[i])
if(set[i]==x){
return f[i];
}
int res=1;
for(int i=2,j;i<=x;i=j+1){
j=x/(x/i);
res=(res+Mod-ll(j-i+1)*M(x/i)%Mod)%Mod;
}
int o=++tot;
set[o]=x,f[o]=res;
next[o]=head[p],head[p]=o;
return res;
}
};
inline void set_IO(){
freopen("math.in","r",stdin);
freopen("math.out","w",stdout);
}
inline int g(int n){
int res=0;
for(int i=1,j;i<=n;i=j+1){
j=n/(n/i);
res=(res+ll(j-i+1)*(n/i)%Mod)%Mod;
}
return res;
}
int n,ans;
int main(){
set_IO();
scanf("%d",&n);
for(int i=1,j,k;i<=n;i=j+1){
j=n/(n/i);
k=g(n/i),k=(ll)k*k%Mod;
k=(ll)k*((Mobius::M(j)+Mod-Mobius::M(i-1))%Mod)%Mod;
ans=(ans+k)%Mod;
}
printf("%d\n",ans);
return 0;
}
help
求 ∑n∈[1,N]∑m∈[1,M]∑k∈[0,m) [(nk+x)/m] (N,M<=500000,x为实数且x∈[0,100000])
0.首先令x = [x]显然是不影响结果的
1.先考虑如何化简∑k∈[0,m) [(nk+x)/m]
首先显然有：[(nk+x)/m] = [(nk%m+x)/m]+[nk/m]
[(nk+x)/m] = [(nk%m+x)/m]+(nk-nk%m)/m
则∑k∈[0,m) [(nk+x)/m] = ∑k∈[0,m) [(nk%m+x)/m] + ∑k∈[0,m) nk/m + ∑k∈[0,m) nk%m/m
令d=gcd(n,m)，则∑k∈[0,m) [(nk%m+x)/m] = d([x/m]+[(x+d)/m]+...+[(x+m-d)/m])
= d([(x/d)/(m/d)]+[((x/d)+1)/(m/d)]+...+[((x/d)+m/d-1)/(m/d)])
= d([(x/d)(m/d)/(m/d)] (注)
= d[x/d]
∑k∈[0,m) nk%m = d(0+d+...+m-d)
= (m-d)/2
∑k∈[0,m) [(nk+x)/m] = d[x/d]+n(m-1)/2+(d-m)/2
= d[x/d]+(n-1)(m-1)/2+(d-1)/2
由此可见：∑k∈[0,m) [(nk+x)/m] = ∑k∈[0,n) [(mk+x)/n]
注：
对于正整数m和实数x, [mx] = [x] + [x + 1/m] + ... + [x + (m-1)/m]
证明：
x = [x] + {x}
设：k = [m{x}]
则k/m <= m{x} < (k+1)/m
[mx] = [m([x] + {x})]
= m[x] + [m{x}]
= m[x] + k
[x]+[x+1/m]+...+[x+(m-1)/m]
= m[x] + [{x}+1/m] + ... + [{x}+(m-1)/m]
= m[x] + k
[mx] = [x] + [x + 1/m] + ... + [x + (m-1)/m]
2. 再考虑如何计算∑n∈[1,N]∑m∈[1,M] d[x/d]+(n-1)(m-1)/2+(d-1)/2 (d=gcd(n,m))
分为3部分计算：
1) ∑n∈[1,N]∑m∈[1,M] d[x/d]
2) ∑n∈[1,N]∑m∈[1,M] (n-1)(m-1)/2
3) ∑n∈[1,N]∑m∈[1,M] (d-1)/2
1) 和 3) 考虑容斥，枚举d , 再枚举d的倍数，复杂度O(nlogn)
2) 可以直接O(1)计算
3.码码码
#include <vector>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#define P 998244353
#define r2 499122177
#define MAXN 500005
using namespace std;
#define g() getchar()
template<class Q>inline void Scan(Q&x){
char c; int f=1;
while(c=g(),c<48||c>57)if(c=='-')f=-1;
for(x=0;c>47&&c<58;c=g())x=10*x+c-48;
x*=f;
}
typedef long long ll;
ll n,m,x;
ll p1,p2,p3;
bool ban[MAXN];
vector<int>p;
int mu[MAXN],tot;
inline void get_mu(int n){
mu[1]=1;
for(int i=2;i<=n;++i){
if(!ban[i]){
++tot;
p.push_back(i);
mu[i]=-1;
}
for(int j=0;j<tot;++j){
if((ll)i*p[j]>n)break;
int x=i*p[j];
ban[x]=1;
if(i%p[j]){
mu[x]=P-mu[i];
}
else{
mu[x]=0;
break;
}
}
}
}
inline void set_IO(){
freopen("help.in","r",stdin);
freopen("help.out","w",stdout);
}
int main(){
set_IO();
Scan(n),Scan(m);
if(n>m)swap(n,m);
get_mu(n);
Scan(x);
p3=(P-n*m%P)%P;
for(int d=1;d<=n;++d)
for(int k=1,D;(D=k*d)<=n;++k){
int a=n/D,b=m/D;
p1=(p1+x/d*d*a%P*b%P*mu[k]%P)%P;
p3=(p3+(ll)d*a%P*b%P*mu[k]%P)%P;
}
p2=((n*(n-1)>>1)%P)*((m*(m-1)>>1)%P)%P;
ll ans=(p1+(p2+p3)*r2%P)%P;
cout<<ans<<endl;
return 0;
}
