回文和数

medie2005

回文和数是指：可由(a*(a+1))/2（a为自然数）表示，并且形式回文的数。 例如： 109*110/2=5995 是回文和数。 2662*2663/2=3544453 是回文和数。 50281*50282/2=1264114621 是回文和数。 能否找出10^40内的所有回文和数？

无心人

你知道袁峰的工作么？ 所以你的要求比较高啊

mathe

10^40以内要求并不高。

medie2005

是的，这个题目现在的世界记录是由袁峰保持，目前最大的回文和数有46位： 8215588527084263951180440811593624807258855128 它可写为 （128184152897963669861552*（128184152897963669861552+1））/2.

mathe

现在有没有人有兴趣解决这个题目？这个同平方回文数的代码几乎可以完全相同。 当然算法时间复杂度很难改善，也就是计算n位需要$O(10^{n/4})$,但是我觉得好的算法可以将大部分乘法淘汰，从而达到较快的速度，应该可以超过袁的速度。下面采用平方回文数作为例子，说明可以采取的比较快的算法： 假设我们现在要计算n位平方回文数，我们可以对大部分数据段采取10000进制进行计算。为此，我们首先需要构造一个4位数的逆序数表rev[x],如rev[1]=1000,rev[1234]=4321等等。 然后，假设n位回文数X，它的平方根为Y,那么Y的位数m为$[{n+1}/2]$. 我们首先只对Y的最前面$4*[m/8]$和$4*[m/8]$位进行搜索，这样中间还会余下$m-8[m*8]$位（最少0位，最多7位）最后特殊处理。 首先我们从Y的最低4位开始，我们可以构造平方表,伪代码如下，需要注意我们要采用10000进制计算而不是普通2进制计算。

1. s=0;
2. step=1;
3. for(i=0;i<10000;i++){
4. sqr_table[i]=s;
5. s += step;
6. step+=2;
7. }

复制代码

