lp3455 POI2007 ZAP-Queries

曾经有一道题,叫做YY的GCD,它求的是这样一个值:
$$\begin{equation}\begin{split}\sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{p\in P}[gcd(x,y)==p] \end{split}\end{equation}$$
然后观察这道题的求和式:
$$\begin{equation}\begin{split}\sum_{x=1}^{n}\sum_{y=1}^{m}[gcd(x,y)==d] \end{split}\end{equation}$$
长得很像对不对?
事实上就是前者的简化版。
剩下的操作就很套路了。
$$\begin{equation}\begin{split}&\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(x,y)==1]\\=&\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{k|gcd(x,y)}\mu(k)\\=&\sum_{x=1}^{\lfloor\frac{n}{dk}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{dk}\rfloor}\mu(k)\\=&\sum_{k=1}^{min(n,m)}\mu(k)\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor \end{split}\end{equation}$$
然后预处理出莫比乌斯函数跑数论分块即可。

#include<iostream>
#include<cstdio>

const int N = 100000+5;

int p[N/10],mu[N],n,m,d;
bool ip[N];
void prpr(){
	p[0]=0;mu[1]=1;ip[1]=1;
	for(int i=2;i<=100000;++i){
		if(!ip[i]){
			p[++p[0]]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=p[0]&&p[j]*i<=100000;++j){
			ip[i*p[j]]=1;
			if(!(i%p[j])){
				mu[i*p[j]]=0;
				break;
			}else{
				mu[i*p[j]]=-mu[i];
			}
		}
	}
	for(int i=2;i<=100000;++i){
		mu[i]+=mu[i-1]; 
	}
}
long long ans;
void init(){
	scanf("%d%d%d",&n,&m,&d);
	n>m?n^=m^=n^=m:0;
	ans=0;
	int k=0;
	for(int i=1;i<=n;i=k+1){
		k=std::min(n/(n/i),m/(m/i));
		ans+=1ll*(n/(i*d))*(m/(i*d))*(mu[k]-mu[i-1]);
	}
	printf("%lld\n",ans);
}

int main(){
	prpr();
	int T;
	scanf("%d",&T);
	while(T--){
		init();
	}
	return 0;
}

lp2257 YY的GCD

