合法字符串——庞果网

最近一直很忙，希望我能够写完。

这道题目并不难，却很复杂。只要你有耐心，根据题解，一定能够做出来。唯一的困难就是，你常常半途而废。能够看完题解的，应该是一个有耐心的人。

这个题目应该是庞果网上目前为止最难的一道题目了，今天在群中大神FancyMouse的提点下，找到一篇名为“Recursive Sequence”的文章，完成了m项递推式的O(m^2logn)时间复杂度的优化，这是这道题目的最终的优化，应该说，完成这个优化，这道题就没有可以优化的地方了，当然，边边角角的优化也是可以提高效率，但终究起不了大作用。

题目将逆元的运用、多项式系数计算、m项递推式优化这三个我认为主要的知识点融合在了一起，不能不说，这是一道很好的题目（个人认为）。

我的想法还是根据群里给出的题解，但是那个写的有些乱，又不详细，喜欢看的人很少，我将按照自己对题解的梳理，叙述这道题目的解法。

题目详情：

用n个不同的字符（编号1 - n），组成一个字符串，有如下2点要求：

1、对于编号为i的字符，如果2 * i > n，则该字符可以作为最后一个字符，但如果该字符不是作为最后一个字符的话，则该字符后面可以接任意字符；

2、对于编号为i的字符，如果2 * i <= n，则该字符不可以作为最后一个字符，且该字符后面所紧接着的下一个字符的编号一定要大于2 * i。

问：有多少长度为M且符合条件的字符串。

例如：N = 2，M = 3。则abb，bab，bbb是符合条件的字符串。剩下的均为不符合条件的字符串。

输入：N，M（2 <= N，M <= 1 * 10^9）

输出：满足条件的字符串的个数，由于数据很大，输出该数Mod 10^9 + 7的结果。

函数头部：int validString(int n , int m){}

题目的叙述简单，但是题目的难度却很大，而难度大的原因，来自于字符的数量和字符创的长度都是十的九次方。这个数量很庞大，为此，我们需要激昂庞大的问题拆解成能够高效完成的小问题。根据题解，不断解决一个小的模块，直到解决整个问题。

理论基础

题解首先给出的是一个定义。【合法链】：从一个字母开始，后边的每一个字母编号都大于等于前一个的2倍，这样尽可能的连接字母，直到不能再连接，即最后一个字母的编号满足2 * i > n。这样，后边的字母必定不满足大于等于前一个字母编号的2倍，就不是这个【合法链】的一部分了。

有了这个定义，我们可以总结如下规律：一个合法的字符串是有若干【合法链】组成的；对于每一个字符串，划分成【合法链】的方法是唯一的。

关于规律我不给出证明，但它是正确的。简单说一下，我们把字母分类，所有编号为i，且2 * i <= n的字母归为A类，其余的归为B类。B类字母的编号i满足2 * i > n。任意一个【合法链】只有一个B类字母，其余的全部是A类字母。任意一个合法字符串，我们按照B类字母将字符串断开，一个B类字母之前的所有A类字母，都属于这个B类字母，这样，一个字符串就划分成【合法链】，这种划分是唯一的，因为，如果从任何一个合法链断开，都不能构成两个【合法链】。于是，一个合法字符串是由若干个【合法链】组成的，我们可以通过拼接【合法链】生成一个合法字符串。

将问题转化为两步：

（1）、G(len)表示长度为len的【合法链】的个数，求出G(len)。

（2）、V(len)表示长度为len的合法字符串个数，求V(len)。

通过以上定义，我们可以很容易的知道，V(len) = G(1)*V(len - 1) + G(2)*V(len - 2) + ......+ G(p)*V(x - p)；其中，p是最长的【合法链】的长度。对于字符数量n来说，最长【合法链】的长度是多少呢？比如，n = 10 的时候，最长的为1248，1249，124(10)，125(10)，最大长度为4，总共有四条。可以总结出，p =[ log(n)] + 1。（注：所有的log(n)在本文中均表示已2为底，不会打下取整的符号，暂时以[]来代替，下同）。

