Luogu P4607 [SDOI2018]反回文串

Link
一个字符串\(S\)的不同的轮换数为\(|\sigma(S)|\)，其中\(\sigma\)表示最小循环节。
直接计算\(\sum|\sigma|\)会把长度为偶数的循环节算重，设\(f(d)=\sum\limits[|\sigma|=d]\)，那么\(ans=\sum\limits_{d|n}\frac{f(d)d}{1+[2|d]}\)。
设\(g(n)\)表示长度为\(n\)的回文串的数量，那么\(g(n)=k^{\lfloor\frac n2\rfloor}\)。
显然\(g=f\times I\)，因此\(f=g\times\mu\)，记\(h(d)=\frac d{1+[2|d]}\)那么

\[\begin{aligned} ans&=\sum\limits_{d|n}(h\cdot(\mu\times g))(d)\\ &=\sum\limits_{d|n}g(d)\sum\limits_{t|\frac nd}\mu(t)h(dt) \end{aligned} \]

当\(2\nmid d\wedge2|t\)时，不难发现\(\sum\limits_{t|\frac nd}\mu(t)h(dt)=0\)。
否则我们有\(h(dt)=th(d)\)，因此

\[\begin{aligned} ans&=\sum\limits_{d|n}g(d)\sum\limits_{t|\frac nd}\mu(t)h(dt)\\ &=((g\cdot h)\times(\mu\cdot id)\times I)(n) \end{aligned} \]

考虑\(((\mu\cdot id)\times I)(n)\)，若\(n=\prod p_i^{e_i}\)，分析贡献不难发现\(((\mu\cdot id)\times I)(n)=\prod(1-p_i)\)。
那么我们dfs枚举\(n\)的约数时维护\(((\mu\cdot id)\times I)(n)=\prod(1-p_i)\)就可以\(O(d(n)\log n)\)地计算\(((g\cdot h)\times(\mu\cdot id)\times I)(n)\)了。

#include<map>
#include<cctype>
#include<cstdio>
#include<random>
#include<algorithm>
using i64=long long;
int tot;i64 n,k,ans,P;
struct prime{i64 p;int e;}a[20];std::map<i64,int>cnt;
std::mt19937 rng(19260817);
i64 read(){i64 x=0;char c=getchar();while(isspace(c))c=getchar();while(isdigit(c))(x*=10)+=c&15,c=getchar();return x;}
void inc(i64&a,i64 b){a+=b-P,a+=a>>63&P;}
i64 mul(i64 a,i64 b,i64 p){return (__int128)a*b%p;}
i64 add(i64 a,i64 b,i64 p){return a+=b-p,a+=a>>63&p;}
i64 pow(i64 a,i64 b,i64 p){i64 r=1;for(;b;b>>=1,a=mul(a,a,p))if(b&1)r=mul(a,r,p);return r;}
int Judge(i64 x)
{
    static int p[7]={2,3,5,7,11,61,24151};
    for(int i=0;i<7;++i) if(x==p[i]) return 1; else if(!(x%p[i])||pow(p[i],x-1,x)!=1) return 0;
    i64 phi=x-1;while(~phi&1) phi>>=1;
    for(int i=0;i<7;++i)
    {
	i64 y=pow(p[i],phi,x);
	if(y==1) continue;
	while(y!=1&&y!=x-1) y=mul(y,y,x);
	if(y==1) return 0;
    }
    return 1;
}
i64 Find(i64 x)
{
    i64 c=rng()%(x-1)+1,las=0,now=0;
    for(int len=1;;las=now,len<<=1)
    {
	i64 prod=1,d;
	for(int i=1;i<=len;++i) if(now=mul(now,now,x)+c,prod=mul(prod,llabs(now-las),x),!(i&127)&&(d=std::__gcd(prod,x))!=1) return d;
	if((d=std::__gcd(prod,x))!=1) return d;
    }
}
void Factor(i64 x,int f)
{
    if(x==1) return ;
    if(Judge(x)) return cnt[x]+=f,void();
    i64 p=Find(x);int c=0;
    while(!(x%p))x/=p,++c;
    Factor(p,f*c),Factor(x,f);
}
void dfs(int c,i64 now,i64 val)
{
    if(c==tot+1)
    {
	if(now&1&&~(n/now)&1) return ;
	return inc(ans,mul(mul(pow(k,(now+1)/2,P),now/(now&1? 1:2),P),(val%P+P)%P,P));
    }
    i64 x=1;for(int i=0;i<=a[c].e;x*=a[c].p,++i) if(i==a[c].e) dfs(c+1,now*x,val); else dfs(c+1,now*x,val*(1-a[c].p));
}
void solve()
{
    cnt.clear(),ans=tot=0,Factor(n=read(),1),k=read(),P=read(),k%=P;
    for(auto it:cnt) a[++tot]={it.first,it.second};
    dfs(1,1,1),printf("%lld\n",ans);
}
int main()
{
    for(int t=read();t;--t) solve();
}