【题解】 [EZEC-4]求和

对于百分之十的数据：随便过。

下面推式子：

\sum_{i = 1}^{n} \sum_{j = 1}^{n} gcd (i, j)^{i + j}

$\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)^{i+j}$

= \sum_{d = 1}^{n} \sum_{i = 1}^{n} \sum_{j = 1}^{n} d^{i + j} [gcd (i, j) = d]

$=\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^nd^{i+j}[\gcd(i,j)=d]$

= \sum_{d = 1}^{n} \sum_{i = 1}^{\frac{n}{d}} \sum_{j = 1}^{\frac{n}{d}} d^{d (i + j)} [gcd (i, j) = 1]

$=\sum_{d=1}^n\sum_{i=1}^\frac{n}{d}\sum_{j=1}^\frac{n}{d}d^{d(i+j)}[\gcd(i,j)=1]$

= \sum_{d = 1}^{n} \sum_{i = 1}^{\frac{n}{d}} \sum_{j = 1}^{\frac{n}{d}} d^{d (i + j)} \sum_{k | gcd (i, j)} μ (k)

$=\sum_{d=1}^n\sum_{i=1}^\frac{n}{d}\sum_{j=1}^\frac{n}{d}d^{d(i+j)}\sum_{k|\gcd(i,j)}\mu(k)$

= \sum_{d = 1}^{n} \sum_{k = 1}^{\frac{n}{d}} μ (k) \sum_{i = 1}^{\frac{n}{k d}} \sum_{j = 1}^{\frac{n}{k d}} d^{k d (i + j)}

$=\sum_{d=1}^n\sum_{k=1}^\frac{n}{d}\mu(k)\sum_{i=1}^\frac{n}{kd}\sum_{j=1}^\frac{n}{kd}d^{kd(i+j)}$

令 $T=kd:$

= \sum_{T = 1}^{n} \sum_{k | T} μ (k) \sum_{i = 1}^{\frac{n}{T}} \sum_{j = 1}^{\frac{n}{T}} [(\frac{T}{k})^{T}]^{i + j}

$=\sum_{T=1}^n\sum_{k|T}\mu(k)\sum_{i=1}^\frac{n}{T}\sum_{j=1}^\frac{n}{T}[(\frac{T}{k})^T]^{i+j}$

现在的问题在于 $\sum_{i=1}^\frac{n}{T}\sum_{j=1}^\frac{n}{T}[(\frac{T}{k})^T]^{i+j}.$

线性递推

以下是@SOSCHINA大佬的思路：

设 $g(n)=\sum_{i=1}^n\sum_{j=1}^n k^s.$

枚举 $s=i+j.$

则有：

g (n) = \sum_{s = 2}^{n + 1} (s - 1) k^{s} + \sum_{s = n + 2}^{2 n} (2 n - s + 1) k^{s}

$g(n)=\sum_{s=2}^{n+1}(s-1)k^s+\sum_{s=n+2}^{2n}(2n-s+1)k^s$

g (n + 1) = \sum_{s = 2}^{n + 1} (s - 1) k^{s} + \sum_{s = n + 2}^{2 n + 2} (s n + 3 - s) k^{s}

$g(n+1)=\sum_{s=2}^{n+1}(s-1)k^s+\sum_{s=n+2}^{2n+2}(sn+3-s)k^s$

g (n + 1) - g (n) = \sum_{s = n + 2}^{2 n} 2 k^{s} + 2 k^{2 n + 1} + k^{2 n + 2}

$g(n+1)-g(n)=\sum_{s=n+2}^{2n}2k^s+2k^{2n+1}+k^{2n+2}$

第三行就是两行相减。

对第一行的解释： $[2,n+1]$ 这里的数，每个数作为 $i+j$ 都出现了 $x-1$ 次。因为 $i$ 可以取遍 $[1,x-1].$ 后面的那一些, $[n+2,2n]$ 会发现 $i$ 最大只能到 $n,$ 不能再取遍 $x-1$ 个值了。此时能取到的应该是 $2n-s+1$ 种。

对于 $g(n+1):$ 这里是把第一个式子的最后一个值移动到了后面那个式子，方便做差。

这时我们可以在小模数的情况下做到 $O(n*mod$ )的预处理。

化简形式

令 $x=(\frac{T}{k})^T.$

则原式为 $\sum_{i=1}^\frac{n}{T}\sum_{j=1}^\frac{n}{T} x^{i+j}.$

像不像一个多项式。