对于步骤（1）这个比较难计算，对于步骤（2）比较好计算，仅仅使用到矩阵快速幂与矩阵乘法的平方复杂度的优化，理论与实现都很简单。在此，步骤（2）的理论叙述我将放到后边，先解决最困难的步骤（1），并给出相应的方法。

步骤（1）的计算技巧：

对于步骤（1），题解中给出了解决的方法，但这个方法有点复杂，我没有采用。对此，群里给出一个相对简单的解决办法。仍然先给出一个定义：【伪合法链】：从一个字母开始，后边每一个连接的字母编号必须大于等于前边的2倍。

名字是我起的，有点难听，不过，数学中也有类似的称呼，并非带有贬义，而是这个叫法能够很好地描述出改定义的一些信息：它不是真正的【合法链】，但它是还没有变成【合法链】的“【合法链】”，或者可以称之为【合法链】的前身。

我们再次观察定义，发现【伪合法链】的定义中，没有了最后一个字符是B类字符的要求，但是其它条件的还是跟【合法链】一样。它有一个特性：首先定义一下，本文中定义G(len)为长度为len的【合法链】个数，那么，再定义g(n , len)为字母个数为n的时候【伪合法链】的个数。于是，有下面一个关系：当字符个数为n的时候，G(len) = g(n , len) - g(n / 2 , len)。这里的n / 2是需要下取整的。

这个关系式的证明，我没有做，只要知道，它是正确的就可以了。下面开始计算【伪合法链】的数量。这跟计算【合法链】的数量基本相同。

假设f(len , x)为一字母x开头的，长度为len的【伪合法链】的个数。于是，有f(len , x) = f(len - 1 , x * 2) + f(len - 1 , x * 2 + 1) + ... + f(len - 1 , n) 。

这样，g(n , len) = f(len , 1) + f(len , 2) + f(len , 3) + ...... + f(len , n)；为了能够直观的理解，以n = 7 为例。

这应该是最简单的一类求f(len , x)的示例了。通过它，我们可以总结出一些规律：

1、对于固定的n，f(len , x)是中的len大为p = [ log(n)] + 1，对于你n = 7，最大值为3。

2、f(len , x)是一个多项式，而要求g(len)，需要对这个多项式求和。

3、对于固定的len，多项式f(len , x)中的x的取值范围并不是固定的1到n，而是，每次都要除以2，下取整。比如，图中的f(2 , x)，取值范围为1到3，其中的3就是7 / 2的结果，再比如，f(3 , x)，其中的x，取值范围为1到1，其中第二个1就是3 / 2下取整的结果。

4、这条比较重要，我用一张图表示：

看到这里，应该能够想到求取g(len)的方法了。对于这个问题，我用最笨的办法，按部就班，没有任何技巧的求出了所有的g(len)。然而，事情还没有完结，因为想要求出多项式，也就是求多项式的系数，你还需要做几项准备。

第一项：求出单项的系数，这个单项求出来的是关于x的多项式。而求这个单项的系数，就需要先求出自然数1到n的i次幂和的系数，这个i次幂和的最高幂次为i+1，幂和的表达式为n的多项式，将2x - 1替换掉n，就可以得出单项的多项式，一个关于x的多项式。这个多项式再乘以之前的系数ki，就可以添加到f(len , x)的系数中。求取所有单项，再加上g(len - 1)，就可以求出新的f(len , x)的多项式的系数。如此递推下去，可以求出所有的f(len , x)，顺便也就有了g(len)。

第二项：等一下，还有如何用2x - 1替换掉n，求出n的几次幂的系数来，我也是写的程序，这个不难，在这里不做讲述。

第三项：如何求出自然数1到n的i次幂和的系数呢？我说过，这是个关于n的i + 1次多项式，最高次数i + 1。对此我才用了网上的一种方法，由于篇幅有限，仅仅把关键的地方贴出来，我看的网上的论文是：《自然数k次方幂和的一种简捷算法》。从百度文库中可以找到。下面是关键的算法：

不多叙述我就是据此求出所有系数的。

