Formula for tiling a chessboard with dominos

[asmeurer]

See http://www.johndcook.com/blog/2016/08/19/how-many-ways-can-you-tile-a-chessboard-with-dominoes/. The function from the post can be adapted for SymPy:

def num_tilings(m, n):
    m, n = sympify((m, n))
    prod = 1
    for k in range(1, m+1):
        for l in range(1, n+1):
            prod *= 2*cos(pi*k/(m+1)) + 2*I*cos(pi*l/(n+1))
    return sqrt(prod)

This is a nice testbed for simplify, especially trigsimp. The 2 x 2 case gives a wrong answer!

In [78]: num_tilings(2, 2)
Out[78]:
  ________   ________   _______   _______
╲╱ -1 - ⅈ ⋅╲╱ -1 + ⅈ ⋅╲╱ 1 - ⅈ ⋅╲╱ 1 + ⅈ

In [79]: simplify(num_tilings(2, 2))
Out[79]: -2

In [80]: num_tilings(2, 2).evalf()
Out[80]: 2.0 + 1.35525271560688e-20⋅ⅈ

The 4 x 4 case takes a lot of work to be simplified. Just simplify doesn't do it. I hit a bug with sqrtdenest:

In [93]: sqrtdenest(simplify(num_tilings(4, 4)))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-93-48450a8958ed> in <module>()
----> 1 sqrtdenest(simplify(num_tilings(4, 4)))

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in sqrtdenest(expr, max_iter)
    130     expr = expand_mul(sympify(expr))
    131     for i in range(max_iter):
--> 132         z = _sqrtdenest0(expr)
    133         if expr == z:
    134             return expr

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in _sqrtdenest0(expr)
    240         args = expr.args
    241         if args:
--> 242             return expr.func(*[_sqrtdenest0(a) for a in args])
    243     return expr
    244

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in <listcomp>(.0)
    240         args = expr.args
    241         if args:
--> 242             return expr.func(*[_sqrtdenest0(a) for a in args])
    243     return expr
    244

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in _sqrtdenest0(expr)
    240         args = expr.args
    241         if args:
--> 242             return expr.func(*[_sqrtdenest0(a) for a in args])
    243     return expr
    244

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in <listcomp>(.0)
    240         args = expr.args
    241         if args:
--> 242             return expr.func(*[_sqrtdenest0(a) for a in args])
    243     return expr
    244

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in _sqrtdenest0(expr)
    229                 if len(args) > 2 and all((x**2).is_Integer for x in args):
    230                     try:
--> 231                         return _sqrtdenest_rec(n)
    232                     except SqrtdenestStopIteration:
    233                         pass

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/simplify/sqrtdenest.py in _sqrtdenest_rec(expr)
    270     if not expr.is_Pow:
    271         return sqrtdenest(expr)
--> 272     if expr.base < 0:
    273         return sqrt(-1)*_sqrtdenest_rec(sqrt(-expr.base))
    274     g, a, b = split_surds(expr.base)

/Users/aaronmeurer/Documents/Python/sympy/sympy/sympy/core/relational.py in __nonzero__(self)
    193
    194     def __nonzero__(self):
--> 195         raise TypeError("cannot determine truth value of Relational")
    196
    197     __bool__ = __nonzero__

TypeError: cannot determine truth value of Relational

To actually simplify it, the only way I found was to kill the outer square root and simplify first.

In [98]: sqrt(simplify(num_tilings(4, 4)**2))
Out[98]: 36

I tried expand_complex but that just produced a huge expression with a bunch of cosines of arctangents.

The 8 x 8 version seems to be too hard. It doesn't know how to reduce the cosines (like cos(pi/9)). Of course, you can get the answer using evalf, but I can't get it to symbolically simplify.

[asmeurer]

The 2 x 3 case also gives a wrong answer. It must be wrong somewhere in auto-simplification:

In [103]: num_tilings(2, 3)
Out[103]:
    ___________   ___________   __________   __________
ⅈ⋅╲╱ -1 - √2⋅ⅈ ⋅╲╱ -1 + √2⋅ⅈ ⋅╲╱ 1 - √2⋅ⅈ ⋅╲╱ 1 + √2⋅ⅈ

In [104]: simplify(num_tilings(2, 3))
Out[104]: 3⋅ⅈ

In [105]: num_tilings(2, 3).evalf()
Out[105]: 1.35525271560688e-20 + 3.0⋅ⅈ

[jksuom]

It seems that the formula is not quite correct. The product can be negative. The absolute value should be taken: sqrt(abs(prod)).

[asmeurer]

Supposedly it doesn't. I'll need to work out an example by hand to check. Regardless, something is wrong for the 2 x 2 case where simplify and evalf disagree.

[asmeurer]

I modified the formula to not take a square root and

In [26]: expand(num_tilings(2, 2))
Out[26]: 4

In [27]: expand(num_tilings(2, 3))
Out[27]: -9

Unless there is some other bug, it seems you are right for the 2 x 3 case. simplify definitely has a bug for the 2 x 2 case, though.

[ylemkimon]

Latest development version seems to give correct answer(2).

[asmeurer]

Indeed, it no longer incorrectly splits the square root

>>> num_tilings(2, 2)
sqrt((-1 - I)*(-1 + I)*(1 - I)*(1 + I))
>>> simplify(num_tilings(2, 2))
2

[asmeurer]

And even simplify(num_tilings(4, 4)) works now.