莫比乌斯反演板子求助
查看原帖
莫比乌斯反演板子求助
237530
rzh123楼主2023/8/23 16:46
#include "cstdio"
#include "string"
#include "algorithm"
using namespace std;
constexpr int N=1e7+5,P{20101009};
int n,m;
int mu[N],mtts[N];
basic_string<int> prime;
void sieve(int n){
	static int d[N];
	mu[1]=1,d[1]=1;
	for(int i{2};i<=n;++i){
		if(!d[i]){
			d[i]=i,mu[i]=-1;
			prime+=i;
		}
		for(int j:prime){
			int t{i*j};
			if(t>n||d[i]<j) break;
			d[t]=j;
			if(i%j==0) break;
			mu[t]=-mu[i];
		}
	}
	for(int t{1};t<=n;++t){
		// printf("mu[%d]=%d\n",t,mu[t]);
		mtts[t]=(mtts[t-1]+1LL*((mu[t]+P)%P)*t%P*t%P)%P;
	}
}
int vsum(int l,int r){
	return 1LL*(l+r)*(r-l+1ll)/2ll%P;
}
signed main(){
	scanf("%d%d",&n,&m);
	sieve(max(n,m));
	int ans{0},mi{min(n,m)};
	for(int l1{1},r1;l1<=mi;l1=r1+1){
		r1=(mi/(mi/l1));
		int sum{0};
		for(int l2{1},r2;l2<=mi/l1;l2=r2+1){
			r2=min((n/l1)/((n/l1)/l2),(m/l1)/((m/l1)/l2));
			sum=(sum+1LL*((mtts[r2]-mtts[l2-1]+P)%P)*vsum(1,n/l1/l2)%P*vsum(1,m/l1/l2)%P)%P;
		}
		ans=(ans+1LL*vsum(l1,r1)*sum%P)%P;
	}
	printf("%d\n",ans);
	return 0;
}

/*
sum(i=1..n)sum(j=1..m)lcm(i,j)
sum(i=1..n)sum(j=1..m)ij/gcd(i,j)
sum(d=1..min)sum(i=1..n)sum(j=1..m)[gcd(i,j)=d]*ij/d
sum(d=1..min)sum(i=1..n/d)sum(j=1..m/d)[gcd(i,j)=1]*i*dj*d/d
sum(d=1..min)d*sum(i=1..n/d)sum(j=1..m/d)[gcd(i,j)=1]*i*j
sum(d=1..min)d*sum(i=1..n/d)sum(j=1..m/d)sum(t|gcd(i,j))mu(t)*i*j

sum(d=1..min)d*sum(i=1..n/d)sum(j=1..m/d)sum(t|gcd(i,j))mu(t)*i*j
sum(d=1..min)d*sum(t=1..min/d)mu(t)*sum(i=1..n/td)sum(j=1..m/td)i*j*t^2
sum(d=1..min)d*sum(t=1..min/d)mu(t)*t^2*sum(i=1..n/td)sum(j=1..m/td)i*j
sum(d=1..min)d*sum(t=1..min/d)mu(t)*t^2*S(1,n/td)*S(1,m/td)
sum(d=1..min)d*sum(t=1..min/d)MD(t)*S(1,n/td)*S(1,m/td)

*/
2023/8/23 16:46
加载中...