[BZOJ5330][SDOI2018]反回文串-Pollard-Rho-莫比乌斯反演

反回文串

Description

“回文串什么的最讨厌了……”
小Q讨厌任何形式的回文串：
（1）如果一个字符串从左往右读和从右往左读是一样的，那么小Q讨厌它；例如aa和aba
（2）对于一个字符串来说，若将某个前缀子串移除并拼接到字符串的尾部，能得到一个小Q讨厌的字符串，
那么小Q也会讨厌原来的这个字符串；例如aab和baa。
现在问题来了，如果任意字符串只可以由k种已知的字符组成（也就是说字符集的大小为k），
那么长度为n的所有字符串里，有多少字符串是小Q讨厌的？
答案可能很大，你只需要给出答案对p取模后的值。

Input

第一行包含一个正整数T，表示有T组测试数据。
接下来T行，每行描述一组测试数据，包含三个正整数n,k和p。
1<=T<=10,1<=n<=10^18,1<=k<=n,10^9<=p<=2^30

Output

对于每组测试数据，输出一行，包含一个整数，表示答案对p取模的值。

Sample Input

10
9311702400 2063322107 1032035101
9311702400 3847844404 1014424217
9311702400 6992683534 1067435323
9311702400 1356819643 1024817347
9311702400 6944668829 1042812553
9311702400 9162670852 1003177793
9311702400 6567188538 1062813337
9311702400 3910301217 1018910293
9311702400 2463242636 1018694167
9311702400 4797739673 1050709027

Sample Output

618961590
28597316
1016219343
977847638
108994801
100559666
694084378
492868358
336126742
176095716

---

还以为是个字符串题才进来的……
虽然练一下几乎没写过的pollard-rho也还行

---

思路:
如果只考虑限制1，显然答案就是$k^{\lfloor \frac{n+1}{2} \rfloor}$，即枚举回文串半边的情况。

考虑限制2的话，一个naive的想法是直接乘以$n$，即$nk^{\lfloor \frac{n+1}{2} \rfloor}$。
然而上面的情况显然是有重复的。
比如aaaaaa，这个串就被算了$6$遍。

考虑上面这个错误出现的原因，事实上是因为字符串的循环节小于字符串长度。
如abaabaabaaba这样循环节出现次数为$4$次的会被在上面的方法中算$4$次，诸如此类。

那么又一个简单的想法是，枚举循环节长度，对于某个循环节长度$d$，合法的起点数量是$d$，然后对于每个$d$计算出现次数，并乘以$d$后求和来得出答案。

然而这样还是有问题。
比如这样的串:
abbaabba
baabbaab
可以发现有与循环节无关的额外的轮换方法。

考虑新的反例出现条件，以及出现后的合法起点数量。
可以发现条件只有一个:循环节长度为偶数。
同时，可以发现当出现反例，合法起点数量除以$2$。
于是设$h(x)$为某个循环节长度$x$的合法起点数量，那么有

$$h(x)=\frac{x}{2-x\%2}$$

那么设循环节长度$x$的出现次数为$f(x)$，有

$$ans=\sum_{d | n}f(d)h(d)$$

可以发现$f(x)$不好求，但咱们可以用反演。
反演利用了下面这个等式:

$$k^{\lfloor \frac{n+1}{2} \rfloor}=\sum_{d|n}f(d)$$
由于 $f(d)$代表的是循环节长度最小为 $d$的回文串的数量，那么就可以简单地加起来得到总方案数。
于是来一发莫比乌斯反演，设$g(x)=k^{\lfloor \frac{x+1}{2} \rfloor}$，有:
$$g(n)=\sum_{d|n}f(d)\ =>\ f(n)=\sum_{d|n}g(d)\mu(\frac{n}{d})$$

于是原式可以写成

$$ans=\sum_{d | n}h(d)\sum_{p|d}g(p)\mu(\frac{d}{p})$$

接下来是激动人心的推式子环节

为了方便令$d$的值整体除以$p$，有

$$ans=\sum_{dp|n}h(dp)\sum_{p|dp}g(p)\mu(d)$$
$$ans=\sum_{p|n}g(p)\sum_{d|\frac{n}{p}}h(dp)\mu(d)$$

式子里的$h(dp)$十分难受，考虑提出来。
考虑$h(dp)==dh(p)$成立的条件。
可以发现，如果$dp$为偶数，那么$p$需要是偶数
如果$dp$为奇数，那么什么都不用管。
于是不成立当且仅当$dp$为偶数且$p$为奇数。
此时，显然$\frac{n}{p}$需要是一个偶数。

于是，观察这种特殊情况下，下面式子的值:

$$p\sum_{d|\frac{n}{p}}h(d)\mu(d)$$
上面由于 $p$是奇数，就事先提了出来。
不考虑$\mu$为 $0$的情况，咱们令所有$d$为偶数的项的 $h$值除以2，考虑$\mu$的正负性，加在一起。
因为显然可以令一个 $d$为奇数的$h(d)\mu(d)$与 $-h(2d)\mu(d)$抵消，所以这个式子最后会得到一个 $0$。

然后级可以直接计算，遇到上述情况直接跳过了。
于是继续化简式子得到:

$$ans=\sum_{p|n}g(p)h(p)\sum_{d|\frac{n}{p}}d\mu(d)$$

考虑如何快速算后半部分，令$p(n)=\sum_{d|n}d\mu(d)$，可以发现:

$$p(n)=\prod_{i=1}^{m}(1-p_i)$$
其中 $p_i$为 $n$的质因数分解。

然后$p(n)$显然可以边枚举 $p$的值边算。
于是枚举约数$p$并暴力运算即可~

分解质因数当然是用Pollard-Rho了。
然后写个搜素，枚举每个质因子用几次，就能找到所有约数了。
如果某个质因子使用的个数还有剩余，那么这个数参与对$p(n)$的贡献。

于是便收工了~

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;

typedef long long ll;
const ll L=109;

inline ll read()
{
    ll x=0;char ch=getchar();
    while(ch<'0' || '9'<ch)ch=getchar();
    while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();
    return x;
}

inline ll mul(ll a,ll b,ll p)
{
    ll tmp=(a*b-(ll)((long double)a/p*b+1e-9)*p);
    return tmp<0?tmp%p+p:(tmp>=p?tmp%p:tmp);
}

inline ll qpow(ll a,ll b,ll p)
{
    ll ret=1ll;
    while(b)
    {
        if(b&1ll)ret=mul(ret,a,p);
        a=mul(a,a,p);b>>=1ll;
    }
    return ret;
}

namespace rho
{
    ll stk[L],c;
    int cnt[L],top;

    inline ll gcd(ll a,ll b)
    {
        return (!b)?a:gcd(b,a%b);
    }

    inline bool check(ll a,ll n)
    {
        ll x=n-1;int t=0;
        for(;!(x&1);x>>=1)t++;
        x=qpow(a,x,n);
        if(x==1 || x==n-1)return 1;
        x=mul(x,x,n);
        for(int i=1;i<=t;i++,x=mul(x,x,n))
        {
            if(x==n-1)return 1;
            if(x==1)return 0;
        }
        return 0;
    }

    int test[]={2,3,5,7,11,19,2333};

    inline bool ispri(ll x)
    {
        if(x==1 || x==2)return 1;
        for(int i=0;i<=6;i++)
            if(!check((test[i]%(x-1)+1),x))
                return 0;
        return 1;
    }

    inline bool rho(ll x);

    inline void dfs(ll x)
    {
        if(ispri(x))
        {
            stk[++top]=x;
            cnt[top]=0;
            return;
        }
        while(!rho(x))c=rand();
    }

    inline bool rho(ll x)
    {
        int cnt=0,k=0;
        ll t1=rand()%x,t2=t1,d;
        while(1)
        {
            t1=(mul(t1,t1,x)+c)%x;
            d=gcd(abs(t1-t2),x);
            if(d!=1 && d!=x)
                return dfs(d),dfs(x/d),1;
            if(t1==t2)return 0;
            if((++cnt)==(1<<k))
                k++,t2=t1;
        }
    }

    inline int div(ll x)
    {
        top=0;dfs(x);ll e=top;
        sort(stk+1,stk+top+1);
        top=0;stk[top]=-1;
        for(int i=1;i<=e;i++)
        {
            if(stk[top]!=stk[i])
                stk[++top]=stk[i];
            cnt[top]++;
        }
        return top;
    }
}

using namespace rho;

int pcnt[L],ptop;
ll pv[L],n,k,p,ans;

inline ll h(ll x)
{
    return ((x&1)?x:(x>>1ll))%p;
}

inline bool dfs(int x,ll d,ll t)
{
    if(x==ptop+1)
    {
        if((d&1) && !((n/d)&1))return 0;
        return (ans+=mul(mul(qpow(k,(d+1)>>1,p),h(d),p),t,p))%=p,0;
    }
    for(ll i=0,v=1;i<=pcnt[x];i++,v*=pv[x])
        dfs(x+1,d*v,t*((i<pcnt[x])?(p+1-pv[x]%p):1)%p);
    return 0;
}

int mina()
{
    n=read();k=read();p=read();k%=p;
    c=rand();ptop=div(n);
    for(int i=1;i<=ptop;i++)
        pv[i]=stk[i],pcnt[i]=cnt[i];
    ans=0;dfs(1ll,1ll,1ll);
    printf("%lld\n",ans);
    return 0;
}

int main()
{
    srand(514);
    for(int T=read();T;T--)
        mina();
    return 0;
}