题解 | #The Pool#

出题人的题解实在是无法令人恭维，特此写一份自己的题解

【大意】

$TT$ 次询问，每次询问给定 $n,m(1\le n,m\le 10^{18})n,\; m(1\backslash leq\; n,m\backslash leq\; 10^\{18\})$ ，问长宽分别为 $n,mn,\; m$ 的矩形顶点摆放在整点后；所有不同摆放方案中，每个方案完全包含的 $1\times 11\backslash times\; 1$ 格子数量的和是多少？

我们认为两个方案不同，当且仅当一个方案的矩形无法仅通过平移变换，得到另一个方案。

【分析】

我们将矩形 $ABCDABCD$ 的一个顶点放置在原点上，设 $AB=n,AD=mAB=n,\; AD=m$ ，则 $BB$ 点的点坐标需要满足 $x_B^2+y_B^2=n^{2}x\_B^2+y\_B^2=n^2$

如果我们能找到所有满足条件的点坐标，那它们显然都是答案的候选（可能有的方案能通过别的方案平移得到、有的方案无法使得 $DD$ 在整点）

因此，当前的问题即是：给定 $ll$ ，如何求解所有满足 $x_B^2+y_B^2=lx\_B^2+y\_B^2=l$ 的点

转化一下问题，由于二维平面上的整点能唯一地对应复平面上，满足 $\mathrm{\Re }z,\mathrm{\Im }z\in Z\backslash Re\; z,\; \backslash Im\; z\backslash in\; Z$ 的所有复数 $zz$ 。因此，我们可以直接以高斯整数 $Z[i]Z[i]$ 的方式描述待求问题：

求解有多少个高斯整数 $z=x_B+i\cdot y_{B}z=x\_B+i\backslash cdot\; y\_B$ ，使得 $z\cdot \bar{z}=lz\backslash cdot\; \backslash bar\; z=l$

我们考虑将 $ll$ 进行质因数分解，显然可以将 $ll$ 分解为 ${\displaystyle l=2^{c_2}\cdot \prod _j{p}_j^{k_j}\cdot \prod _j{q}_j^{t_j}}\backslash displaystyle\; l=2^\{c\_2\}\backslash cdot\; \backslash prod\_j\; p\_j^\{k\_j\}\; \backslash cdot\; \backslash prod\_j\; q\_j^\{t\_j\}$ 的形式，其中 $c_2\ge 0,k_j,t_j>0,q_{j}c\_2\backslash geq\; 0,\; k\_j,\; t\_j>0,\; q\_j$ 为 $4n+34n+3$ 型质数，$p_{j}p\_j$ 为 $4n+14n+1$ 型质数

根据高斯质数的规律，$2=(1+i)(1-i),q_{j}2=(1+i)(1-i),\; q\_j$ 无法继续分解，而 $p_{j}p\_j$ 能进一步分解为两个共轭的高斯质数 $(a+bi),(a-bi)(a+bi),(a-bi)$ 的乘积

为了对 $p_{j}p\_j$ 进行高斯质数意义上的质因数分解，我们可以随机一个 $t\mathrm{\ne }0(\mathrm{m}\mathrm{o}\mathrm{d}{p}_j)t\backslash neq\; 0\backslash pmod\; \{p\_j\}$ ，再令 ${\displaystyle k=t^{\frac{p_j-1}{4}}\mathrm{m}\mathrm{o}\mathrm{d}{p}_j}\backslash displaystyle\; k=t^\{p\_j-1\backslash over\; 4\}\backslash bmod\; p\_j$ 。根据二次剩余的理论，$k^2\equiv \pm 1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_j)k^2\backslash equiv\; \backslash pm\; 1\backslash pmod\; \{p\_j\}$ ，取值为 $-1-1$ 的概率恰为 $\frac{1}{2}\{1\backslash over\; 2\}$ 。

所以，我们期望两次随机，就能选出一个 $kk$ ，使得 $k^2\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_j)k^2\backslash equiv\; -1\backslash pmod\; \{p\_j\}$ ，即 $p_j\mid (k^2+1)p\_j\backslash mid\; (k^2+1)$ 。由于在 $Z[i]Z[i]$ 上，$k^2+1=(k+i)(k-i)k^2+1=(k+i)(k-i)$ ，故有 $p_j\mid (k+i)(k-i)p\_j\backslash mid\; (k+i)(k-i)$

