LG3768 简单的数学题

P3768 简单的数学题

题目描述

输入一个整数n和一个整数p，你需要求出$(\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))~mod~p$，其中gcd(a,b)表示a与b的最大公约数。

输入输出格式

输入格式：

一行两个整数p、n。

输出格式：

一行一个整数$(\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))~mod~p$。

输入输出样例

输入样例#1： 复制

998244353 2000

输出样例#1： 复制

883968974

说明

对于20%的数据，$n \leq 1000$。

对于30%的数据，$n \leq 5000$。

对于60%的数据，$n \leq 10^6$，时限1s。

对于另外20%的数据，$n \leq 10^9$，时限3s。

对于最后20%的数据，$n \leq 10^{10}$，时限6s。

对于100%的数据，$5 \times 10^8 \leq p \leq 1.1 \times 10^9$且p为质数。

题解

先推一波式子。

\[\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\ =\sum_{i=1}^n\sum_{j=1}^nij \sum_{d|i\wedge d|j}\varphi(d)\\ =\sum_{d=1}^n\varphi(d)\sum_{d|i}\sum_{d|j}ij\\ =\sum_{d=1}^n\varphi(d)d^2(\sum_{i=1}^{\lfloor\frac nd\rfloor}i)^2\\ =\sum_{d=1}^n\varphi(d)d^2S(\lfloor\frac nd\rfloor,3) \]

其中\(S\)表示自然数幂和。所以只要能求出\(f(n)=\varphi(n)n^2\)的前缀和\(F(n)=\sum_{i=1}^nf(i)\)，这个式子就能整数分块做。

考虑\(f(n)=\varphi(n)n^2\)它的形式，\(f=\varphi\cdot id^2\)，符合杜教筛处理的形式。于是构造\(g(n)=n^2\)，则

\[(f*g)(n)=\sum_{d|n}\varphi(d)d^2(\frac nd)^2=n^3 \]

而\(g,f*g\)的前缀和就是\(S(n,2),S(n,3)\)，所以这个问题就解决了。最后列一下求和式：

\[F(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)F(\lfloor \frac ni\rfloor)\\ =\sum_{i=1}^ni^3-\sum_{i=2}^ni^2F(\lfloor \frac ni\rfloor) \]

时间复杂度\(O(n^\frac 23+n^\frac 12)\)。话说这题数据真的水，我减号写成加号都有60分。

#include<bits/stdc++.h>
#define il inline
#define co const
template<class T>T read(){
    T data=0,w=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-') w=-w;
    for(;isdigit(ch);ch=getchar()) data=data*10+ch-'0';
    return data*w;
}
template<class T>il T read(T&x) {return x=read<T>();}
typedef long long LL;
using namespace std;

co int N=4641589;
int mod,i6;
int pri[N],tot,phi[N];
int add(int a,int b){
	return (a+=b)>=mod?a-mod:a;
}
int mul(int a,int b){
	return (LL)a*b%mod;
}
int fpow(int a,int b){
	int ans=1;
	for(;b;b>>=1,a=mul(a,a))
		if(b&1) ans=mul(ans,a);
	return ans;
}
void init(){
	i6=fpow(6,mod-2);
	pri[1]=phi[1]=1;
	for(int i=2;i<N;++i){
		if(!pri[i]) pri[++tot]=i,phi[i]=i-1;
		for(int j=1;j<=tot&&i*pri[j]<N;++j){
			pri[i*pri[j]]=1,phi[i*pri[j]]=phi[i]*phi[pri[j]];
			if(i%pri[j]==0){
				phi[i*pri[j]]=phi[i]*pri[j];
				break;
			}
		}
	}
	for(int i=2;i<N;++i) phi[i]=add(phi[i-1],mul(phi[i],mul(i,i)));
}
int sum_s2(int n){
	return mul(i6,mul(n,mul(n+1,2*n+1)));
}
int sum_s3(int n){
	return n&1?mul(mul(n,(n+1)/2),mul(n,(n+1)/2)):mul(mul(n/2,n+1),mul(n/2,n+1));
}
map<LL,int> sf;
int sum_f(LL n){
	if(n<N) return phi[n];
	if(sf.count(n)) return sf[n];
	int ans=sum_s3(n%mod);
	for(LL l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=add(ans,mod-mul(add(sum_s2(r%mod),mod-sum_s2((l-1)%mod)),sum_f(n/l)));
	}
	return sf[n]=ans;
}
int solve(LL n){
	int ans=0;
	for(LL l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=add(ans,mul(add(sum_f(r),mod-sum_f(l-1)),sum_s3(n/l%mod)));
	}
	return ans;
}
int main(){
	cerr<<(sizeof(pri)+sizeof(phi))/1024.0/1024<<endl;
	read(mod);
	init();
	printf("%d\n",solve(read<LL>()));
	return 0;
}