这样做的好处是不需要任何乘法运算。 然后我们对于sqr_table[ i ]中每项（10000进制两位数）的低位，可以构造一个平方逆表 如sqr_table[1]=0*10000+1,我们可以得到inv_sqr[1]=1 (1^2最低位是1） 需要注意这个逆表格不是单射，也就是可以有多个数字平方低位是相同的。 让后我们需要对Y的最高四个十进制位进行处理，这个处理方法根据n是奇数还是偶数有所不同，而通常来说，Y的最高四位可以决定X的最高四位数。我们看看这个关系如何，假设Y的最高4位为u,X的最高4位为v 那么我们知道 $u*10^{m-4}<=Y<(u+1)*10^{m-4}$, $v*10^{n-4}<=X<(v+1)*10^{n-4}$ 所以我们可以得到关系式$u^2*10^{2m-8}<(v+1)*10^{n-4}, v*10^(n-4)<(u+1)^2*10^{2m-8}$ 其中如果n是偶数n=2m,如果n是奇数n=2m-1,我们需要分别处理两种情况。 如果n=2m,同样我们将$u^2$和$(u+1)^2$根据其平方表高位(10000进制),可以直接得出v的范围，也就是 sqr_table[ u ].hi=v 由此，我们可以知道，对于指定的Y的最高4为u可以得到一个可以使用的X的最高四位数的范围，对于这个范围内每个数v,我们通过rev表格查询到v的倒序数，然后根据inv_sqr表格查询处Y最低4位的可能情况。 需要注意的是对于第一层计算，会出现很多数据的inv_sqr表格查询结果不唯一，但是对于另外一些数据，查询结果不存在。

mathe

同样对于n是奇数的情况，我们也有类似的结果，只是这时候我们需要的是$u^2/10$和$(u+1)^2/10$,这里牵涉到除10运算。由于第一层计算所占比例不高，这里就不做进一步优化了。 下面开始第二层数据计算，也就是对于每种Y的最低4位和最高4位确定的情况，（这时X的最高4位和最低4位也确定了） 这时假设最低4位为x1,最高4位为u都已经知道，同样对应X的最低4位和最高4位已经确定 我们现在继续假设Y的次低4位的值位x2,次高4位为u2. 同样我们采用类似上面的计算方法去计算$(x2*10000+x1)^2$,伪码如下

1. s=sqr_table[x1];
2. step=(x1+x1)*10000;///乘10000相当于移位运算
3. for(x2=0;x2<10000;x2++){
4. sqr_table_level_2[x2]=s;
5. s+=step;
6. step+=2*10000*10000;
7. }

复制代码

而且上面加法过程中我们可以知道的是sqr_table_level_2的最低10000进制位都相同，所以计算过程中这部分实际上可以不参与，所以上面是一个3位10000进制数的累加过程。

mathe

然后对于高位，我们需要计算的是 $(u*10^4+u2)^2*10^{2m-16}<=Y^2=X<(v*10^4+v2+1)*10^{n-8}$ $<(v*10^4+v2)*10^{n-8}=X=Y^2<(u*10^4+u2+1)^2*10^{2m-16}$ 这里我们发现对于第二层，我们不能像第一层那样只产生一个平方表了，我们需要两个了 所以通过类似方法，我们可以产生平方表 $(u*10^4+u2)^2$,其中u是常数，u2是变量。同样这个过程中不需要任何乘法运算。只需要O(10000)次3位10000进制加法。（我们还需要判断第三位计算是否产生进位，如果进位说明越界了。 同样，有了这个关系以后，我们就可以得到u2到v2的映射表格。这时这个表格是几乎一一映射？（应该除了少数边缘上的数据，都是一一映射）。 现在我们查看sqr_table_level_2的倒数第二位（10000进制位），同样我们可以建立这个的逆表格（这个表格是单射？？？） 然后对于每个u2，我们可以查出对应的v2,然后通过rev表格算出X对应位置的4位，然后通过第二层逆平方表格查出Y的次低4位x2.

mathe

然后，在现在的数据范围之内，有可能我们还需要进行第三层计算。对于每个计算出来的(u,u2),(x,x2),我们需要重复一个类似上面的过程，得出数据的第三层，于是我们有候选数据(u,u2,u3),(x,x2,x3).这时我们已经确定了Y的24位数据（如果平方数总位数低于48位，那么停止在第二层计算，已经确定的位数为16位) 然后我们查看余下的数据位数。如果余下数据位数为0或1，我们跳过这一步。 如果余下数据位数大于2位，（如果奇数位，中间一位总是跳过这一步） 我们同样可以将余下数据位数对称分成两个部分，同样进行类似前面的计算，不过这时，多出来的数据部分采用的进制将低于10000（比如是1000，100或10），同样我们可以事先制作好一个位数较短的逆序表格用于辅助计算。 通过上面的方法，可以穷举Y中最多除了一位（最中间一位）以外的所有数据。

mathe

接下去，我们需要对候选数据Y进行检验，是否$Y^2$是一个真正的平方数。 我们现在已经拥有的数据是 $U^2=(u*10^{8+h}+u2*10^{4+h}+u3*10^h+u4)^2$ (假设u4是h位) 和 $Z^2=(x+x2*10^4+x3*10^8+x4*10^12)^2$ 我们需要计算的是$(U*10^m+Z)^2$或$(U*10^(m+1)+Z)^2$(对于后一种情况，我们最后还需要根据最中间那位数的值，加9次去枚举所有情况。这步需要我们还要进行一次额外大数乘法$U*Z$,然后还可能需要乘上一个$10^h$. 总之，这部分计算相对前面计算由于存在乘法会比较费时。 所以我们应该想办法去尽量减少这一部分运算。

mathe

一种可以试验的方法是在h比较大的情况，比如h=6或h=7, 我们可以继续前面的累加过程（需要注意的是这个时候，前后的数据将有一部分重叠）多计算一层数据，然后再判断这个算出来的数据是否合法（合法性检查同样可以通过一个事先构着的表格来判断）。 但是对于h<=5这样计算应该对速度不会有好处。 另外还有一个可以使用的信息是在n是偶数时，Y必须是11的倍数，这个在最后一层搜索过程中可以减少计算量。