Newcoder 147 J.Juggernaut（polya）

Description

用$0,1$填一个$n\times m$的矩阵，使得每行每列元素异或和为$0$，且对于两个矩阵$A,B$，若存在$x\in [0,n),y\in [0,m)$使得$\forall 0\le i<n,0\le j<m$都有$A_{i,j}=B_{(i+x)\%n,(j+y)\%m}$，则$A,B$视为一种方案，问方案数

Input

输入两个整数$n,m(1\le n,m\le 10^9)$

Output

输出方案数，结果模$998244353$

Sample Input

4 4

Sample Output

48

Solution

$Polya$，有$nm$种变换群，先不考虑行列异或和的限制，对于平移$(x,y),1\le n\le n,1\le y\le m$，求出在该变换下不动的方案数，第一维需要$\frac{n}{gcd(x,n)}$次平移才能回到起点，第二维需要$\frac{m}{gcd(y,m)}$次平移才能回到起点，那么一个循环的长度即为$lcm(\frac{n}{gcd(x,n)},\frac{m}{gcd(y,m)})$，故此时染色方案数为$2^{\frac{nm}{lcm(\frac{n}{gcd(x,n)},\frac{m}{gcd(y,m)})}}$，枚举$\frac{n}{gcd(x,n)}$和$\frac{m}{gcd(y,m)}$，则有
$$ans=\frac{1}{nm}\sum _{a|n}\sum _{b|m}\phi (a)\phi (b)2^{\frac{nm}{lcm(a,b)}}$$
现在考虑行列异或和的限制，要修改的即为指数部分$\frac{nm}{lcm(a,b)}$，对于位置$(i,j)$和一个平移$(\frac{n}{a},\frac{m}{b})$，第一维需要$a$次平移回到起点，第二维需要$b$次平移回到起点，那么值需要与$(i,j)$相同且第一维为$i$的有$\frac{lcm(a,b)}{a}$个，第二维为$j$的有$\frac{lcm(a,b)}{b}$个，记$l=lcm(a,b)$

1.若$\frac{l}{a}$为奇数，$\frac{l}{b}$为奇数，由于每$\frac{n}{a}$行的状态相同，每$\frac{m}{b}$列的状态相同，故只要前$\frac{n}{a}$行和前$\frac{m}{b}$列的异或和为$0$即可，考虑第一列前$\frac{n}{a}$个元素和第一行前$\frac{m}{b}$个元素，显然这$\frac{n}{a}+\frac{m}{b}-1$个元素取值相互独立且依赖于右下角$(\frac{n}{a}-1)\times (\frac{m}{b}-1)$这个子矩阵的取值，故此时修改为$\frac{nm}{l}-(\frac{n}{a}+\frac{m}{b}-1)$

2.若$\frac{l}{a}$为奇数，$\frac{l}{b}$为偶数，则每个第二维为$j$的值都会出现偶数次，也即每列异或和必然为$0$，但是行和需要做同样的考虑，即第一列前$\frac{n}{a}$个元素的取值依赖于右边$\frac{n}{a}\times (\frac{m}{b}-1)$这个子矩阵的取值，此时修改为$\frac{nm}{l}-\frac{n}{a}$

3.若$\frac{l}{a}$为偶数，$\frac{l}{b}$为奇数，则每个第一维为$i$的值都会出现偶数次，也即每行异或和必然为$0$，但是列和需要做同样的考虑，即第一行前$\frac{m}{b}$个元素的取值依赖于下边$(\frac{n}{a}-1)\times \frac{m}{b}$这个子矩阵的取值，此时修改为$\frac{nm}{l}-\frac{m}{b}$

由于结果要模$p$，当$n=m=p$时很麻烦，但是用例已经给出；当$n,m$中有一个为$p$时，不妨设$n=p$，那么以$p^2$为模数先得到答案，之后除以$p$再乘上$m$在模$p$意义下逆元即可；当$n,m$均不等于$p$时，直接在模$p$意义下求解即可

Code

#include<cstdio>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef pair<int,int>P;
ll mod=998244353;
ll add(ll x,ll y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
ll mul(ll x,ll y)
{
	return (x*y-(ll)(x/(ld)mod*y+1e-3)*mod+mod)%mod;
}
ll Pow(ll x,ll y)
{
	ll ans=1;
	while(y)
	{
		if(y&1)ans=mul(ans,x);
		x=mul(x,x);
		y>>=1;
	}
	return ans;
}
vector<P>f[2];
vector<int>fact;
int phi(int n)
{
	int ans=n;
	for(int i=0;i<fact.size();i++)
		if(n%fact[i]==0)ans=ans/fact[i]*(fact[i]-1);
	return ans;
}
void deal(int n,int id)
{
	fact.clear();
	int nn=n;
	for(int i=2;i*i<=nn;i++)
		if(nn%i==0)
		{
			fact.push_back(i);
			while(nn%i==0)nn/=i;
		}
	if(nn>1)fact.push_back(nn);
	for(int i=1;i*i<=n;i++)
		if(n%i==0)
		{
			f[id].push_back(P(i,phi(i)));
			if(i*i!=n)f[id].push_back(P(n/i,phi(n/i)));
		}
}
int gcd(int a,int b)
{
	return b?gcd(b,a%b):a;
}
int main()
{
	int n,m;
	scanf("%d%d",&n,&m);
	if(n==mod&&m==mod)
	{
		printf("295980207\n");
		return 0;
	}
	if(m==mod)swap(n,m);
	if(n==mod)mod=mod*n;
	ll ans=0;
	deal(n,0),deal(m,1);
	for(int i=0;i<f[0].size();i++)
		for(int j=0;j<f[1].size();j++)
		{
			ll res=mul(f[0][i].second,f[1][j].second);
			int a=f[0][i].first,b=f[1][j].first;
			ll l=(ll)a*b/gcd(a,b);
			ll t=(ll)n*m/l;
			if(((l/a)&1)&&((l/b)&1))t-=n/a+m/b-1;	
			else if((l/a)&1)t-=n/a;
			else if((l/b)&1)t-=m/b;
			ans=add(ans,mul(res,Pow(2,t)));
		}
	if(n==998244353)ans/=n,mod/=n;
	else ans=mul(ans,Pow(n,mod-2));
	ans=mul(ans,Pow(m,mod-2));
	printf("%lld\n",ans);
	return 0;
}