题解 P4607 【[SDOI2018]反回文串】

反演神题…… 首先你得看出来这道题是一个反演题 如果你这道题卡94了检查一下你的Miller-Rabin有没有在判定大质数的时候返回奇奇怪怪的值…… _______________ ## 本题题解 题目的要求是求下列符合条件的字符串个数(长度为n,字符集大小为k) ## 这个串的一个轮换是一个回文串 两个字符串不同当且仅当这个字符串中有一个位置的字符不同 然后此时我们看了一样数据范围$n<=10^{18}$非常魔幻 这已经不是dp可以解决的问题了，所以让我们来考虑一下如何列式子 ### 找规律 首先我们想一个非常trival的问题，长度为n字符集大小为k的字符串集里有多少个回文串？显然是 ## $k^{\lfloor\frac{n+1}{2}\rfloor}$ 但是我们的直觉告诉我们答案显然不是 ## $nk^{\lfloor\frac{n+1}{2}\rfloor}$ 比如像aaaaaa……这种串的所有轮换都是同一个串，然后显然是不可以算n遍的 那么我们说我们刨掉这种情况不就好咯，然后我们发现有这样的反例 abc abc abc abc abc abc 轮换3次之后陷入了循环还是不能算n遍 那么我们考虑一件事情，一个串到底有多少个互不相同的轮换？ 发现刚好等于它的最小循环节长度(这里要求的最小循环节是不允许重叠的，换句话说最小循环节必须是这个串的约数) 那么我们是不是可以枚举每一个回文串的最小循环节长度d，然后每个回文串对答案产生d的贡献呢？ 然后我们发现会有这样的反例 abcd abcd abcd abcd 轮换两次之后变成了 cd abcd abcd abcd ab 我们惊讶的发现这个串也是一个回文串，换句话说"abcdabcdabcdabcd"和"cdabcdabcdabcdab"这两个串共同产生了4的贡献，因此每个串平摊下来仅仅产生了2的贡献 进一步思考，如果最小循环节长度d为偶数的话，我们会发现在它的d种不同的轮换字符串当中恰好有两个或者零个回文串，因此我们枚举的回文串最小循环节长度d如果是偶数的话仅仅会产生$d/2$的贡献，因为另一半被和它"配对"的回文串摊掉了 然后让我们来想一下d为奇数的时候会出现什么，万幸的是在它的d种不同的轮换字符串中最多有一个回文串，所以如果d为奇数那么刚好产生d的贡献 那么我们似乎可以用一个分段函数h来描述我们最小循环节长度为d的字符串对答案的贡献 ## $h(d)=(dmod2)?d:d/2$ 那么我们现在终于可以愉快的枚举循环节长度d了，什么？怎么求最小循环节长度为d的回文串数量?别急，让我们设最小循环节长度为d的回文串数量为$f(d)$ 那么我们现在终于可以列出式子来了 ## $Ans=\sum_{d|n}h(d)f(d)$ 问题来了，怎么求$f(d)$?你不会求我也不会求…… 但是万幸的是这里有一个比较显然的恒等式($g(m)$代表长度为m的回文串总数) ## $g(m)=k^{\lfloor\frac{m+1}{2}\rfloor}=\sum_{d|m}f(d)$ 也就是说我们枚举了长度为m的所有回文串的最小循环节,然后把每种回文串的数量加到了一起(显然两个串最小循环节长度不同那么这两个串不同，因此最小循环节为d的回文串属于互斥事件，用加法原理进行合并)那么我们就会得到长度为m的所有回文串的数量了 (下面开始反演了，如果还对莫比乌斯反演不是相当熟练的话还先去学习下反演再来看这道题吧) 如果你对狄利克雷卷积那一套相当熟练的话你会立即想到 ## $g(m)=\sum_{d|m}f(d)1(m/d)$ 其中$1()$是常函数,1(x)=1，其中x是任意值 也就是说 ## $g=f×1$ 那么莫比乌斯反演公式告诉我们 ## $g=f×1 \leftrightarrow f=g×\mu$ 换句话说我们可以推出 ## $f(m)=\sum_{d|m}g(d)\mu(m/d)$ 好了那么我们现在可以计算f函数的值了，那么我们把它代入原来答案的表达式会发生什么…… ## $Ans=\sum_{d|n}h(d)\sum_{p|d}g(p)\mu(d/p)$ 好了下面是神奇的交换$\Sigma$的环节了！ 我们令$d/=p$,则有$d=dp$,因为$d|n$所以$dp|n$所以$d|n/p$ 我们接下来就全部用$dp$替换原来式子当中的d,得到了 ## $Ans=\sum_{p|n}g(p)\sum_{d|\frac{n}{p}}h(dp)\mu(d)$ 那么我们到这一步其实已经不可以反演了 所以让我们来考虑一下h函数的性质，我们尝试着把$h(dp)$中的$d$提出来 换句话说，我们考虑这个式子在什么条件下成立 ## $h(dp)=dh(p)$ 发现当且仅当d为偶数且p为奇数的时候这个式子不会成立 但是让我考虑下什么时候d为偶数且p为奇数呢？ d必须是$\frac{n}{p}$的约数，然而如果$\frac{n}{p}$是奇数的话d不可能是偶数 那么这个式子唯一不成立的情况就是p为奇数且$\frac{n}{p}$为偶数了 那么我们考虑此时这个式子的值是多少 ## $\sum_{d|\frac{n}{p}}h(dp)\mu(d)$ 由于p是奇数，所以可以写成 ## $p\sum_{d|\frac{n}{p}}h(d)\mu(d)$ 那么我们考虑后边式子展开之后的项,$\mu$值为0的项我们不管，考虑$\mu$值为1和-1的项 那么我们对d是奇数的项把它乘2，偶数的的项把它除2 那么显然可以把这些$\mu$值为1和-1的项两两配上对 那么我们发现乘2或者除2的过程中$h(d)$的值并没有变而$\mu(d)$的值刚好取反 因此每个对的和恰好是0 所以我们证明了什么？刚才的式子的和为0，因此这种情况遇到之后直接continue掉就好了 那么剩下的情况我们认为$h(dp)=dh(p)$这个式子都是成立的 那么式子可以化简成 ## $Ans=\sum_{p|n}g(p)h(p)\sum_{d|\frac{n}{p}}d\mu(d)$ 那么我们记$t(m)=\sum_{d|m}d\mu(d)$ 那么我们看一看它等于什么，还是老样子我们考虑这个式子的展开式，$\mu$值为0的我们不用管，考虑1和-1的项 那么我们可以把m质因数分解，于是上边的式子和下边的式子应该是等价的 ## $t(m)=\prod_{i=1}^{q}(1-p_{i})$ 其中$p_{1}……p_{q}$是m的质因数分解 因为我们把下面的式子展开之后应该每一项是一堆p乘起来，而如果这个项里有奇数个质数的话那么系数就是-1反之就是1，那么刚好和莫比乌斯函数的定义对应起来，而这些质数乘起来的值恰好就是d,由于每个质因子只出现一次因此$\mu$值为0的项不会出现，所以上边和下面的式子是等价的 所以最后我们的式子就变成了 ## $Ans=\sum_{p|n}g(p)h(p)t(n/p)$ 最后一个问题，如何生成一个数的所有约数? 我们可以把这个数进行质因数分解，那么就可以dfs出所有的约数了 问题来了怎么质因数分解？ ___________________ ### 前置芝士:Pollard-Rho算法 期望复杂度$O(\sqrt{\sqrt{n}})$的神奇算法 简单来说，我们不停的随机两个数，然后取他们的差，和要分解的数取gcd，如果gcd不是1的话递归分解gcd和n/gcd，直到遇到了质数，此时我们就可以把这个质因数放到set里了 但是其实我们写不出随机函数来，我们的随机函数是用一些奇奇怪怪的式子生成的，比如$x=(x^2+1)modn$这个东西最后一定会陷入死循环，此时如果我们还是没有猜对我们就凉了 那么我们就考虑如何判掉这个环，一个比较古老的方式是Floyd2倍速判断法，这里介绍一个比较机智的判环法，Brent判环法 我们先钦定一个初始的种子x1 然后计算$x_{2}=x_{1},x_{1}=(x_{1}^{2}+1)modn$ 然后我们不停的计算$x_{i}$如果i是2的整次幂那么令$x_{2}=x_{1}$ 每次返回$x_{2}-x_{1}$作为我们函数的返回值，当返回值是0的时候表明出现了环直接重新选择种子 听起来挺玄学的但是你自己画一个$\rho$型的环你会发现这个东西真的可以把环判出来 然后有人证明了这个东西的复杂度是期望$O(\sqrt{\sqrt{n}})$的，于是我们就可以对这个数做质因数分解了 ## 实现 我们dfs出所有的约数值，t函数的值可以在dfs的时候顺便求出，然后我们直接按照推出来的式子枚举约数暴力计算即可，t函数的值用unordered_map存一下就行了，Pollard-Rho的质数判定采用Miller-Rabin随机判定法即可 最后一个小问题，我们可能要做膜数是longlong范围内的乘法，该怎么办呢？ ### 强行转成long double手动模拟膜法即可/int_128也可 上代码~ ```C // luogu-judger-enable-o2 #include<cstdio> #include<algorithm> #include<set> #include<cstdlib> #include<tr1/unordered_map> using namespace std::tr1; using namespace std;const int N=1e5+1e4+10;typedef long long ll;typedef long double ld; ll n;ll k;ll mod;set <ll> s;int T;ll tims; inline ll mul(ll a,ll b,ll md){ll tp=a*b-(ll)((ld)a/md*b+1.0e-8)*md;return tp<0?tp+md:tp;}\\转longdouble的乘法 inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;}//快速幂 namespace MiiR//miiler-rabin { const ll bs[20]={2,3,5,7,11,61,24151};//偷懒少选了几个base inline ll po(ll a,ll p,ll mod){ll r=1;for(;p;p>>=1,a=mul(a,a,mod))if(p&1)r=mul(r,a,mod);return r;} inline bool tst(ll x,int t)//二次探测 {for(ll p=x-1;;p>>=1){ll g=po(bs[t],p,x)%x;if(g==x-1)return 0;if(g!=1)return 1;if(p&1)return 0;}} inline bool ck(ll x){for(int i=0;i<=6&&bs[i]<x;i++){if(tst(x,i))return false;}return true;} } namespace PollR { ll st; inline ll mabs(ll x,ll y){return (x>y)?x-y:y-x;} inline ll gcd(ll a,ll b){ll c;while(b){c=a%b;a=b;b=c;}return a;} inline ll rd(ll x,ll md){return (mul(x,x,md)+st)%md;} inline ll bigr(){ll ret=0;for(int i=0;i<=62;i++)ret+=(1LL<<i)*(rand()%2);return ret;}//这里写了个暴力rand，大家不要学我 inline void solve(ll x)//递归分解，质数直接返回 { if(x==1)return;if(MiiR::ck(x)){s.insert(x);return;} while(1) { ll x1=bigr()%x;st=rand()%x+1; ll g=gcd(x,x1);if(g!=1&&g!=x){solve(g);solve(x/g);return;} ll x2=x1;x1=rd(x1,x);ll k=0;ll tr=1; for(;x1-x2;x1=rd(x1,x),k++)//brent判圈法 { g=gcd(x,mabs(x1,x2));if(g!=1&&g!=x){solve(g);solve(x/g);return;} if(k==tr){x2=x1;tr<<=1;} } } } } ll ys[N];int hd; unordered_map <ll,ll> fuc;//存t函数的map void dfs(ll d,ll f,set <ll>::iterator it)//dfs找质数 { if(it==s.end()){ys[++hd]=d;fuc[d]=f;return;} set <ll>:: iterator it1=++it;--it;dfs(d,f,it1); for(ll mi=*it;;mi*=(*it))//注意别爆longlong {dfs(d*mi,f*(mod-*it%mod+1)%mod,it1);if((n/mi)%(*it)!=0)break;} } inline void solve() { scanf("%lld%lld%lld",&n,&k,&mod);ll ret=0;k%=mod;//小心爆longlong PollR::solve(n);dfs(1,1,s.begin());//分解之后找约数 for(int i=1;i<=hd;i++)//枚举约数暴力计算即可 { ll d=ys[i];if(d%2==1&&(n/d)%2==0){continue;} (ret+=((d%2)?d:d/2LL)%mod*po(k,(d+1)/2)%mod*fuc[n/d]%mod)%=mod; }printf("%lld\n",ret); } inline void clear(){set <ll> emp;unordered_map <ll,ll> ep;swap(s,emp);swap(fuc,ep);hd=0;tims=0;}//O(1)的clear int main(){srand(66623366);scanf("%d",&T);for(int z=1;z<=T;z++){solve();clear();}return 0;}//拜拜程序~ ```