Solution -「LOCAL 28731」「重庆市 2021 中学友谊赛」Rainyrabbit 爱求和

Description

Link.

(operatorname{Rainyrabbit}) 是一个数学极好的萌妹子，近期他发现了一个可爱的函数：

[f(n,m,k)=sum_{d=1}^n d^klfloordfrac{n}{operatorname{lcm}(d,m)} floor ]

其中 (operatorname{lcm}(d,m)) 表示 (d) 和 (m) 的最小公倍数。

她觉得只算这一个函数太单调了于是想要对这个函数求和，给出 (a,b,c)，求：

[sum_{i=1}^asum_{j=1}^bsum_{k=0}^cf(i,j,k) ]

同时 (operatorname{I.w.rabbit}) 固定 (c) 的值不变，给出 (T) 组 (a,b)，请回答此时式子对 (998244353) 取模后的值。

算是概括过了额。

Solution

先看 (f) 怎么把那个烦死的整除去掉。

[egin{aligned} f(n,m,k)&=sum_{d=1}^{n}d^{k}lfloorfrac{n}{ ext{lcm}(d,m)} floor \ &=sum_{d=1}^{n}d^{k}lfloorfrac{frac{n}{d}}{frac{m}{gcd(d,m)}} floor \ &=sum_{d=1}^{n}d^{k}lfloorfrac{frac{n}{d}}{frac{m}{gcd(d,m)}} floor \ end{aligned} ]

来看 (lfloorfrac{frac{n}{d}}{frac{m}{gcd(d,m)}} floor)，因为 (lfloorfrac{a}{b} floor=sum_{i=1}^{a}[b|i])，所以 (lfloorfrac{frac{n}{d}}{frac{m}{gcd(d,m)}} floor=sum_{i=1}^{lfloorfrac{n}{d} floor}[frac{m}{gcd(d,m)}|i]=sum_{i=1}^{lfloorfrac{n}{d} floor}[m|id])。所以

[egin{aligned} f(n,m,k)&=sum_{d=1}^{n}d^{k}lfloorfrac{frac{n}{d}}{frac{m}{gcd(d,m)}} floor \ &=sum_{d=1}^{n}d^{k}sum_{i=1}^{lfloorfrac{n}{d} floor}[m|id] \ &=sum_{m|i}^{n}sum_{d|i}d^{k} \ &=sum_{m|i}^{n}sigma_{k}(i) \ end{aligned} ]

然后代回原式。

[egin{aligned} sum_{i=1}^{a}sum_{j=1}^{b}sum_{k=0}^{c}sum_{j|d}^{i}sigma_{k}(d)&=sum_{k=0}^{c}sum_{d=1}^{a}sigma_{k}(d)(a-d+1)sum_{j=1}^{b}[j|d] \ end{aligned} ]

考虑如何计算 (sum_{k=0}^{c}sigma_{k}(d))。把函数又拆回去：

[egin{aligned} sum_{k=0}^{c}sigma_{k}(d)&=sum_{w|d}sum_{k=0}^{c}w^{k}=sum_{w|d}frac{w^{c+1}-1}{w-1} end{aligned} ]

最后一步是等比数列求和，然后你就可以调和级数预处理了。具体来说就是线筛的时候筛一下 (w^{c+1})，这东西是个完全积性函数，你乱筛就行了。

设这玩意儿为 (s(d)=sum_{w|d}frac{w^{c+1}-1}{w-1})，原式改写为：

[sum_{d=1}^{a}s(d)(a-d+1)sum_{j=1}^{b}[j|d] \ ]

然后后面那个 sigma 你也可以反过来直接调和级数。还有就是 (b>a) 的时候没有贡献，所以可以取个 (min)，这样能多几分。

来看看怎么屮多测。

和 数表 那道题一样，我们把询问离线下来，以 (b) 为关键字排序后树状数组。

把中间那个系数拆出来，变成：

[sum_{d=1}^{a}(a+1)s(d)-d imes s(d)sum_{j=1}^{b}[j|d] \ ]

前面那个好说，直接来；后面就在树状数组修改时乘上系数即可。

综上，维护两个树状数组即可。

#include<cstdio>
#include<algorithm>
using namespace std;
void read(long long &x)
{
	x=0;
	char c=getchar();
	while(c<'0'||c>'9')	c=getchar();
	while(c>='0'&&c<='9')
	{
		x=(x<<3)+(x<<1)+(c^'0');
		c=getchar();
	}
}
void write(long long x)
{
	if(x>9)	write(x/10);
	putchar((x%10)^'0');
}
const long long mod=998244353;
struct node
{
	long long a,b,ID;
}nodes[200010];
struct fenwick
{
	#define lowbit(x) ((x)&-(x))
	long long fen[1000010],mx;
	void ins(long long x,long long y)
	{
		while(x<=mx)
		{
			fen[x]=(fen[x]+y)%mod;
			x+=lowbit(x);
		}
	}
	long long find(long long x)
	{
		long long res=0;
		while(x)
		{
			res=(res+fen[x])%mod;
			x^=lowbit(x);
		}
		return res;
	}
}onefe,anofe;
long long t,a,b,c,tag[1000010],prime[1000010],cnt,fu[1000010],exfu[1000010],power[1000010],cur=1,mx,ans[200010];
bool cmp(node one,node ano)
{
	return one.b<ano.b;
}
long long cqpow(long long bas,long long fur)
{
	long long res=1;
	while(fur)
	{
		if(fur&1)	res=res*bas%mod;
		bas=bas*bas%mod;
		fur>>=1;
	}
	return res;
}
void search(long long x)
{
	tag[1]=power[1]=1;
	for(long long i=2;i<=x;++i)
	{
		if(!tag[i])
		{
			prime[++cnt]=i;
			power[i]=cqpow(i,c+1);
		}
		for(long long j=1;j<=cnt&&prime[j]*i<=x;++j)
		{
			tag[prime[j]*i]=1;
			power[prime[j]*i]=power[prime[j]]*power[i]%mod;
			if(i%prime[j]==0)	break;
		}
	}
	fu[1]=(c%mod+1)%mod;
	for(long long i=2;i<=x;++i)	fu[i]=((power[i]-1+mod)%mod)*cqpow(i-1,mod-2)%mod;
	for(long long i=1;i<=x;++i)
	{
		for(long long j=i;j<=x;j+=i)	exfu[j]=(exfu[j]+fu[i])%mod;
	}
}
int main()
{
	read(t);
	read(c);
	search(1000000);
	for(long long i=1;i<=t;++i)
	{
		read(nodes[i].a);
		read(nodes[i].b);
		nodes[i].ID=i;
		nodes[i].b=min(nodes[i].a,nodes[i].b);
		mx=max(mx,nodes[i].a);
	}
	sort(nodes+1,nodes+t+1,cmp);
	onefe.mx=anofe.mx=mx;
	for(long long i=1;i<=mx&&cur<=t;++i)
	{
		for(long long j=i;j<=mx;j+=i)
		{
			onefe.ins(j,exfu[j]);
			anofe.ins(j,exfu[j]*j%mod);
		}
		while(i==nodes[cur].b)
		{
			ans[nodes[cur].ID]=((onefe.find(nodes[cur].a)*(nodes[cur].a+1))%mod-anofe.find(nodes[cur].a)+mod)%mod;
			cur++;
		}
	}
	for(long long i=1;i<=t;++i)
	{
		write(ans[i]);
		putchar('
');
	}
	return 0;
}