而 ，能推出不存在 $z\in Z[i]z\backslash in\; Z[i]$ 使得 $p_j\cdot z=k+ip\_j\backslash cdot\; z=k+i$ ，从而 $p_{j}p\_j$ 并不整除 $k+ik+i$ ，同理 $p_{j}p\_j$ 并不整除 $k-ik-i$ ；所以 $(a+bi),(a-bi)(a+bi),\; (a-bi)$ 其中一者整除 $k+ik+i$ ；又因为两者均是高斯质数，不可再分解；故直接求解 $\gcd (p_j,k+i)\backslash gcd(p\_j,\; k+i)$ 即可得到两者之一，再用 $p_{j}p\_j$ 除即可再得到另一个解

根据上述理论，我们可以将 $ll$ 进一步分解为 ${\displaystyle l=(1+i{)}^{c_2}\cdot (1-i{)}^{c_2}\cdot \prod _j(a+bi{)}^{k_j}(a-bi{)}^{k_j}\cdot \prod _j{q}_j^{t_j}}\backslash displaystyle\; l=(1+i)^\{c\_2\}\backslash cdot\; (1-i)^\{c\_2\}\backslash cdot\; \backslash prod\_j\; (a+bi)^\{k\_j\}(a-bi)^\{k\_j\}\backslash cdot\; \backslash prod\_j\; q\_j^\{t\_j\}$ ，我们需要求解有多少个高斯整数 $zz$ 使得 $z\cdot \bar{z}=lz\backslash cdot\; \backslash bar\; z=l$

简单起见，我们不妨先考虑不含因子 $22$ 的情况：

对于 $(a+bi)(a+bi)$ 与 $(a-bi)(a-bi)$ 型的高斯质数，我们显然每次需要将其中一个分配给 $zz$ ，另一个分配给 $\bar{z}\backslash bar\; z$ ，才能使得两者保持共轭的性质；因此我们可以枚举 $(a+bi)(a+bi)$ 中，有 $t(0\le t\le k_j)t(0\backslash leq\; t\backslash leq\; k\_j)$ 个分配给了 $zz$ ，剩下的分配给了 $\bar{z}\backslash bar\; z$ ；因此 $zz$ 中的方案为 $(a+bi{)}^t\cdot (a-bi{)}^{k_j-t}(a+bi)^t\backslash cdot\; (a-bi)^\{k\_j-t\}$ ，方案数为 $k_j+1k\_j+1$

对于 $q_j^{t_j}q\_j^\{t\_j\}$ 型高斯质数，当且仅当 $t_{j}t\_j$ 为偶数时可以平分给 $zz$ 与 $\bar{z}\backslash bar\; z$ ，产生唯一的方案；否则分配是不平均的，方案数为 $00$

由于式子 $z\cdot \bar{z}=m=n^{2}z\backslash cdot\; \backslash bar\; z=m=n^2$ ，故方案数必然不为 $00$ ；然而对于每个求解出的合法解 $zz$ ，我们乘上单位数 $1,i,-1,-i1,\; i,\; -1,\; -i$ ，并向 $\bar{z}\backslash bar\; z$ 乘上相应单位数的共轭，都能得到一组互不相同的解（相伴解）；因此我们即可得出所有的合法方案，方案数为 ${\displaystyle 4\prod _j(k_j+1)}\backslash displaystyle\; 4\backslash prod\_j\; (k\_j+1)$

对于含 $22$ 因子时，同样我们可以考虑分配多少个 $(1+i)(1+i)$ 给 $zz$ ，剩下的给 $\bar{z}\backslash bar\; z$ ；然而两者交换一个 $(1+i),(1-i)(1+i),(1-i)$ 后，我们发现答案实际上只相差一个 $ii$ ，这个的贡献其实在后续乘单位数中已经进行了消除；为此，我们不妨将所有的 $(1+i)(1+i)$ 分配给 $zz$ ，因为其不影响方案数

根据题解，方案数为 ${\displaystyle 4\prod _j(k_j+1)}\backslash displaystyle\; 4\backslash prod\_j(k\_j+1)$ ，最大值为 $4\times 2952454\backslash times\; 295245$ ，在 $n=495229111954868525n=495229111954868525$ 时取得

现在，我们已经求解所有满足 $x_B^2+y_B^2=lx\_B^2+y\_B^2=l$ 的点，需要在这些候选点中，选择无法通过别的方案平移得到、$DD$ 在整点的方案

