【Burnside引理】【Pólya计数法】【Pollard's Rho】【JZOJ 5457】 项链

Description

现在有m 种颜色的珠子。定义一个长度为n的项链为一个顺次连接n个珠子的环, 将所有旋转和翻转看作是等价的。
比如说, [1,2,3,4]通过旋转等价于[2,3,4,1],[3,4,1,2], [4,1,2,3]; [1, 2,3,4] 通过翻转等价于[1,4,3,2], [3,2,1,4], [2,1,4,3],[4,3,2,1]。
同时, 你还可以进行一种颜色转换操作. 这种操作会将所有珠子的颜色编号加1, 特别地, 对于所有颜色编号为m的珠子, 它们的颜色编号会变为1。
如果一个项链A在经过任意的旋转, 翻转, 颜色转换之后变为了项链B,则称A和B是等价的。
现在你要统计有多少个本质不同的项链, 对$998244353$取模。
对于100% 的数据，$1<=T<=20; 3<=n<=10^{18}, 2<=m<=10^{18}, 998244353 ∤ n,m$
T为数据组数

n是珠子个数，m是颜色数

群论

分析题目，有旋转、翻转、颜色转换三种不同类型的置换
旋转可以旋1~n格
翻转相当于要么从正面看，要么从反面看
颜色转换可以整体加1~m
那么，置换的总数为$2nm$

Burnside引理

我们考虑对于每种置换求出其不动点（置换前后相同的序列）个数，求和再除以2nm就是答案
由于翻转比较特殊，要么从正面看（翻转偶数次，相当于不翻转）要么从反面看（翻转奇数次，相当于翻转一次），考虑从翻转入手

不翻转-旋转i格-颜色转换加d

先不要考虑颜色转换，纯粹一点，只有旋转i格？
相信你可以感受出来它的不动点个数为$m^{(i,n)}$
如果你感受不出来也没关系，用Pólya定理解释，置换“旋转i格”可以分解成若干循环的乘积。点x，x+i必然处于同一个循环，相当于从x开始沿着环走，步长为i，经过的点都在同一个循环。由于回到原点走的步数为$\frac{lcm(i,n)}{i}=\frac{n}{(i,n)}$，所以单个循环的长度就是$\frac{n}{(i,n)}$，循环节数就是$(i,n)$
由Pólya定理，不动点个数为$m^{(i,n)}$
加上颜色转换？
当旋转i格的时候，颜色转换哪些d是可行的？显然不同的循环本质是一样的，我们只需要针对某一个循环思考
满足$d\frac{n}{(i,n)}≡0(mod\ m)$的$d$是可行的（转一圈回到自己颜色相等），自行感受
将和式写出来

$$\sum_{i=1}^nm^{(i,n)}\sum_{d=1}^m[d\frac{n}{(i,n)}≡0(mod\ m)]$$

$$=\sum_{i=1}^nm^{(i,n)}\sum_{d=1}^m[\frac{m}{(d,m)}|\frac{n}{(i,n)}]$$
$$=\sum_{x|n}m^{n/x}\varphi(x)\sum_{d=1}^m[\frac{m}{(d,m)}|x]$$
考虑对第二个 $\sum$化简
$$\sum_{d=1}^m[\frac{m}{(d,m)}|x]$$
$$=\sum_{i|x}\sum_{d=1}^m[\frac{m}{(d,m)}=i]$$
$$=\sum_{i|x,i|m}\sum_{d=1}^m[(d,m)=m/i]$$
$$=\sum_{i|(x,m)}\varphi(i)$$
$$=(x,m)$$
所以，原式化为
$$=\sum_{x|n}m^{n/x}\varphi(x)(x,m)$$
至此，不翻转部分不动点总数已经求完
但是一个细节是由于要枚举 $n$的因子，所以要对n分解质因数。而n可能有$10^{18}$，暴力分解会T。我们需要 Pollard’s Rho来分解。

翻转-旋转i格

翻转了之后不动点个数应该比较好算了
不论n是奇数还是偶数，不同翻转数为n，沿轴对称的一对点要相等

n为奇数

肯定是某个点不动，其它的翻转
显然它怎么颜色转换（除了转m次）都不可能相等
$m^{n/2+1}*n$

n为偶数

