51nod 1034 骨牌覆盖 V3 | URAL 1594 Aztec Treasure

题目：
http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1034
http://acm.timus.ru/problem.aspx?space=1&num=1594

题意：
给定一个 n×m 的棋盘，用 1×2 的骨牌去覆盖，问有多少种放法能够用骨牌盖满所有位置。
骨牌覆盖 V3 ：求方案数精确值
Aztec Treasure ：求方案数 mod(109+7)
1≤n,m≤100

题解：
有解的情况实际上只有 3775 种， URAL 上的题目可以打表，或者写完之后猜一下二次剩余的符号。

这里提供一个算具体方案数（不取模）但是不用高精度有理数的做法，会用到高精度整数的只有CRT和开方的部分，其他部分可以用 __int128_t 完成，如果只是需要打表，可以用 C++ 完成主要部分， python 完成大数部分。

这个平面图比较特殊，所以FKT Algorithm计算完美匹配数平方的那个矩阵可以不用写成 nm×nm 的，而是用一个关于 n×n 的矩阵的 m 次多项式表示。

为了避开我不会的内容，我们从Wikipedia上给出的公式开始介绍。

令 rn(i)=2cos(πin+1) ，则公式可以表示为

∏j=1n∏k=1m(r2n(j)−i2r2m(k))−−−−−−−−−−−−−−−−−−⎷4,

其中 i=−1−−−√ 。

注意到 n,m 之间没有什么限制，我们可以将公式表示成 m 个多项式的乘积，即

∏k=1mF(irm(k))F(−irm(k))−−−−−−−−−−−−−−−−−−−−√4,

其中

F(x)=∏j=1n(x−rn(j)).

不难发现 −rm(k)=rm(m+1−k) ，所以公式其实是

∣∣∣∏k=1mF(irm(k))∣∣∣−−−−−−−−−−−−⎷,

注意这里有绝对值符号。

对于这个 F(x) ，我们可以在Wikipedia找到它和第二类切比雪夫多项式 Un(x) 的关系，可以根据 Un(x) 零点恰好是 rn(j)2 (j=1,2,⋯,n) 以及 Un(x) 的最高次项系数为 2n 得知 F(x)=Un(x2) ，和切比雪夫多项式拥有相同的特征——它们都是整系数多项式。

然而上面那个乘积还是并不好算的，为此我们需要更棒的性质，在MathWorld可以找到它的矩阵特征是

Un(x)=det(2xIn+An),

其中 In 是 n 阶单位矩阵， An 是只在主对角线两侧的对角线上值为 1 的 n 阶矩阵，也即 An(i,j)=[|i−j|=1] 。

然后你带进去就可以发现

∏k=1mF(irm(k))=∏k=1mdet(irm(k)In+An)=(−i)nmdet(∏k=1m(iAn−rm(k)In)),

机智的你一定发现这里出现了这样一个多项式

G(x)=∏k=1m(x−rm(k))=Um(x2),

所以公式变成了

∣∣(−i)nmdet(G(iAn))∣∣−−−−−−−−−−−−−−−−−√=|det(G(iAn))|−−−−−−−−−−−√.

根据Wikipedia，有

U0(x)=1,U1(x)=2x,Un(x)=2xUn−1(x)−Un−2(x) (n>1),

那么就有

Um(iAn2)=iAnUm−1(iAn2)−Um−2(iAn2),

仔细观察可以发现 Um(iAn2) 的每个元素要么实部为 0 要么虚部为 0 ，所以我们可以维护 imUm(iAn2) 从而避免复数运算，而一个矩阵乘以 An 的复杂度实际上是 O(n2) 的，所以我们可以 O(n2m) 递推求出 Um(iAn2) 。

剩余的内容是对一个整系数矩阵求行列式并开方的过程，我采取的方法是取多个大质数计算模意义下的行列式，例如取 T 个质数的时候复杂度是 O(Tn3)。由于行列式的值一定是一个平方数，可以利用二次剩余的技巧确定行列式的正负号（对于 4k+3 型素数， x 和 −x 在模意义下勒让德符号不同）。

开方和求矩阵的部分都不是最慢的部分，主要是求行列式比较慢，用 C++ 实现较好。

代码就不给了，比较短，而且内容已经在上面描述完毕了。

2017.09.06 更新：注意到矩阵具有对称性，最终消元得到的 trace 将满足奇数位置的乘积等于偶数位置的乘积，利用该性质可以避免开根号。