不妨设\(n<m\)
首先将原提问形式化,可以得到原式为:
$$\sum_{x=1}^{n}\sum_{y=1}^{m}\sum_{p\in P}[gcd(x,y)==p] $$
将所有的\(x,y\)都除以\(gcd(x,y)\),则原式可转化为:
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}[gcd(x,y)==1]$$
然后,根据莫比乌斯函数的基本性质,我们有:
$$\sum_{d|n}\mu(d)=[n==1]$$
替换入原式可以得到:
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}\sum_{d|gcd(x,y)}\mu(d)$$
将\(d|gcd(x,y)\)一式化简可得:
$$d|gcd(x,y)=d|x\&\&d|y$$
代入原式:
$$\sum_{x=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{p\in P}\sum_{d|x\&\&d|y}\mu(d)$$
将所有的\(x,y\)都除以\(d\),则原式可转化为:
$$\sum_{x=1}^{\lfloor \frac{n}{pd} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{pd} \rfloor}\sum_{p\in P}\mu(d)$$
由于\(x,y\)对后半部分式子无影响,且\(pd\le n\)原式可变式为:
$$\sum_{p\in P}\sum_{d=1}^{\lfloor \frac{n}{p} \rfloor}\mu(d)\lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor$$
为了减少枚举的数量,不妨设\(q=pd\),那么我们又可以将原式转化为:
$$\sum_{q=1}^{n}\lfloor \frac{n}{q} \rfloor \lfloor \frac{m}{q} \rfloor \sum_{p\in P\&\&pd==q}\mu(d)$$
又,根据狄利克雷卷积的基本式
$$g=f u \Rightarrow g(x)=\sum_{a b=x} f(a) u(b)$$
我们发现,上式靠右的如下部分可以作出变换:
$$\sum_{p\in P\&\&pd==q}\mu(d)=g(q)=\sum_{pd=q}[p\in P]\mu(d)$$
若令:
$$f(p)=[p\in P]$$
则我们可以发现:
$$ g(q)=f(p) \mu(d)$$
故而这一部分的值可以用狄利克雷卷积来\(O(nln_n)\)预处理。
一看数据范围,完了:\(10^7\),\(O(nln_n)\)的复杂度显然是不够优秀的。我们需要将它优化到\(O(n)\)

我们考虑对这个函数\(g\)进行线性筛。
我们发现,上式可以进行如下变换:
$$g(q)=\sum_{pd=q\&\&p\in P}\mu(d)=\sum_{p\in P\&\&p|q}\mu(\frac{q}{p})$$
本质上就是枚举其质因子之和并求其除以该质因子的莫比乌斯函数值之和。
我们分三类讨论。
对于\(g(x)\)
如果\(x\in P\),显然\(g(x)=\mu(1)=1\)
如果\(x=a^vb,a\in P,v>1\),则显然枚举其他的质数是没有意义的。\(g(x)=\mu(b)\)
如果\(x=ab,a\in P\),则,\(g(x)\)的值,等价于在计算\(g(b)\)的和的时候每一次计算都多乘上一个\(a\)的贡献,也就是乘上负一。然后将这个值再加上枚举\(a\)的情况,也就是\(\mu(b)\),即\(g(x)=-g(b)+\mu(b)\)
然后大力上一个线性筛即可。

对于左半部分,上一个数论分块即可。

#include<iostream>
#include<cstdio>

const int N = 10000000+5;

int p[N/10],g[N],mu[N],n,m;
bool ip[N];
void prpr(){
	p[0]=0;mu[1]=1;ip[1]=1;
	for(int i=2;i<=10000000;++i){
		if(!ip[i]){
			p[++p[0]]=i;
			mu[i]=-1;
			g[i]=1;
		}
		for(int j=1;j<=p[0]&&p[j]*i<=10000000;++j){
			ip[i*p[j]]=1;
			if(!(i%p[j])){
				g[i*p[j]]=mu[i];
				mu[i*p[j]]=0;
				break;
			}else{
				mu[i*p[j]]=-mu[i];
				g[i*p[j]]=-g[i]+mu[i];
			}
		}
	}
	for(int i=2;i<=10000000;++i){
		g[i]+=g[i-1]; 
	}
}
long long ans;
void init(){
	scanf("%d%d",&n,&m);
	n>m?n^=m^=n^=m:0;
	ans=0;
	int k=0;
	for(int i=1;i<=n;i=k+1){
		k=std::min(n/(n/i),m/(m/i));
		ans+=1ll*(n/i)*(m/i)*(g[k]-g[i-1]);
	}
	printf("%lld\n",ans);
}

int main(){
	prpr();
	int T;
	scanf("%d",&T);
	while(T--){
		init();
	}
	return 0;
}

lp2260 清华集训2012 模积和

题目大意:求式子:
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j) (i \ne j)$$
的值。

那么容斥以后等价于求:
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j)-\sum_{i=1}^{min(n,m)}(n\ mod\ i)*(m\ mod\ i)$$

首先求:
$$\sum_{i=1}^{n}\sum_{j=1}^{m}(n\ mod\ i)*(m\ mod\ j)$$

首先,我们有定理:
$$ n\ mod\ p=n-p*\lfloor \frac{n}{p} \rfloor$$
然后配合上\(\Sigma\)的运算法则,我们可以把原式转化为:
$$n^{2}*m^{2}-m^{2}\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor-n^{2}\sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor-\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor \sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor$$
我们设:
$$ P=\sum_{i=1}^{n}i\lfloor \frac{n}{i} \rfloor,Q=\sum_{j=1}^{m}j\lfloor \frac{m}{j} \rfloor$$
所以我们可以用数论分块来分别计算\(P,Q\)的值。

而,对于第二部分:
$$R=\sum_{i=1}^{min(n,m)}(n\ mod\ i)*(m\ mod\ i)$$
我们不妨令\(n<m\)
那么将原式展开,可得:
$$R=\sum_{i=1}^{n}n*m-n*i\lfloor \frac{m}{i} \rfloor-m*i\lfloor \frac{n}{i} \rfloor+i^{2}\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor$$
这个式子显然也是可以数论分块的。但是我们需要计算\(i^{2}\)对于式子的贡献。
用一些简单的证明方法可以得到:
$$\sum_{i=1}^{n}i^{2}=\frac{n(2n+1)(n+1)}{6}$$
那么它的区间和即是:
$$\sum_{i=1}^{r}i^{2}-\sum_{i=1}^{l-1}i^{2}=\frac{(2r^2+2rl+2l^2-r-3l+2)(r-l+1)}{6}$$
于是我们就可以统计整个结果了。