它就等于 $(x+x^2+...x^\frac{n}{T})^2.$

于是我们可以等比数列求和解出。

剩下的，可以做到 $O(n\log n\log mod)$ 处理出整个式子。

#include<bits/stdc++.h>
using namespace std;
const int MAXN=1500001;
int mod,TT;
bitset<MAXN+1>vis;
int p[MAXN+1],mu[MAXN+1],T[MAXN+1],cnt,n,Ans;
inline int Mod(long long x){
	if(x<0)return x+mod;
    if(x>=mod)return x%mod;
    return x;
}
inline int add(int x,int y) {return Mod(1ll*x+1ll*y+1ll*mod);}
inline int mul(int x,int y) {return Mod(1ll*x*y);}
inline int qpow(int a,int b) {
	if(!b)return 1;
	if(a<=1||b==1)return a;
	a %= mod; 
	int res=1;
	while(b) {
		if(b&1)res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
inline int calc(int x,int y){
	if(y==1)return x;
	if(x==1)return y;
	int ans=x;
	int inv=qpow((1-x+mod)%mod,mod-2);
	int fm=(1-qpow(x,y)+mod)%mod;
	ans=mul(ans,mul(fm,inv));
	return ans;
}
inline int Calc(int x,int y){int ans=calc(x,y);return mul(ans,ans);} 
int main() {
	scanf("%d",&TT);
	mu[1]=1;
	int N=MAXN;
	for(register int i=2; i<=N; ++i) {
		if(!vis[i])p[++cnt]=i,mu[i]=-1;
		for(register int j=1; j<=cnt&&i*p[j]<=N; ++j) {
			vis[i*p[j]]=1;
			if(i%p[j]==0)break;
			mu[i*p[j]]=-mu[i];
		}
	}
	while(TT--) {
		scanf("%d%d",&n,&mod);
		N=n;Ans=0;
		for(register int i=1; i<=N; ++i) {
			for(register int j=i,k,x; j<=N; j+=i) {
				k=i;if(!mu[k])continue;
				x=qpow(j/k,j);
				T[j]=add(T[j],mul(mu[k],Calc(x,n/j)));
			}
		}
		for(register int i=1; i<=n; ++i)Ans=add(Ans,T[i]),T[i]=0;
		printf("%d\n",Ans);
	}

	return 0;
}

由于这里是 $5*10^5$ 的数据，所以略微卡常，但笔者通过非常不精湛的卡常技术跑到了 $3s$ 以内，所以这里的时间限制我开了 $3.2s$ .

对等比数列进行精细处理，可以做到 $O(n\log^2n)$ 的复杂度。

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int MAXN=1500000;
int mod,TT;
bitset<MAXN<<1>vis;
int p[MAXN<<1],mu[MAXN<<1],T[MAXN<<1],cnt,n,Ans;
inline int Mod(long long x){
    if(x>=mod)return x%mod;
    return x;
}
inline int add(int x,int y) {
	return Mod(x+y+mod);
}
inline int mul(int x,int y) {
	return Mod(1ll*x*y);
}
inline int qpow(int a,int b) {
	if(!b)return 1;
	if(a<=1||b==1)return a;
	a %= mod; 
	int res=1;
	while(b) {
		if(b&1)res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
inline int calc(int x,int y){
	if(y==1)return x;
		int res=calc(x,y/2);
	res=add(res,mul(res,qpow(x,y/2)));
	if(y&1)res=add(res,mul(x,qpow(x,y-1)));
	return res;
}
inline int Calc(int x,int y){int ans=calc(x,y);return mul(ans,ans);} 
signed main() {
	scanf("%lld",&TT);
	mu[1]=1;
	int N=MAXN;
	for(int i=2; i<=N; ++i) {
		if(!vis[i])p[++cnt]=i,mu[i]=-1;
		for(int j=1; j<=cnt&&i*p[j]<=N; ++j) {
			vis[i*p[j]]=1;
			if(i%p[j]==0)break;
			mu[i*p[j]]=-mu[i];
		}
	}
	while(TT--) {
		scanf("%lld%lld",&n,&mod);
		N=n;
		Ans=0;
		for(int i=1; i<=N; ++i) {
			for(int j=i; j<=N; j+=i) {
				int k=i;
				int x=qpow(j/k,j);
				if(!mu[k])continue;
				T[j]=add(T[j],mul(mu[k],Calc(x,n/j)));
			}
		}
		for(int i=1; i<=n; ++i)Ans=add(Ans,T[i]),T[i]=0;
		cout<<Ans<<endl;
	}

	return 0;
}

由于常数等原因，这分代码可以拿到 $50$ 分的好成绩。但我们可以通过另一种做法将常数/复杂度降低。

另一种做法

观察：

\sum_{d = 1}^{n} \sum_{k = 1}^{\frac{n}{d}} μ (k) \sum_{i = 1}^{\frac{n}{k d}} \sum_{j = 1}^{\frac{n}{k d}} d^{k d (i + j)}

$\sum_{d=1}^n\sum_{k=1}^\frac{n}{d}\mu(k)\sum_{i=1}^\frac{n}{kd}\sum_{j=1}^\frac{n}{kd}d^{kd(i+j)}$

= \sum_{d = 1}^{n} \sum_{k = 1}^{\frac{n}{d}} μ (k) (d^{k d} + d^{2 k d} + . . . + d^{\frac{n}{k d} * k d = n})^{2} .

$=\sum_{d=1}^n \sum_{k=1}^\frac{n}{d} \mu(k)(d^{kd}+d^{2kd}+...+d^{\frac{n}{kd}*kd=n})^2.$

这里同样观察式子发现可以直接算。前一部分是 $O(n\ln n)$ 的 $n$ 倍调和级数的复杂度，后面带上一个 $O(\log n)$ 精细处理的等比数列求求和复杂度。

（代码中的优化即使不加也是可以过的）

#define __AVX__ 1
#define __AVX2__ 1
#define __SSE__ 1
#define __SSE2__ 1
#define __SSE2_MATH__ 1
#define __SSE3__ 1
#define __SSE4_1__ 1
#define __SSE4_2__ 1
#define __SSE_MATH__ 1
#define __SSSE3__ 1
#pragma GCC optimize("Ofast,no-stack-protector,unroll-loops,fast-math")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4.1,sse4.2,avx,avx2,popcnt,tune=native")
#include <immintrin.h>
#include <emmintrin.h>
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <string>
#include <bitset>
using namespace std;
const int MAXN=1.5e6+10;
int mod,T;
bitset<MAXN+1>vis;
int p[MAXN+1],cnt,mu[MAXN+1],N;
inline int Mod(long long a, int pp){
    return a>=pp ? a%pp : a>=0 ? a : a+pp;
}
inline int add(int x,int y){return Mod( (1ll+x+y+mod-1ll),mod);}
inline int mul(int x,int y){return Mod(1ll*x*y,mod);}
void pretreatment(){
	mu[1]=1;
	for(int i=2;i<=MAXN;++i){
		if(!vis[i])p[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*p[j]<=MAXN;++j){
			vis[i*p[j]]=1;
			if(Mod(i,p[j])==0)break;
			mu[i*p[j]]=-mu[i];
		}
	}
}
inline int qpow(int a,int b){
	if(!b)return 1;
	if(a<=1||b==1)return a;
	int res=1;
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);b>>=1;
	}
	return res;
}
inline int calc(int x,int y){
	if(y==1)return x;
	int res=calc(x,y>>1);
	res=add(res,mul(res,qpow(x,y>>1)));
	if(y&1)res=add(res,mul(x,qpow(x,y-1)));
	return res;
}

inline int Calc(int x,int y){int ans=calc(x,y);return mul(ans,ans);} 
int ssolve(int n,int d){
	int res=0;
	for(register int l=1;l<=n;++l){
		if(!mu[l])continue;
		res=add(res,mul(mu[l],Calc(qpow(d,l),n/l)));
	}
	return res;
}
int solve(int n){
	int ans=0;
	for(register int l=1;l<=n;l++){
		ans=add(ans,ssolve(n/l,qpow(l,l)));
	}
	return ans;
}
signed main(){
	scanf("%lld",&T);
	pretreatment();
	for(;T;T--){
		scanf("%lld%lld",&N,&mod);
		printf("%lld\n",solve(N));
	}
	return 0;
}

可以用整除分块减少循环中乘法的使用，对代码速度可能有一定的提升。

#include<bits/stdc++.h>
using namespace std;
const int MAXN=1.5e6+10;
int mod,T;
bitset<MAXN+1>vis;
int p[MAXN+1],cnt,mu[MAXN+1],N;
inline int Mod(long long a, int pp){return a>=pp ? a%pp : a>=0 ? a : a+pp;}
inline int add(int x,int y){return Mod( (1ll+x+y+mod-1ll),mod);}
inline int mul(int x,int y){return Mod(1ll*x*y,mod);}
inline int qpow(int a,int b){
	if(!b)return 1;
	if(a<=1||b==1)return a;
	a=Mod(a,mod);
	int res=1;
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);b>>=1;
	}
	return res;
}
inline int calc(int x,int y){
	if(y==1)return x;
	int res=calc(x,y>>1);
	res=add(res,mul(res,qpow(x,y>>1)));
	if(y&1)res=add(res,mul(x,qpow(x,y-1)));
	return res;
}
inline int Calc(int x,int y){int ans=calc(x,y);return mul(ans,ans);} 
int ssolve(int n,int d){
	int res=0;
	for(register int l=1,r;l<=n;l=r+1){
		r=(n/(n/l));
		int D=n/l;
		for(int i=l;i<=r;++i){
			if(!mu[i])continue;
			res=add(res,mul(mu[i],Calc(qpow(d,i),D)));
		}
	}
	return res;
}
int solve(int n){
	int ans=0;
	for(register int l=1,r;l<=n;l=r+1){
		r=(n/(n/l));
		int D=n/l;
		for(int i=l;i<=r;++i)ans=add(ans,ssolve(D,qpow(i,i)));
	}
	return ans;
}
int main(){
	scanf("%d",&T);
	mu[1]=1;
	for(register int i=2;i<=MAXN;++i){
		if(!vis[i])p[++cnt]=i,mu[i]=-1;
		for(register int j=1;j<=cnt&&i*p[j]<=MAXN;++j){
			vis[i*p[j]]=1;
			if(Mod(i,p[j])==0)break;
			mu[i*p[j]]=-mu[i];
		}
	}
	for(;T;T--){
		scanf("%d%d",&N,&mod);
		printf("%d\n",solve(N));
	}
	return 0;
}