对于 $n=mn=m$ 的情况，显然，角 $\mathrm{\angle }BAx\in [0,\frac{\pi }{2})\backslash angle\; BAx\backslash in\; [0,\; \{\backslash pi\backslash over\; 2\})$ 都对应互不相同的矩形。我们从候选点中，选择那些 $\mathrm{\Im }z\ge 0\wedge \mathrm{\Re }z>0\backslash Im\; z\backslash geq\; 0\backslash wedge\; \backslash Re\; z>0$ 的解即可。

而对于 $n\mathrm{\ne }mn\backslash neq\; m$ 的情况，角 $\mathrm{\angle }BAx\in [0,\pi )\backslash angle\; BAx\backslash in[0,\; \backslash pi)$ 都对应互不相同的矩形；同理可以从候选点中，选择 $\mathrm{\Im }z>0\vee (\mathrm{\Im }z=0\wedge \mathrm{\Re }z>0)\backslash Im\; z>0\backslash vee\; (\backslash Im\; z=0\backslash wedge\; \backslash Re\; z>0)$ 的解。

那如何确认候选点 $(a+bi)(a+bi)$ 是否满足 $DD$ 也位于整点呢？

根据几何关系，我们直接将候选点乘上单位数 $ii$ 即可旋转 $\frac{pi}{2}\{pi\backslash over\; 2\}$ ，再乘上 $\frac{m}{n}\{m\backslash over\; n\}$ 即可放缩至 $DD$ 点；因此只需要验证 $(a+bi)\cdot i\cdot \frac{m}{n}(a+bi)\backslash cdot\; i\backslash cdot\; \{m\backslash over\; n\}$ 是否为整点即可

终于来到了最终答案求解的部分，已知 $BB$ 位于 $(a^{\mathrm{\prime }}+b^{\mathrm{\prime }}i)(a\text{'}+b\text{'}i)$ ，$DD$ 位于 $(c^{\mathrm{\prime }}+d^{\mathrm{\prime }}i)(c\text{'}+d\text{'}i)$ ，如何求解答案包含的 $1\times 11\backslash times\; 1$ 格子数量的和？

为了方便，我们不妨假设 $a=\mathrm{\mid }{a}^{\mathrm{\prime }}\mathrm{\mid },b=\mathrm{\mid }{b}^{\mathrm{\prime }}\mathrm{\mid },c=\mathrm{\mid }{c}^{\mathrm{\prime }}\mathrm{\mid },d=\mathrm{\mid }{d}^{\mathrm{\prime }}\mathrm{\mid }a=|a\text{'}|,\; b=|b\text{'}|,\; c=|c\text{'}|,\; d=|d\text{'}|$ ，它们有几何关系如下：

我们按下述方法分别统计每个三角形对答案的贡献

红色三角形的贡献是 $y=\frac{b}{a}x,x\in [0,a]y=\{b\backslash over\; a\}x,x\backslash in[0,\; a]$ 内的格子数，蓝色的是 $y=\frac{d}{c}x,x\in [0,c]y=\{d\backslash over\; c\}x,x\backslash in[0,\; c]$ 内的格子数，很显然，有一部分会被重复计算或遗漏，需要再扣除一个 $(c-a)(d-b)(c-a)(d-b)$

而关于三角形内部的贡献，我们考虑直线的斜率是不少于 $00$ 的，因此只要格子的左上端点位于直线 $y=\frac{b}{a}xy=\{b\backslash over\; a\}x$ 下方（或线上）则能被统计到

故答案变为：除最下行、最右列外，直线下方的整点数

显然是可以由类欧或万欧解决，这里介绍一下 pick 定理的解法：

对于这类所有顶点都在整点上的多边形，称为格点多边形，其面积满足公式：$S=\frac{1}{2}B+I-1S=\{1\backslash over\; 2\}B+I-1$ ，其中 $BB$ 为边缘点数，$II$ 为多边形内部点数

该直线面积显然为 $\frac{ab}{2}\{ab\backslash over\; 2\}$ ，边缘点数可以发现是 $(a+1)+(b+1)+(\gcd (a,b)+1)-3=a+b+\gcd (a,b)(a+1)+(b+1)+(\backslash gcd(a,b)+1)-3=a+b+\backslash gcd(a,b)$ ，故内部点数为 $I=\frac{ab-a-b-\gcd (a,b)}{2}+1I=\{ab-a-b-\backslash gcd(a,b)\backslash over\; 2\}+1$

