2019.03.07【SDOI2018】【洛谷P4607】【BZOJ5330】反回文串（莫比乌斯反演）（Pollard-Rho）

洛谷传送门

BZOJ传送门

---

解析：

首先稍有常识的人都知道这道题绝对不可能是字符串题。

这种只给串长和字符集大小的题目只可能是计数。

而计数方式有很多啊，DP，群论，生成函数，甚至这道题的做法，莫比乌斯反演。

看数据$DP$和群论基本上就告别正解了，生成函数推不出通项公式也是白搭（而且字符串轮换循环的情况让生成函数也告别正解了）。

这TM谁一眼看得出是莫反啊。

我们先忘记我们知道这道题需要用莫比乌斯反演这件事，来尝试推理一遍。

什么？你不会莫比乌斯反演？敢来淦反演毒瘤省SDOI的题你说你不会莫比乌斯反演？

没办法了，洗洗睡吧

---

推♂倒：

我们需要求有多少个长度为$n$，字符集大小为$k$的字符串的某一个轮换是一个回文串。

首先，字符集大小为$k$，长度为$n$的总回文串个数，很显然是$k^{\lfloor \frac{n+1}{2}\rfloor }$

但是我们也很清楚，并不是所有回文串的每一个轮换都不同，比如$aaaaa$就只能算作一个串。

考虑每个回文串有多少种不同的轮换，答案显然就是最小循环节长度。

但是每个回文串的所有轮换都要考虑一次吗？

我们发现$abbaabba$的所有轮换和$baabbaab$都相同。。。

也就是说，可能出现某个回文串的某个轮换是另一个不同的回文串。

日，感觉分析不下去了。

但是仔细一看，这种情况只有当回文串的循环节长度为偶数的时候才会发生。就是循环节转到一半的情况。

而且只可能是两个回文串平摊贡献，不然循环节长度还能再小。

换句话说，如果考虑所有回文串的不同轮换对答案的贡献，设某个回文串的最小循环节长度为$d$，则设贡献函数$h(d)=(d\%2==0)?d/2:d$

现在令$f(d)$表示所有最小循环节长度为$d$的回文串个数，则答案为$$Ans=\sum _{d\mid n}h(d)f(d)$$

那么现在问题变成了求$f$。

注意这里要求是最小循环节，即不能有更小的循环节了。

日，感觉分析不下去了。

但是你回去想一想为什么我们会推到最小循环节就会发现：

$$g(N)=\sum _{d\mid N}f(d)=k^{\lfloor \frac{N+1}{2}\rfloor }$$

其中$g(N)$表示具有长度为$N$的循环节的回文串总数。

你知道我要说什么吧

莫比乌斯反演！

$$f(N)=\sum _{d\mid N}\mu (d)g(\frac{N}{d})$$

代回去得到：

$$\begin{array}{rlrl}Ans& =& & \sum _{d\mid n}h(d)\sum _{t\mid d}\mu (d)g(\frac{d}{t})\\ & =& & \sum _{d\mid n}g(d)\sum _{td\mid n}h(td)\mu (t)\end{array}$$

我们现在想把$h(td)$中的一个常数$d$提出来得到$th(d)$，因为${\sum }_{t\mid \frac{n}{d}}t\mu (t)$有一个很好用的性质，一会儿会用到。

现在考虑这个式子$h(td)=th(d)$什么时候不成立

很显然就是$t$为偶数但是$d$为奇数数的时候。但是这时候可以把$d$提出来。

考虑这时候这个式子的值是什么：$$d\sum _{t\mid \frac{n}{d}}h(t)\mu (t)$$

我们发现对于所有的偶数$t$，如果$\mu (t)!=0$，都有$h(t)\mu (t)=-h(\frac{t}{2})\mu (\frac{t}{2})$。

换句话说，$\sum$里面的所有奇数和偶数凑成一对后和为$0$。

特判掉所有这种情况就行了。

去掉上面的情况后我们直接来提出$t$：
$$Ans=\sum _{d\mid n}g(d)h(d)\sum _{t\mid \frac{n}{d}}t\mu (t)$$

设$S(n)={\sum }_{t\mid n}t\mu (t)$，这个函数有显然的递推式，设$n$的唯一分解为$n=\prod _{i=1}^t{p}_i^{k_i}$，则$S(n)=\prod _{i=1}^n(1-p_i)$

利用Pollard-Rho筛出所有质因子，dfs找出所有因数和$S$函数值就行了。

---

代码（很丑，丑到我都不想看）：

