Sciorz'Blog

hihocoder1456 Rikka with Lattice

  • Dec 17, 2018

有一个 $n × m$ 的点阵，他想要从这 $n × m$ 个点中选出三个点 ${A, B, C}$，满足 1.三角形 $ABC$ 面积不为0且其内部不存在整点。 2.边 $AB, BC, CA$ 上不存在除了端点以外的整点。 现在勇太想要知道有多少种不同的选取方案满足条件。 当然，这个问题对于萌萌哒六花来说实在是太难了，你可以帮帮她吗？ 注意 ${A, B, C}$ 与 ${B, A, C}$ 视为同一种方案。 $1\leq n,m \leq 5×10^9$

题目中限制三角形内部没有整点并且边上没有除了端点之外的点，那么根据皮克定理

三角形面积只会是$0.5$。假设三角形的包围盒大小为$a×b$，必定有$gcd(a,b)=1$，我们假设三角形的另外一个定点为$(x,y)$，三角形面积满足$S=\frac{\mid ay-bx\mid }{2}=0.5$，所以$ay-bx=\pm 1$。

根据扩展欧几里得，这个方程在满足$0\leq x\leq a \land 0\leq y \leq b$时，等式右边取-1或者1各仅有一组解。换一下包围盒中对角线的位置，又可以得到2个解，所以每个包围盒对应了四个满足条件的三角形。

所以只需要计算

反演一下，原式$=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{\lfloor n/d \rfloor}\sum_{j=1}^{\lfloor m/d \rfloor}(n-i\times d)(m-j\times d)$， 然后就可以分块求了，需要用杜教筛维护一下$\mu(n),\mu(n)\times n,\mu(n)\times n\times n$的前缀和。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD=998244353;
const int maxn=2e6+50;
bool isp[maxn];
int prime[maxn],pz,mob[maxn];
ll sum[maxn],sum1[maxn],sum2[maxn];
void init() {
	memset(isp,true,sizeof isp);
	isp[0]=isp[1]=false;
	mob[1]=sum[1]=sum1[1]=sum2[1]=1;
	for(int i = 2; i < maxn; ++i) {
		if(isp[i]) prime[pz++]=i,mob[i]=-1;
		sum[i]=sum[i-1]+mob[i];
		sum1[i]=(sum1[i-1]+mob[i]*i)%MOD;
		sum2[i]=(sum2[i-1]+(ll)mob[i]*i%MOD*i)%MOD;
		for(int j = 0; j < pz &&(ll)i*prime[j]<maxn; ++j) {
			isp[i*prime[j]]=false;
			if(i%prime[j]==0) {
				mob[i*prime[j]]=0;
				break;
			}
			else mob[i*prime[j]]=-mob[i];
		}
	}
}
map<ll,ll> mp;
ll djs(ll n) {
	if(n<maxn) return sum[n];
	if(mp.count(n)) return mp[n];
	ll ans=1;
	for(ll i = 2; i <= n; ++i) {
		ll nxt=n/(n/i);
		ans-=(nxt-i+1)*djs(n/i);
		ans%=MOD;
		i=nxt;
	}
	return mp[n]=ans;
}
ll inv2,inv6;
ll djs1(ll n) {
	if(n<maxn) return sum1[n];
	if(mp.count(n)) return mp[n];
	ll ans=1;
	for(ll i = 2; i <= n; ++i) {
		ll nxt=n/(n/i);
		ans-=(i+nxt)%MOD*((nxt-i+1)%MOD)%MOD*inv2%MOD*djs1(n/i);
		ans%=MOD;
		i=nxt;
	}
	return mp[n]=ans;
}
ll G(ll n) {
	n%=MOD;
	return n*(n+1)%MOD*(n*2+1)%MOD*inv6%MOD;
}
ll djs2(ll n) {
	if(n<maxn) return sum2[n];
	if(mp.count(n)) return mp[n];
	ll ans=1;
	for(ll i = 2; i <= n; ++i) {
		ll nxt=n/(n/i);
		ans-=(G(nxt)-G(i-1))*djs2(n/i);
		ans%=MOD;
		i=nxt;
	}
	return mp[n]=ans;
}
ll mulmod(ll x,ll y) {
	return (x%MOD)*(y%MOD)%MOD;
}
ll getInv(ll x) {
	return x==1?x:getInv(MOD%x)*(MOD-MOD/x)%MOD;
}
int main() {
	init();
	inv2=getInv(2),inv6=getInv(6);
	ll n,m;
	scanf("%lld%lld", &n,&m);
	if(n>m) swap(n,m);
	ll ans=0;
	for(ll i = 1; i <= n; ++i) {
		ll nxt=min(n/(n/i),m/(m/i));
		ans+=(djs(nxt)-djs(i-1))*mulmod(mulmod(n,n/i),mulmod(m,m/i));
		ans%=MOD;
		i=nxt;
	}
	mp.clear();
	for(ll i = 1; i <= n; ++i) {
		ll nxt=min(n/(n/i),m/(m/i));
		ans-=(djs1(nxt)-djs1(i-1))*(mulmod(mulmod(n,n/i),mulmod(m/i%MOD*inv2,m/i+1))+mulmod(mulmod(m,m/i),mulmod(n/i%MOD*inv2,n/i+1)));
		ans%=MOD;
		i=nxt;
	}
	mp.clear();
	for(ll i = 1; i <= n; ++i) {
		ll nxt=min(n/(n/i),m/(m/i));
		ans+=(djs2(nxt)-djs2(i-1))*mulmod(mulmod(n/i%MOD*inv2,n/i+1),mulmod(m/i%MOD*inv2,m/i+1));
		ans%=MOD;
		i=nxt;
	}
	cout<<(ans*4%MOD+MOD)%MOD<<endl;
	return 0;
}
Search