而可以发现，格子的左上端点仅由内部点、斜线上的非端点构成；因此格子数为 $I+\gcd (a,b)-1=\frac{ab-a-b+\gcd (a,b)}{2}I+\backslash gcd(a,b)-1=\{ab-a-b+\backslash gcd(a,b)\backslash over\; 2\}$

综上，我们可以先对 $nn$ 进行 Pollard_Rho 算法进行质因数分解；暴力求解所有高斯整数后，选出合法点，并利用上述公式求解答案。

【代码】

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define mp make_pair
#define de(x) cout << #x <<" = "<<x<<endl
#define dd(x) cout << #x <<" = "<<x<<" "
mt19937 rnd(time(0));
const int P=998244353;
inline __int128 absl(__int128 x) { return x<0?-x:x; }
inline ll kpow(ll a, ll x, ll p) { ll ans=1; for(;x;x>>=1, a=(__int128)a*a%p) if(x&1) ans=(__int128)ans*a%p; return ans; }
struct Zi {
	__int128 re, im;
	inline Zi (ll re_=0, ll im_=0):re(re_), im(im_) {}
	inline Zi& operator += (const Zi &rhs) { re+=rhs.re; im+=rhs.im; return *this; }
	inline Zi& operator -= (const Zi &rhs) { re-=rhs.re; im-=rhs.im; return *this; }
	inline Zi& operator *= (const Zi &rhs) { tie(re, im)=mp(re*rhs.re-im*rhs.im, re*rhs.im+im*rhs.re); return *this; }
	inline Zi& operator /= (const Zi &rhs) {
		__int128 r=rhs.re, i=rhs.im, l=r*r+i*i;
		tie(r, i)=mp(2*(re*r+im*i)+l, 2*(im*r-re*i)+l); l*=2;
		if(r<0) re=-((-r+l-1)/l);
		else re=r/l;
		if(i<0) im=-((-i+l-1)/l);
		else im=i/l;
		return *this;
	}
	inline Zi& operator %= (const Zi &rhs) { Zi tmp=*this; tmp/=rhs; tmp*=rhs; return *this-=tmp; }
	
	inline Zi operator + (const Zi &rhs) const { Zi lhs=*this; return lhs+=rhs; }
	inline Zi operator - (const Zi &rhs) const { Zi lhs=*this; return lhs-=rhs; }
	inline Zi operator * (const Zi &rhs) const { Zi lhs=*this; return lhs*=rhs; }
	inline Zi operator / (const Zi &rhs) const { Zi lhs=*this; return lhs/=rhs; }
	inline Zi operator % (const Zi &rhs) const { Zi lhs=*this; return lhs%=rhs; } 
};
inline Zi kpow(Zi a, ll x) { Zi ans=1; for(;x;x>>=1, a=a*a) if(x&1) ans=ans*a; return ans; }
inline Zi norm(Zi x) {
	__int128 r=x.re, i=x.im;
	if(r==0) return Zi(absl(i), 0);
	switch((r<0)<<1|(i<0)) {
		case 0: return Zi(r, i); break;
		case 1: return Zi(-i, r); break;
		case 2: return Zi(i, -r); break;
		case 3: return Zi(-r, -i); break;
	}
	return Zi();
}
inline Zi gcd(Zi a, Zi b) {
	while(b.re||b.im)
		swap(a=norm(a%b), b);
	return a;
}
inline Zi GaussPrime(ll p) {
	if(p==2)
		return Zi(1, 1);
	if(p%4!=1)
		return Zi(p, 0);
	ll k=0;
	while(((__int128)k*k+1)%p!=0) k=kpow(rnd(), p>>2, p);
	return gcd(Zi(p, 0), Zi(k, 1));
}