第四项：解决减法运算、分数的方法：这个方法就是运用逆元。关于逆元我也不做解释，可以百度之，有一大堆资料。这里数一下为什么要解决减法运算和分数，关于分数，我们是不好存储的，由于题目中是求解模掉10^9 + 7的结果，所以，运用10^9 + 7的对于各个数值的逆元，能够让所有的数都变成整数，便于我们取模和运算。而对于减法，由于取模的缘故，也是不好计算的，不如全部都编程加法，求取逆元就可以了。下面说一下逆元运用到的一些性质（其中x、y、z都是整数，以gbc(x)代表x关于10^9 + 7的逆元，设MOD = 10^9 + 7）：

（1）、x * （y / z） mod MOD = x * y * gbc(z) mod MOD；

（2）、(x + y) mod MOD = (x mod MOD + y mod MOD) mod MOD；

这两个式子基本在取模运算中够用了。下面是我的求逆元运算的函数：

/** * 扩展欧几里德算法:求a模b的逆元 * @param a * @param b * @return */public static long extendGCD(long a, long b) {long res = b;long x1 = 1, x2 = 0;long y1 = 0, y2 = 1;long x = 0, y = 0;long r = a % b;long q = (a - r) / b;while (r > 0) {x = x1 - q * x2;y = y1 - q * y2;x1 = x2;y1 = y2;x2 = x;y2 = y;a = b;b = r;r = a % b;q = (a - r) / b;}if (x < 0) {res += x;} else {res = x;}return res;}

有关运算中用得到的还有求取组合数，网上已经有很多解答，感觉篇幅很大了，不再赘述。

有了求逆元的算法，有了以上四项求多项式的规律，应该可以求出【伪合法链】的数量g(n , len)了，于是利用公式G(len) = g(n , len) - g(n / 2 , len)。这里需要求数量为n的【伪合法链】的数量，还要求字符数量为n / 2的【伪合法链】的数量。我没有研究过这种求法会不会降低效率，或许直接求【合法链】数量中求两个多项式比这个更复杂，就不得而知了。

到此，我们完成了开头讲到的步骤（1）的全部工作，最难得已经完成了。

步骤（2）的计算技巧：

步骤（2）是比较容易计算的。它就是m项递推式，回顾一下开头给出的公式：V(len) = G(1)*V(len - 1) + G(2)*V(len - 2) + ......+ G(p)*V(x - p)；我们如下构造矩阵：

这就是p项递推式的矩阵。右边的初始向量，我们可以先用公式求出来，这样，乘以一次中间的矩阵，就可以获得V(p)，继续下去，可以求得任意的V(m)。这里求矩阵需要用到快速幂递推，快速幂我并不做介绍，因为这个网上找到很容易。单纯使用快速幂，其复杂度为O(p^3 * log(m))。对于庞果网上的题目，时间复杂度已经足够低了，但是如果你要挑战51nod上的同样的题目”字符串数量V3“，那么，你还需要使用矩阵相乘的优化，使得时间复杂度成为O(p^2 * log(m))。

而矩阵相乘的优化，只做简单叙述，我在开始提到一篇文章，优化的详细方法在里边。假设中间的矩阵为M，那么，我们可以很容易的知道M的特征多项式

根据Cayley-hamilton定理，可以知道

注意，这个公式可以将M的k+i次幂变成次数比它小的一系列M的方幂的线性组合。从而，我们可以得到一个结论：M的任意次方，我们都能够将它表示成E、M^1、M^2、M^3......M^p-1的线性表达式。

于是，我们求M的m次方，可以先求出M的m次方的E、M^1、M^2、M^3......M^p-1线性表达式。根据快速幂运算的理论，我们仅仅需要求出M的1次方，M的2次方，M的4次方等等，指数全部为2的幂。于是，我们有了如下求M矩阵的快速幂的步骤：

（1）、向量（m0 , m1 , m2 , m3......mp-1）分别表示M的下标次方的系数。初始化的时候m1 = 1。这正好表示M的1次方。