现在我们来考虑数论分块。

数论分块是一种处理形如
$$ \sum_{i=1}^{n}\lfloor \frac{n}{i} \rfloor$$
的公式的求值的方法。

而对于函数:

$$f(x)=\lfloor \frac{a}{x} \rfloor$$

可以发现,对于连续的一段\(x\)的区间,它的值是相同的。
具体来说是,令值为\(k\)则,令区间的左端点为\(l\),则整个区间的值为:
$$k=\lfloor \frac{n}{l} \rfloor$$
那么右端点\(r\)很显然是满足方程:
$$\lfloor \frac{n}{x} \rfloor=k$$
的最大值。
依据整除的性质得到的引理一,我们发现:
$$r=\lfloor \frac{n}{k} \rfloor$$
从而,区间为:
$$[l,\lfloor \frac{n}{\lfloor \frac{n}{l} \rfloor} \rfloor$$

附:对于引理一的证明:

引理一:求证方程:
$$\lfloor \frac{n}{x} \rfloor=k$$
的解的最大值\(t\)为
$$t=\lfloor \frac{n}{k} \rfloor$$
我们令方程的一个合法解\(x\)满足:
$$x*k+r=n\ (0\le r < x)\ (1)$$
并令:
$$t=x+d$$
那么,有:
$$\lfloor \frac{n}{x+d} \rfloor=k$$
则:
$$(x+d)*k+r’=n\ (0\le r < x+d)\ (2)$$
联立方程\((1),(2)\)可得:
$$d*k+r’=r$$
由整除的定义可得:
$$d=\lfloor \frac{r}{k} \rfloor$$
很显然,当
$$x=\lfloor \frac{n}{k} \rfloor$$
时,
$$r=n\ mod\ k,r<k,\lfloor \frac{r}{k} \rfloor=0$$
此时\(x\)取到最大值。
证毕。

另:尽量用模块化编程,而不是去化简式子。这样可以用复杂度换取正确率。

#include<iostream>
#include<cstdio>
const long long MOD = 19940417;
const long long inv2 = 9970209;
const long long inv6 = 3323403;
long long n,m;
inline long long sm1(long long A){
    return (A*(A+1)%MOD)*inv2%MOD;
}
inline long long sm2(long long A){
    return ((A*(A+1)%MOD)*(2*A+1)%MOD)*inv6%MOD;
}
void init(){
    scanf("%lld%lld",&n,&m);
    n>m?n^=m^=n^=m:0; 
    long long P,Q,R;
    int l,r;
    r=0,P=0;
    while(r<n){
        l=r+1;
        r=(n/l)?(n/(n/l)):n;
        P+=(((r-l+1)*(n/l)%MOD)*(l+r)%MOD*inv2)%MOD;
        P%=MOD;
//        printf("P:%lld\n",P);
    }
    P%=MOD;
    r=0,Q=0;
    while(r<m){
        l=r+1;
        r=(m/l)?(m/(m/l)):m;
        Q+=(((r-l+1)*(m/l)%MOD)*(l+r)%MOD*inv2)%MOD;
        Q%=MOD;
//        printf("Q:%lld\n",Q);
    }
//    printf("%lld %lld\n",P,Q);
    Q%=MOD;
    R=(n*m%MOD)*n%MOD;
    long long ans=(n*m%MOD)*(n*m%MOD)%MOD-(n*n%MOD)*Q%MOD-(m*m%MOD)*P%MOD+P*Q%MOD;
    ans%=MOD;
    ans+=MOD;
    ans%=MOD;
    r=0;
    while(r<n){
        l=r+1;
        r=std::min((n/l)?(n/(n/l)):n,(m/l)?(m/(m/l)):m);
        CNT1=m/l;
        CNT2=n/l;
        R-=((((sm1(r)-sm1(l-1))%MOD)*(m/l)%MOD)*n%MOD+(((sm1(r)-sm1(l-1))%MOD)*(n/l)%MOD)*m%MOD);
        R+=(((sm2(r)-sm2(l-1))%MOD)*((m/l)*(n/l)%MOD))%MOD;
        R%=MOD;
    }
    ans-=R;
    ans%=MOD;
    ans+=MOD;
    ans%=MOD;
    printf("%lld",ans);
    
}
int main(){
    init();
    return 0;
}