[JZOJ5457]【NOIP2017提高A组冲刺11.6】项链

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。

Solution

巨坑慎入。。。

有个叫burnside引理的东西。。

此处不再赘述，附上链接。
http://blog.csdn.net/hzj1054689699/article/details/78480961
总置换数肯定是2*N*M,就是翻转与否，颜色转换多少次，旋转多少次

方便讨论，将翻转次数设为0到1，旋转次数设为1到N（就是不转），颜色转换次数设为1到M（就是不转换）

那么推式子，把旋转+颜色转换（可以是m次），翻转+旋转，旋转+翻转+颜色转换（不为m次）看作置换分开讨论

求的是每种置换不动点的个数，此处我们将一个序列看成一个点，置换不动点也就是经过置换序列不变

翻转的很好弄
分偶数和奇数讨论

下面的公式请自行感受。。这题搞得我倦生了。。

n种旋转，n/2对数要相同

奇数就是$m^{n/2+1}*n$
因为奇数翻转有一个数会不动，那么它怎么颜色转换（除非转换M次）都不会等于自己

偶数比较麻烦
偶数翻转有可能全部动，也有可能有两个数作为对称轴不动
就是$m^{n/2}*(n/2)+m^{n/2+1}*(n/2)$

偶数还要考虑颜色转换
只有全部动可以
并且转换两次要等于自己，那么m必须得是偶数，只有转换m/2次可行
就是$m^{n/2}*(n/2)$

重点在于旋转
枚举旋转多少次
$x,x+i,x+2i,x+3i...x+ki$
直到$x+ki\equiv x\pmod n$

方便起见$(i,n)=\gcd(i,n)$

可以很容易算出$k=n/(i,n)$
那么总共就有$(i,n)$组这样的环

每种环染色，那么就是$m^{(i,n)}$

考虑颜色转换
设转换了$j$次，那么转换$k*j$要回到自己
即$j*n/(i,n)\equiv 0\pmod m$

那么总的式子就是$\sum\limits_{i=1}^{n}m^{(i,n)}\sum\limits_{j=1}^{m}[j*n/(i,n)\equiv 0\pmod m]$

考虑化简
设$d=n/(i,n)$
那么$\sum\limits_{i=1}^{n}[n/(i,n)=d]=φ(d)$
这个可以自己推导一下

那么原式化为$\sum\limits_{d|n}m^{n/d}φ(d)\sum\limits_{j=1}^{m}[jd\equiv 0]\pmod m$

解同余方程
$jd\equiv 0\pmod m$
$j*d/(d,m)\equiv 0 \pmod {{m\over(d,m)}}$
那么$j\equiv 0\pmod {{m\over(d,m)}}$
$j=k*{m\over(d,m)},k\in N^*$

又有$1\leq j\leq m$
所以$j$总共有$(d,m)$种取值

最终的式子就是
$\sum\limits_{d|n}m^{n/d}φ(d)*(d,m)$

可以用pollard_rho算法分解质因数，然后枚举约数即可

关于pollard_rho算法Miller_Rabin测试，参看
http://blog.csdn.net/hzj1054689699/article/details/78483767

Code

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <ctime>
#include <map>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define mo 998244353
#define LL long long
using namespace std;
int d[5]={2,3,5,7,10007};
LL pr[70],p1[70][2],p2[70][2],l1,l2,ans,n,m;
map<LL,LL> h;
LL gcd(LL x,LL y)
{
    return(!y)?x:gcd(y,x%y);
}
LL ti(LL x,LL y,LL m)
{
    if(y==0||x==0) return 0;
    LL s=ti(x,y/2,m);
    return (y%2)?((s+s)%m+x)%m:(s+s)%m;
}
LL ksm(LL k,LL n,LL m)
{
    LL s=1;
    for(;n;n>>=1,k=ti(k,k,m)) if(n&1) s=ti(s,k,m);
    return s;
}
LL rd(LL k)
{
    return ((rand()*rand()%k*rand()%k*rand()%k)+rand())%k;
}
bool pd(LL k)
{
    fo(i,0,4)
    {
        if(k==d[i]) return 1;
        if(k%d[i]==0) return 0;
        LL t,n=k-1;
        while((n&1)==0) n>>=1;
        t=ksm(d[i],n,k);
        while(n<k-1&&t!=1&&t!=k-1) t=ti(t,t,k),n<<=1;
        if(!((n&1)||t==k-1)) return 0;
    }   
    return 1;
}
LL get(LL k)
{
    LL c=rd(1e9)+1,r,x,y,i=1,n=2;
    x=y=rd(k-1)+1;
    while(1)
    {
        i++;
        x=(ti(x,x,k)+c)%k;
        if(x==y) return k;
        r=gcd(abs(x-y),k);
        if(1<r&&r<k) 
        {
            return r;
        }
        if(i==n) y=x,n<<=1;
    }
}
void fd(LL k)
{
    if(k==1) return;
    if(pd(k))
    {
        pr[++pr[0]]=k;
        return;
    }
    LL p=k;
    while(p==k&&p!=1) p=get(k);
    fd(k/p),fd(p);
}
LL km(LL k,LL n)
{
    LL s=1;
    k%=mo;
    for(;n;k=k*k%mo,n>>=1) if(n&1) s=s*k%mo;
    return s;
}
void dfs(int k,LL s,LL g,LL phi)
{
    if(k>l1)
    {
        (ans+=ti(phi,g,mo)*km(m,n/s)%mo)%=mo;
        return;
    }
    dfs(k+1,s,g,phi);
    LL l=(h.count(p1[k][0])==0)?0:h[p1[k][0]];
    fo(i,1,p1[k][1])
    {
        s*=p1[k][0];
        if(i<=l) g*=p1[k][0];
        if(i==1)phi*=(p1[k][0]-1);
        else phi*=p1[k][0];
        dfs(k+1,s,g,phi);
    }
}
int main()
{
    int t;
    srand(time(0));
    cin>>t;
    while(t--)
    {
        n,m;
        scanf("%lld%lld",&n,&m);
        LL cnt=(LL)2*(n%mo)*(m%mo);
        ans=0;
        pr[0]=0;
        l1=l2=0;
        h.clear();
        fd(n);
        sort(pr+1,pr+pr[0]+1);
        l1=1;
        p1[1][0]=pr[1],p1[1][1]=1;
        fo(i,2,pr[0])
        {
            if(pr[i]!=pr[i-1]) p1[++l1][0]=pr[i],p1[l1][1]=1;
            else p1[l1][1]++;
        }
        pr[0]=0;
        fd(m);
        sort(pr+1,pr+pr[0]+1);
        l2=1;
        p2[1][0]=pr[1],p2[1][1]=1;
        fo(i,2,pr[0])
        {
            if(pr[i]!=pr[i-1]) p2[++l2][0]=pr[i],p2[l2][1]=1;
            else p2[l2][1]++;
        }
        fo(i,1,l2) h[p2[i][0]]=p2[i][1];
        dfs(1,1,1,1);
        if(n%2)
        {
            ans=(ans+km(m,n/2+1)*(LL)(n%mo))%mo;
        }
        else
        {
            ans=(ans+km(m,n/2)*(LL)((n/2)%mo)%mo+km(m,n/2+1)*(LL)((n/2)%mo)%mo)%mo;
            if(m%2==0) ans=(ans+km(m,n/2)*(LL)((n/2)%mo)%mo)%mo;
        }
        ans=ans*km(cnt,mo-2)%mo;
        printf("%lld\n",ans);
    }   
}