可能是沿两个数为对称轴，也可能是沿两个空隙为对称轴
$m^{n/2}*(n/2)+m^{n/2+1}*(n/2)$
如果考虑颜色转换，只有沿两个空隙为对称轴可以
并且转换两次要等于自己，那么m必须得是偶数，只有转换m/2次可行
$m^{n/2}*(n/2)$

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<ctime>
#include<cmath>
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define fd(i,b,a) for(ll i=b;i>=a;i--)
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define mset(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
const ll mo=998244353;
ll qmul(ll x,ll y,ll m)
{
    ll t=0;
    for(;y;y>>=1,x=(x<<1)%m)
        if(y&1) t=(t+x)%m;
    return t;
}
ll qmi(ll x,ll n)
{
    ll t=1;
    for(x%=mo;n;n>>=1,x=x*x%mo)
        if(n&1) t=t*x%mo;
    return t;
}
ll qpow(ll x,ll n,ll m)
{
    ll t=1;
    for(x%=m;n;n>>=1,x=qmul(x,x,m))
        if(n&1) t=qmul(t,x,m);
    return t;
}
ll gcd(ll x,ll y)
{
    if(x%y==0) return y;
    else return gcd(y,x%y);
}
ll n,m,ans,num,a[20],b[20];
bool check(ll n,ll p,ll u,ll t)
{
    ll x=qpow(n,t,p);
    fo(i,1,u)
    {
        ll nxt=qmul(x,x,p);
        if(nxt==1 && x!=1 && x!=p-1) return 0;
        x=nxt;
    }
    return x==1;
}
bool miller_rabin(ll p)
{
    if(p==2) return 1;
    if(p<2 || !(p&1)) return 0;
    ll u=0,t=p-1;
    for(;t%2==0;t/=2) u++;
    fo(i,1,10)
    {
        ll n=rand()%(p-2)+2;
        if(!check(n,p,u,t)) return 0;
    }
    return 1;
}
ll pollard_rho(ll n,ll c)
{
    int i=1,k=2;
    ll x=rand()%n;ll y=x;
    for(;;)
    {
        i++;
        x=(qmul(x,x,n)+c)%n;
        ll d=gcd(abs(x-y),n);
        if(d!=1 && d!=n) return d;
        if(y==x) return n;
        if(i==k) y=x,k<<=1;
    }
}
void find(ll n)
{
    if(n==1) return;
    if(miller_rabin(n))
    {
        a[++num]=n;
        return;
    }
    ll d=n;
    while(d>=n) d=pollard_rho(n,rand()%(n-1)+1);
    find(d);
    while(n%d==0) n/=d;
    find(n);
}
void getpr(ll n)
{
    num=0;
    find(n);
    sort(a+1,a+num+1);
    num=unique(a+1,a+num+1)-a-1;
    fo(i,1,num)
        for(b[i]=0;n%a[i]==0;n/=a[i]) b[i]++;
}
void dfs(int x,ll d,ll phi)
{
    if(x>num)
    {
        ans=(ans+gcd(d,m)%mo*phi%mo*qmi(m,n/d)%mo)%mo;
        return;
    }
    ll t=1,xphi=phi*((a[x]-1)%mo)%mo*qmi(a[x],mo-2)%mo;
    fo(i,0,b[x])
    {
        if(!i) dfs(x+1,d*t,phi);
        else dfs(x+1,d*t,t%mo*xphi%mo);
        t=t*a[x];
    }
}
ll calc()
{
    ans=0;
    getpr(n);
    dfs(1,1,1);
    if(n&1) ans=(ans+n%mo*qmi(m,n/2+1)%mo)%mo;
    else
    {
        ans=(ans+(n/2)%mo*(qmi(m,n/2+1)+qmi(m,n/2)%mo)%mo)%mo;
        if(n%2==0 && m%2==0) ans=(ans+(n/2)%mo*qmi(m,n/2)%mo)%mo;
    }
    ans=ans*qmi(2*n%mo*(m%mo)%mo,mo-2)%mo;
    return ans;
}
int main()
{
    srand(time(0));rand();
    freopen("necklace.in","r",stdin);
    freopen("necklace.out","w",stdout);
    int T;
    for(scanf("%d",&T);T;T--)
    {
        scanf("%lld %lld",&n,&m);
        printf("%lld\n",calc());
    }
    return 0;
}