经由大佬 @C3H5ClO 大佬证明，上面这份代码实际上是 $O(n\log n)$ 的。

这里借用一下 @C3H5ClO 大佬的证明：

f (n) = \sum_{d = 1}^{n} \sum_{e = 1}^{⌊ \frac{n}{d} ⌋} μ (e) (\sum_{i = 1}^{⌊ \frac{n}{d e} ⌋} d^{d e})^{2}

$f(n)=\sum_{d=1}^n\sum_{e=1}^{\lfloor \frac{n}{d} \rfloor}\mu(e)(\sum_{i=1}^{\lfloor \frac{n}{de}\rfloor}d^{de})^2$

T (n) = O (\sum_{d = 1}^{n} \sum_{e = 1}^{⌊ \frac{n}{d} ⌋} \log ⌊ \frac{n}{d e} ⌋)

$T(n)=O(\sum_{d=1}^n\sum_{e=1}^{\lfloor \frac{n}{d} \rfloor}\log\lfloor \frac{n}{de}\rfloor)$

T (n) = O (\sum_{i = 1}^{n} \int_{0}^{\frac{n}{i}} \log \frac{n}{i x} d x)

$T(n)=O(\sum_{i=1}^n\int_0^{\frac{n}{i}}\log\frac{n}{ix}\mathrm{d}x)$

