ProjectEuler-Problem 625. Gcd Sum -题解

【题目地址】

题意简述

求当$n=10^{11}$时的如下式子在$\mathrm{m}\mathrm{o}\mathrm{d}998244353$后的值。

$$\sum _{i=1}^n\sum _{j=1}^{i}gcd(i,j)$$

---

我们可以枚举$gcd$，将式子转化为：

$$\begin{gathered} \sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^i[gcd(i,j)=d] \\ \sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^i[gcd(i,j)=1] \end{gathered}$$

根据定义可知，后面的${\sum }_{j=1}^i[gcd(i,j)=1]=\phi (i)$

所以原式可以继续化简为:

$$\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\phi (i)$$

然后我们对于后面的$\phi$的前缀和用杜教筛，前面分块就可以在$O(n^{\frac{2}{3}})$时间内解决，大概$n=10^{11}$的时候大约5s左右可以跑出来。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#define ll long long
using namespace std;
const ll M=3e7+10;
const ll Mod=998244353ll;
ll phi[M],prime[M],cnt,inv_2;
bool vis[M];
map<ll,ll> mp;
ll fpow(ll a,ll b){
	ll ans=1ll;
	for(;b;b>>=1ll,a=(a*a)%Mod){
		if(b&1ll)ans=(a*ans)%Mod;
	}
	return ans;
}
void init(){
	phi[1ll]=1ll;
	inv_2=fpow(2ll,Mod-2ll);
	for(ll i=2ll;i<M;i++){
		if(!vis[i]){prime[++cnt]=i;phi[i]=i-1ll;}
		for(ll j=1,v;j<=cnt&&i*prime[j]<M;j++){
			v=i*prime[j];
			vis[v]=1;
			if(!(i%prime[j])){
				phi[v]=phi[i]*prime[j];
				break;
			}
			phi[v]=phi[i]*phi[prime[j]];
		}
	}
	for(ll i=2ll;i<M;i++)(phi[i]+=phi[i-1ll])%=Mod;
}
#include<cmath>
ll S(ll a){if(a>=Mod)a%=Mod;return ((a*(a+1))%Mod)*inv_2%Mod;}
ll calc(ll x){
	if(x<M) return phi[x];
	if(mp.count(x)) return mp[x];
	ll ans=S(x);
	for(ll i=2ll,j;i<=x;i=j+1ll){
		j=x/(x/i);
		ans=((ans-((j-i+1)%Mod*calc(x/i)%Mod))%Mod+Mod)%Mod;
	}
	return mp[x]=ans;
}
ll solve(ll x){
	ll ans=0,ls;
	for(ll i=1,j;i<=x;i=j+1ll){
		j=(x/(x/i));
		ls=(S(j)-S(i-1ll))%Mod;
		if(ls<0)ls+=Mod;
		ans=(ans+ls*(calc(x/i))%Mod)%Mod;
	}
	return (ans)%Mod;
}
int main(){
	init();
	ll x=1e11;
	printf("%lld\n",solve(x));
	return 0;
}