（2）、向量乘以自身一次，其复杂度为O(p^2)，则可以求得M的平方，这是，向量的长度扩大，我们得到了M的平方，E、M^1、M^2、M^3......M^2p-2的系数，由于上面展式的作用，我们可以将M^p、M^p+1、M^p+2、......M^2p-2依次降低其指数的大小，注意，顺序需要从2p-2次幂开始降低，这样的算法复杂度为O(p^2)。

（3）、重复步骤（2）可以依次获取我们需要的矩阵。而这个过程时间复杂度为O(p^2)。

到此，所有的优化全部解决，我自己的代码经过这几个大的优化，已经可以通过51nod上的”字符串数量V3“，虽然优化的余地还有很多，但达到目的即可。

断断续续写了快一天了，越写越觉得自己对问题的表述不是很好，希望看到这里的人，能够看懂我粗糙的表述，并干掉这道题目！下面是我的代码的一小部分：

大素数：

/** * 大素数,求模 */public static long MOD = 1000000007;

计算排列组合的代码：

/** * 计算排列组合(取模) *  * @param n * @param m * @return */public static long C(long n, long m) {long i, a, b, p;if (n < m) {return 0; // n不能小于m}p = 1;i = n - m;a = i < m ? i : m;b = i > m ? i : m;long middle = 0;for (i = 1; i <= a; i++) {middle = p * b;middle %= MOD;middle *= extendGCD(i, MOD);middle %= MOD;p += middle;p %= MOD;}return p;}

计算(2x-1)^n的系数(幂次从0到n)：

/** * 计算(2x-1)^n的系数(幂次从0到n) *  * @param n */public static long[] getFactor(int n) {long fac[] = new long[n + 1];long k = 1; // 2的幂(开始为2的0次幂)long count = n; // -1的幂次int sign = 0;for (int i = 0; i <= n; ++i) {sign = count-- % 2 == 0 ? 1 : -1; // -1提供的系数fac[i] = (C(n , i) % MOD);// 组合数提供的系数fac[i] *= k; // 计算2的幂次fac[i] %= MOD;if(sign == -1){fac[i] = MOD - fac[i];}k *= 2;k %= MOD;}return fac;}

获取1到x的k次方幂的系数：

/** * 获取1到x的k次方幂的系数(结果数组一次代表: 幂次从1到(n + 1)) * @param k * @return */public static long[] getFactorOfKPow(int k){long fac[] = new long[k + 1];fac[0] = extendGCD(k + 1 , MOD);long mid = 0 , val = 0;for(int i = 1 ; i <= k ; ++ i){fac[i] = extendGCD(k - i + 1, MOD);mid = C(k , i);for(int j = 0 ; j < i ; ++ j){val = C(k - j + 1 , i - j + 1);val *= (MOD - fac[j]);val %= MOD;mid += val;mid %= MOD;}fac[i] *= mid;fac[i] = (fac[i] + MOD) % MOD;}//系数倒转,编程1次幂,2次幂...k次幂k ++;long res[] = new long[k];for(int i = 0 ; i < k ; ++ i){res[k - 1 - i] = fac[i];}return res;}

获取1到(2 * x - 1)的k次幂和系数：

/** * 获取1到(2 * x - 1)的k次幂和系数: 结果数组一次代表幂次从0到k+1 * @param k * @return */public static long[] getFactorOfPolynomial(int k){long res[] = new long[k + 2];for(int i = 0 ; i < res.length ; ++i){ // 初始化置零res[i] = 0;}long powK[] = getFactorOfKPow(k);// 获取1到(k + 1)的系数long powF[] = null;long mid = 0;for(int i = 0 ; i < powK.length ; ++ i){powF = getFactor(i + 1);for(int j = 0 ; j < powF.length ; ++ j){mid = powF[j];mid *= powK[i];mid %= MOD;res[j] += mid;res[j] %= MOD;}}return res;}

到此为止了，通过这些代码，和之前求逆元的代码，可以获得求取1到2x - 1的幂和多项式的系数，有了这个，计算【合法链】就不是很困难了。代码中所有的都经过了模10^9 + 7。

最后，附上几个结果，以便能够进行测试：

n = 8 , m = 8 结果为1530445;

n = 1000000000 , m = 1000000000 结果为917875839;