友谊赛 Rainyrabbit 爱求和

这篇博客介绍了如何解决一个关于等比数列求和的问题,并通过公倍数分析将其转化为更易处理的形式。通过莫比乌斯函数和预处理,作者提供了一个时间复杂度优化的解决方案,适用于大规模查询。关键步骤包括利用模运算、莫比乌斯反演和分块维护技巧。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题面描述

已知定值 ccc,给出 a,ba,ba,b 的值,求 ∑i=1a∑j=1b∑k=0c∑d=1idk[ilcm(d,j)]\sum_{i=1}^a \sum_{j=1}^b \sum_{k=0}^c \sum_{d=1}^i d^k [\frac{i}{lcm(d,j)}]i=1aj=1bk=0cd=1idk[lcm(d,j)i]998244353998244353998244353的值。共 ttt 次询问。
1⩽t⩽200000,1⩽a⩽106,1⩽b⩽1018,0⩽c⩽10181 \leqslant t \leqslant 200000,1 \leqslant a \leqslant 10^6,1 \leqslant b \leqslant 10^{18},0 \leqslant c \leqslant 10^{18}1t200000,1a106,1b1018,0c1018
时间限制3s,空间限制512MB。

题解

不难看出原式中的等比数列求和,将原式变为:
∑i=1a∑j=1b∑d=1i[ilcm(d,j)]g(d)\sum_{i=1}^a \sum_{j=1}^b \sum_{d=1}^i [\frac{i}{lcm(d,j)}] g(d)i=1aj=1bd=1i[lcm(d,j)i]g(d)
其中 g(n)=∑i=0cnig(n)=\sum_{i=0}^c n^ig(n)=i=0cni,即:
g(n)={c+1(d=1)dc+1−1d−1(d>1)g(n)=\begin{cases} c+1(d=1) \\ \frac{d^{c+1}-1}{d-1}(d > 1) \end{cases}g(n)={c+1(d=1)d1dc+11(d>1)
考虑到 [ilcm(d,j)][\frac{i}{lcm(d,j)}][lcm(d,j)i] 即为不超过 iii 且同时为 d,jd,jd,j 倍数的数的个数,枚举 d,jd,jd,j 的公倍数,得到:
∑i=1a∑j=1b∑d=1ig(d)∑k=1i[j∣k][d∣k]\sum_{i=1}^a \sum_{j=1}^b \sum_{d=1}^i g(d) \sum_{k=1}^i [j|k][d|k]i=1aj=1bd=1ig(d)k=1i[jk][dk]
交换求和符号顺序:
∑i=1a∑k=1i∑j∣kb∑d∣kig(d)\sum_{i=1}^a \sum_{k=1}^i \sum_{j|k}^b \sum_{d|k}^i g(d)i=1ak=1ijkbdkig(d)
f(n)=∑d∣ng(d)f(n) = \sum_{d|n} g(d)f(n)=dng(d),则有:
∑i=1a∑k=1i∑j∣kbf(i)\sum_{i=1}^a \sum_{k=1}^i \sum_{j|k}^b f(i)i=1ak=1ijkbf(i)
考虑每个 f(i)f(i)f(i) 被计算的次数,得到等价的式子:
∑i=1a∑j∣ibf(i)(a−i+1)\sum_{i=1}^a \sum_{j|i}^b f(i)(a-i+1)i=1ajibf(i)(ai+1)
b2=min⁡(a,b)b_2=\min(a,b)b2=min(a,b),得到最终的式子:
(a+1)∑i=1a∑j∣ib2f(i)−∑i=1a∑j∣ib2f(i)i(a+1) \sum_{i=1}^a \sum^{b_2}_{j|i} f(i) - \sum_{i=1}^a \sum_{j|i}^{b_2} f(i)i(a+1)i=1ajib2f(i)i=1ajib2f(i)i
对询问的 b2b_2b2 排序,用分块维护式子的值即可。预处理 fff 函数的值。
时间复杂度 O(aln⁡a+ta)O(a \ln a + t \sqrt{a})O(alna+ta),空间复杂度 O(a)O(a)O(a)

#include<stdio.h>
#include<algorithm>
#define R register int
#define L long long
#define I inline
#define N 1000001
#define P 998244353
int ans[200000],g[N],inv[N],f[N],c1[1000488],c2[1000488],sum1[977],sum2[977],prime[78498];
I void Add(int&x,int y){
	x+=y;
	if(x>=P){
		x-=P;
	}
}
I int PowMod(int x,int y){
	int res=1;
	while(y!=0){
		if((y&1)==1){
			res=(L)res*x%P;
		}
		y>>=1;
		x=(L)x*x%P;
	}
	return res;
}
struct Query{
	int A,B,Id;
	I friend bool operator<(Query x,Query y){
		return x.B<y.B;
	}
}q[200000];
I void Modify(const int x){
	const int d1=f[x],d2=(L)x*f[x]%P;
	Add(c1[x],d1);
	Add(sum1[x>>10],d1);
	Add(c2[x],d2);
	Add(sum2[x>>10],d2);
}
I int GetAns(int x){
	int res1=0,res2=0,b=x>>10;
	for(R i=0;i!=b;i++){
		Add(res1,sum1[i]);
		Add(res2,sum2[i]);
	}
	for(R i=b<<10;i<=x;i++){
		Add(res1,c1[i]);
		Add(res2,c2[i]);
	}
	return((1ll+x)*res1-res2+P)%P;
}
int main(){
	inv[1]=1;
	int t,a=0;
	scanf("%d",&t);
	L c;
	scanf("%lld",&c);
	c++;
	g[1]=c%P;
	c%=P-1;
	for(R i=2;i!=N;i++){
		inv[i]=(L)(P-P/i)*inv[P%i]%P;
		if(f[i]==0){
			prime[a]=i;
			g[i]=PowMod(i,c);
			a++;
		}
		for(R j=0;prime[j]*i<N;j++){
			int t=i*prime[j];
			f[t]=-1;
			g[t]=(L)g[i]*g[prime[j]]%P;
			if(i%prime[j]==0){
				break;
			}
		}
	}
	for(R i=2;i!=N;i++){
		g[i]=(g[i]-1ll)*inv[i-1]%P;
	}
	for(R i=1;i!=N;i++){
		for(R j=i;j<N;j+=i){
			if(f[j]==-1){
				f[j]=0;
			}
			Add(f[j],g[i]);
		}
	}
	for(R i=0;i!=t;i++){
		q[i].Id=i;
		scanf("%d",&q[i].A);
		L b;
		scanf("%lld",&b);
		q[i].B=q[i].A>b?b:q[i].A;
	}
	std::sort(q,q+t);
	a=1;
	for(R i=0;i!=t;i++){
		while(a<=q[i].B){
			for(R j=a;j<N;j+=a){
				Modify(j);
			}
			a++;
		}
		ans[q[i].Id]=GetAns(q[i].A);
	}
	for(R i=0;i!=t;i++){
		printf("%d\n",ans[i]);
	}
	return 0;
}

此题限于式子结构不可莫比乌斯反演。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值