用1*2的骨牌覆盖m*n的矩形的问题

[绿色夹克衫]

最近想写个程序算算，用1*2的骨牌覆盖m*n的矩形有多少种不同的方式，Wiki上的公式也很漂亮。
http://en.wikipedia.org/wiki/Domino_tiling#CITEREFKasteleyn1961

但这种计算怎么保证精度呢？还是可以在计算前先把某些项合并？

这是我现在写的一个C#程序，用的double，算稍微大一些的m*n时，就不准了。
http://www.51nod.com/question/index.html#!questionId=316

[mathe]

发现那个公式可以利用第二类切皮雪夫多项式变成一个n阶高斯整数矩列式问题，可以定点计算

[mathe]

高斯整数矩阵行列式，上面漏俩字

[mathe]

动态规划，从m*(k-1)行到m*k行，记录填满m*(k-1)行时有多少个第k行的格子被覆盖作为状态即可。
可以参考一个更加复杂题目的解答:
http://tieba.baidu.com/p/395424512?pid=4076687686&cid=0#4076687686

[mathe]

可以参考链接围棋问题的82楼

[绿色夹克衫]

百度贴吧似乎有问题，这个帖子看不了呀！

[Quote=引用 18 楼 的回复:]

可以参考链接围棋问题的82楼
[/Quote]

[mathe]

好像是出问题了，将出问题后IE地址栏中前面部分都删除，仅留url=后面的东西即可

[绿色夹克衫]

嗯，看到了，学习一下。

[Quote=引用 20 楼 的回复:]
好像是出问题了，将出问题后IE地址栏中前面部分都删除，仅留url=后面的东西即可
[/Quote]

[绿色夹克衫]

wayne 用 Mathematica 计算了100*200的情况，是个10053位数，比较好奇Mathematica是怎么解决这个问题的。

http://bbs.emath.ac.cn/thread-4579-2-1.html

[绿色夹克衫]

mathe，我能够把你这个答案转走么？会标明是你的原创的

[Quote=引用 30 楼 的回复:]

其实变换以后的公式从复杂度上可能还高一些，好处是不需要浮点运算。
原公式是计算所有关于j,k时
4(cos^2((j*pi)/(m+1))+cos^2((k*pi)/(n+1)))
的乘积的1/4次方
上面可以写成||2(i*cos((j*pi)/(m+1))-cos((k*pi)/(n+1)))||^2
我们查看第二类切皮雪夫多项式
U_n(x)=2^n \prod (x - co……
[/Quote]

[mathe]

[Quote=引用 36 楼 的回复:]

mathe，我能够把你这个答案转走么？会标明是你的原创的

引用 30 楼 的回复:

其实变换以后的公式从复杂度上可能还高一些，好处是不需要浮点运算。
原公式是计算所有关于j,k时
4(cos^2((j*pi)/(m+1))+cos^2((k*pi)/(n+1)))
的乘积的1/4次方
上面可以写成||2(i*cos((j*pi)/(m+1))-cos((k*pi)/(n+1……
[/Quote]
没问题

[mathe]

其实变换以后的公式从复杂度上可能还高一些，好处是不需要浮点运算。
原公式是计算所有关于j,k时
4(cos^2((j*pi)/(m+1))+cos^2((k*pi)/(n+1)))
的乘积的1/4次方
上面可以写成||2(i*cos((j*pi)/(m+1))-cos((k*pi)/(n+1)))||^2
我们查看第二类切皮雪夫多项式
U_n(x)=2^n \prod (x - cos((k*pi)/(n+1))
所以结果可以写成
sqrt(||\Prod_{j=1}^m U_n( i*cos((j*pi)/(m+1)))||）
我们记n阶矩阵B_n为除了主对角线边上元素为1，其余为0的矩阵，也就是Bn_ij在|i-j|=1时为1，其余时候为0
那么根据链接第二类切皮雪夫多项式又可以写成det(2x*I_n+B_n),其中I_n为n阶单位阵
所以结果可以写成
sqrt(||det(\Prod_{j=1}^m (2i*cos((j*pi)/(m+1))I_n+B_n))||)
=sqrt(||det(2^m*\Prod_{j=1}^m(i*B_n/2-cos((j*pi)/(m+1))I_n))||)
=sqrt(||det(U_m(i*B_n/2))||)
所以只要我们先计算整系数多项式U_m(x/2),然后将i*B_n代入得出n阶矩阵，然后求行列式的绝对值再开平方即可。

[mathe]

根据上面算法写了个Pari/Gp程序

MB(n)=
{
   local(B);
   B=matrix(n,n);
   for(u=1,n-1,
      B[u,u+1]=1;
      B[u+1,u]=1
   );
   B
}
MI(n)=
{
    local(B);
    B=matrix(n,n);
    for(u=1,n, B[u,u]=1);
    B
}

V(m,n)=
{
    local(M1,M2,M3,B);
    B=MB(n)*I; 
    M1=MI(n);
    M2=B;
    for(u=2,m,
      M3=B*M2-M1;
      M1=M2;M2=M3
    );
    sqrtint(abs(matdet(M2)))
}

计算V(100,100)还是比较容易的

[mathe]

不过计算V(200,100)我的Pari/Gp受不了了,Pari/Gp还是肉了一些。
另外对于偶数m,最后那个矩阵必然是整数矩阵

[绿色夹克衫]

感觉这个方法我花些时间学习一下，应该能够看懂

[Quote=引用 30 楼 的回复:]

其实变换以后的公式从复杂度上可能还高一些，好处是不需要浮点运算。
原公式是计算所有关于j,k时
4(cos^2((j*pi)/(m+1))+cos^2((k*pi)/(n+1)))
的乘积的1/4次方
上面可以写成||2(i*cos((j*pi)/(m+1))-cos((k*pi)/(n+1)))||^2
我们查看第二类切皮雪夫多项式
U_n(x)=2^n \prod (x - co……
[/Quote]

[cnmhx]

啥叫带方向的截断？
就是四舍五入roundoff，只有舍或者入（对于你的问题，只有舍去）！

[绿色夹克衫]

mathe大牛的这个方法，我表示暂时也看不懂。

[Quote=引用 27 楼 的回复:]

发现那个公式可以利用第二类切皮雪夫多项式变成一个n阶高斯整数矩列式问题，可以定点计算
[/Quote]

[绿色夹克衫]

原来是这个意思，明白了。

[Quote=引用 34 楼 的回复:]

啥叫带方向的截断？
就是四舍五入roundoff，只有舍或者入（对于你的问题，只有舍去）！
[/Quote]

[jediael_lu]

数学，又是数学。

[mathe]

上面的程序有个缺点，需要先计算出结果的平方然后开方，由于大数计算数据越大计算量越大，所以比较浪费。现在做了一个更新版本，可以直接计算结果而不需要开方:
(其中，假设m是偶数而且m=2s)
V2(s,n)=
{
local(M1,M2,M3,S,B);
B=MB(n)*I;
M1=B;
S=MI(n)+M1;
M2=B*B-2*MI(n);
for(u=2,s,
S=S+M2;
M3=M2*B-M1;
M1=M2;M2=M3
);
matdet(S)
}

[FancyMouse]

直接用数值方法计算是3次方，行列式目测要5次方……
开方一次倒是小问题，不增加复杂度的。
不过Chebyshev polynomial of 2nd kind这方法还真不错。能把答案弄出个漂亮的形式。