inline bool Fermat(ll n,int p){
	if( kpow(p,n-1,n)!=1 ) return 0;
	ll k=n-1;
	k>>=__builtin_ctzll(k);
	ll t=kpow(p, k, n);
	if(t==1) return 1;
	for(;k<n-1;k<<=1, t=(__int128)t*t%n)
		if(t==n-1) return 1;
	return 0;
}
inline bool Miller_Rabin(ll n){
	static int p[]={2, 3, 5, 7, 11, 13, 17, 19, 23, 61, 233}, siz=sizeof(p)/sizeof(p[0]);
	if(n==1) return 0;
	for(int i=0;i<siz;++i)
		if(n%p[i]==0) return n==p[i];
	else if( !Fermat(n,p[i]) ) return 0;
	return 1;
}
inline ll Rho(ll n,ll c){
	static ll lim=128, tmp[128];
	ll buf=1, tot=0;
	for(ll f1=(c==n-1?0:n+1), f2=((__int128)f1*f1+c)%n; f1!=f2;
		f1=((__int128)f1*f1+c)%n, f2=((__int128)f2*f2+c)%n, f2=((__int128)f2*f2+c)%n){
		tmp[tot++]=f1-f2;
		if(tmp[tot-1]<0) tmp[tot-1]+=n;
		buf=(__int128)buf*tmp[tot-1]%n;
		if(tot==lim){
			if( __gcd(buf,n)>1 ) break;
			tot=0;
		}
	}
	for(int i=0;i<tot;++i){
		buf=__gcd(tmp[i], n);
		if(buf>1) return buf;
	}
	return n;
}
void Pollard_Rho(ll n, vector<pair<ll, int> > &pf, ll cnt){
	static int p[]={2, 3, 5, 7, 11, 13, 17, 19, 23, 61, 233}, siz=sizeof(p)/sizeof(p[0]);
	if(n==1) 
		return ;
	if( Miller_Rabin(n) ) {
		pf.emplace_back(n, cnt);
		return ;
	}
	ll d=n;
	int cntd=0;
	for(int i=0;i<siz;++i) if(n%p[i]==0){
		d=p[i];
		cntd=0;
		while(n%d==0) n/=d, ++cntd;
		pf.emplace_back(d, cnt*cntd);
	}
	while(d==n) d=Rho(n, rand()%(n-1)+1);
	cntd=0;
	while(n%d==0) n/=d, ++cntd;
	Pollard_Rho(d, pf, cnt*cntd);
	Pollard_Rho(n, pf, cnt);
}
inline ll tri2(ll a, ll b) {
	if(a<=1||b==0) return 0;
	return ((a%P)*(b%P)-a-b+__gcd(a, b))%P;
}

vector<pair<ll, int> > pf;
vector<Zi> tmp;
inline void solvePrime(ll n, vector<Zi> &v) {
	pf.clear();
	Pollard_Rho(n, pf, 1);
	
	v.clear();
	v.push_back(Zi(1, 0));
	Zi q, qk, r;
	for(auto &[p, k] : pf) {
		k*=2;
		if(p==2) {
			q=kpow(Zi(1, 1), k);
			for(auto &e : v)
				e*=q;
			continue;
		}
		if(p%4!=1) {
			q=kpow(Zi(p, 0), k/2);
			for(auto &e : v)
				e*=q;
			continue;
		}
		
		qk=kpow(q=GaussPrime(p), k);
		r=Zi(p, 0)/q;
		tmp=v;
		v.clear();
		v.reserve((k+1)*tmp.size());
		for(int i=0; i<=k; ++i) {
			for(const auto &e : tmp)
				v.push_back(e*qk);
			qk=qk/q*r;
		}
	}
	for(auto &e : v)
		e=norm(e);
	sort(v.begin(), v.end(), [](const Zi &a, const Zi &b) {
		if(a.re!=b.re) return a.re<b.re;
		else return a.im<b.im;
		});
	int cnt=1;
	for(int i=1, I=v.size(); i<I; ++i)
		if(v[i].re!=v[i-1].re||v[i].im!=v[i-1].im)
			v[cnt++]=v[i];
	v.erase(v.begin()+cnt, v.end());
}
vector<Zi> a;
inline ll calc(ll a, ll b, ll c, ll d) {
	return (tri2(a, b)+tri2(c, d)-((a-c)%P)*((b-d)%P))%P;
}
inline ll ans() {
	ll n, m;
	cin>>n>>m;
	solvePrime(n, a);
	if(n!=m) {
		int siz=a.size();
		a.resize(siz*2);
		for(int i=0, I=siz; i<I; ++i)
			a[i+siz]=a[i]*Zi(0, 1);
	}
	
	ll res=0;
	Zi c;
	for(auto z : a) {
		c=z*Zi(0, 1);
		if((__int128)c.re*m%n!=0||(__int128)c.im*m%n!=0)
			continue;
		c=Zi((__int128)c.re*m/n, (__int128)c.im*m/n);
		res=(res+calc(absl(z.re), absl(z.im), absl(c.re), absl(c.im)))%P;
	}
	if(res<0) res+=P;
	return res;
}
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	int T; cin>>T;
	while(T--)
		cout<<ans()<<"\n";
	cout.flush();
	return 0;
}