#include<bits/stdc++.h>
#include<tr1/unordered_map>
using namespace std;
#define ll long long
#define re register
#define gc get_char
#define cs const

inline ll mul(cs ll &a,cs ll &b,cs ll &mod){
	return (a*b-(ll)((long double)a/mod*b)*mod+mod)%mod;
}

inline ll quickpow(ll a,ll b,cs ll &mod){
	ll res=1;
	while(b){
		if(b&1)res=mul(res,a,mod);
		a=mul(a,a,mod);
		b>>=1;
	}
	return res;
}

inline ll gcd(ll a,ll b){
	static ll tmp;
	while(b){
		tmp=a%b;
		a=b;
		b=tmp;
	}
	return a;
}

cs int P=1e7+7;
int prime[P],pcnt,minpri[P];
bool mark[P];
inline void linear_sieves(int len=P-7){
	mark[1]=true;
	for(int re i=2;i<=len;++i){
		if(!mark[i])prime[++pcnt]=i,minpri[i]=i;
		for(int re j=1;i*prime[j]<=len;++j){
			mark[i*prime[j]]=true;
			minpri[i*prime[j]]=minpri[i];
			if(i%prime[j]==0)break;
		}
	}
}

inline bool isprime(ll x){
	if(x<=P-7)return !mark[x];
	ll t=x-1,s=0;
	while(!(t&1))t>>=1,++s;
	for(int re i=1;i<=5;++i){
    	ll p=prime[rand()%pcnt+1]%x;
        ll num=quickpow(p,t,x),pre=num;
        if(x%p==0)return false;
        for(int re j=0;j<s;++j){
            num=mul(num,num,x);
            if(num==1&&pre!=1&&pre!=x-1)return false;
            pre=num;
        }
        if(num!=1)return false;
    }
    return true;
}

inline ll Pollard_Rho(ll x){
	if(x%2==0)return 2;
	if(x%3==0)return 3;
	if(x%5==0)return 5;
	if(x%7==0)return 7;
	ll n=0,m=0,q=1,t=1,c=rand()%(x-1)+1;
	for(ll re k=2;;k<<=1,m=n,q=1){
		for(ll re i=1;i<=k;++i){
			n=(mul(n,n,x)+c)%x;
			q=mul(q,abs(n-m),x);
		}
		if((t=gcd(q,x))>1)return t;
	}
}

ll pfact[60];
int pfcnt;

inline void sieve(ll x){
	if(x==1)return ;
	if(x<=P-7){
		while(x>1){
			pfact[++pfcnt]=minpri[x];
			int p=minpri[x];
			while(x%p==0)x/=p;
		}
		return ;
	}
	if(isprime(x)){
		pfact[++pfcnt]=x;
		return ;
	}
	ll p=x;
	while(p==x)p=Pollard_Rho(p);
	sieve(p);
	while(x%p==0)x/=p;
	sieve(x);
}

ll n,k,mod;
ll fact[P];
int fcnt;
tr1::unordered_map<ll,ll>f;

inline ll quickpow(ll a,ll b){ll res=1;while(b){if(b&1)res=mul(res,a,mod);a=mul(a,a,mod);b>>=1;}return res;}
inline ll mul(cs ll &a,cs ll &b){return (a*b-(ll)((long double)a/mod*b)*mod+mod)%mod;}
inline void dfs(ll now,ll nowf,int pos){
	if(pos>pfcnt){
		fact[++fcnt]=now,f[now]=nowf;
		return ;
	}
	dfs(now,nowf,pos+1);
	for(ll p=pfact[pos];;p*=pfact[pos]){
		dfs(now*p,mul(nowf,(mod-pfact[pos]%mod+1)%mod),pos+1);
		if((n/p)%pfact[pos])return ;
	}
}

inline void solve(){
	scanf("%lld%lld%lld",&n,&k,&mod);k%=mod;
	ll ans=0;
	sieve(n);
	sort(pfact+1,pfact+pfcnt+1);
	pfcnt=unique(pfact+1,pfact+pfcnt+1)-pfact-1;
	dfs(1,1,1);
	for(int re i=1;i<=fcnt;++i){
		ll d=fact[i];
		if((d%2)==1&&(n/d)%2==0)continue;
		ans=(ans+mul(f[n/d],mul(quickpow(k,(d+1)/2),(d%2?d:d/2))))%mod;
	}
	cout<<ans<<"\n";
}

inline void init(){
	pfcnt=0;
	fcnt=0;
	f.clear();
}

int T;
signed main(){
	srand(time(0));
	linear_sieves();
	scanf("%d",&T);
	while(T--){
		init();
		solve();
	}
	return 0;
}