全部
博文收录
Ada助手
问答

用12的骨牌覆盖mn的矩形的问题

绿色夹克衫 2012-08-27 08:32:04

最近想写个程序算算，用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

...全文

给本帖投票

1491 40 打赏收藏转发到动态举报

写回复

用AI写文章

单链表法编写malloc运行时出现段错误，请教各位高手，帮忙查出错误！！！>>

40 条回复

切换为时间正序

请发表友善的回复…

发表回复

FancyMouse 2012-09-04

打赏
举报

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

mathe 2012-09-03

打赏
举报

上面的程序有个缺点，需要先计算出结果的平方然后开方，由于大数计算数据越大计算量越大，所以比较浪费。现在做了一个更新版本，可以直接计算结果而不需要开方:
(其中，假设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)
}

jediael_lu 2012-09-01

打赏
举报

数学，又是数学。

绿色夹克衫 2012-08-31

打赏
举报

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

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

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

绿色夹克衫 2012-08-31

打赏
举报

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

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

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

cnmhx 2012-08-31

打赏
举报

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

绿色夹克衫 2012-08-31

打赏
举报

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

[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 2012-08-31

打赏
举报

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

mathe 2012-08-31

打赏
举报

根据上面算法写了个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 2012-08-31

打赏
举报

其实变换以后的公式从复杂度上可能还高一些，好处是不需要浮点运算。
原公式是计算所有关于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 2012-08-31

打赏
举报

[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]
没问题

绿色夹克衫 2012-08-31

打赏
举报

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]

绿色夹克衫 2012-08-30