$\int\log\frac{n}{x}\mathrm{d}x=(\log n+1)x-x\log x+C$
$\int_0^n\log\frac{n}{x}\mathrm{d}x=n$

T (n) = O (\sum_{i = 1}^{n} \frac{n}{i}) = O (n \ln n)

$T(n)=O(\sum_{i=1}^n\frac{n}{i})=O(n\ln n)$

（蒟蒻不会微积分惨被教育.jpg）

出这题的本意其实是想看看有没有吊打 $\text{std}$ 的做法的，笔者推了很久并没有找到线性的做法。

posted @ 2020-09-06 18:41 Refined_heart 阅读(604) 评论(3) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· TypeScript + Deepseek 打造卜卦网站：技术与玄学的结合
· Manus的开源复刻OpenManus初探
· .NET Core 中如何实现缓存的预热？
· 三行代码完成国际化适配，妙~啊~
· 阿里巴巴 QwQ-32B真的超越了 DeepSeek R-1吗？

公告

昵称： Refined_heart
园龄： 5年11个月
粉丝： 13
关注： 6

+加关注

2025年3月

日

一

二

三

四

五

六

Refined_heart

人生如棋，落子无悔。

【题解】 [EZEC-4]求和

另一种做法

公告

搜索

常用链接

最新随笔

我的标签

积分与排名

随笔档案

友链

阅读排行榜

评论排行榜

推荐排行榜

最新评论