ShinriiTin's Blog - 博主是sb

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]

           = ki∈[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<ii∈[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;
}