Phyllotaxis and Fibonacci

Dobble Connect

 

 

I was given a new version of Dobble for Xmas called Dobble Connect.  Dobble consists of stack of cards with pictures on it where each pair of cards has one match, and each pair of symbols appears on one card.  This is a new version which has 91 symbols in total, 10 symbols per card, and 91 cards.  The fingerprint of the original Dobble was (57 symbols, 8 symbols, 57 cards).

An additional innovation in Dobble Connect is that the cards are hexagonal and grouped into a colour per player.  Players build a tiling on the tabletop, placing cards next to each other as soon as they spot a matching symbol between their top card and one of the cards already placed.  The first player to get four of their own colour in a row wins.

How is it possible for each pair of cards to have one and only one symbol in common, and for each pair of symbols to appear on one and only one card?  Does that have something to do with the strange number 91?  The answer is "yes" and it turns out this is only possible for packs of $p^{2n} + p^n + 1$ cards, where $p$ is prime, and $n \ge 1$.  Let's see why.

Oddly, it's easiest to begin with an infinite version of the game.  Cards have an infinite number of symbols, but each pair has only one in common.  How do we construct such a pack of cards?  We shall start by making each card a two-dimensional subspace of $\mathbb{R}^3$

 


Each pair of cards then intersects in a unique one-dimensional subspace of $\mathbb{R}^3$.  So, if we choose to interpret 1D subspaces as "symbols" then each card contains an infinite number of symbols, but each pair of cards have exactly one in common.  Furthermore, any pair of symbols is present in one and only one card, since the basis vectors of each 1D subspace combine to form a basis of a 2D subspace.  But how do we go from here to Dobble?

The argument above used $\mathbb{R}^3$ which is a 3 dimensional vector space over the infinite field $\mathbb{R}$.  However, none of the arguments used above were field dependent.  (That the intersection of two 2D subspaces of a 3D vector space is a 1D subspace, and that each N-dimensional subspace can be defined by N linearly independent vectors.)  Let's see what happens if we use a finite field.

For each prime $p$ the integers modulo $p$ form a field, namely $\mathbb{F}_p = \{0, 1, ..., p-1\}$.  It is easy to see that addition (modp) and multiplication (modp) meet most of the axioms required of a field.  The only non-obvious one is the existence of a multiplicative inverse.  This is where the fact that $p$ is prime comes in.  To find the inverse of $a$ we must use Euclid's algorithm to find the greatest common divisor $r$ of $a$ and $p$, as well as integers $x$ and $y$ satisfying $xa + yp = r$.  We know $r = 1$ because $p$ is prime, which means $x$ must be the inverse of $a$ modulo $p$.

So, what happens if we choose $\mathbb{F}_p$ as our field?  The number of symbols is the number of one-dimensional subspaces of $\mathbb{F}_p^3$.  There are $p^3 - 1$ options for the basis vectors (we can't use zero) but we want to avoid double counting caused by the fact that $v$ and $\lambda v$ span the same 1D subspace.  So let's only consider vectors whose leading non-zero component is $1$.  This gives us $p^2 + p + 1$ basis vectors all of which span unique 1D subspaces:

$$
\begin{align}
(1,x,y) &\longleftarrow p^2\ \text{vectors} \\
(0,1,x) &\longleftarrow p\  \text{vectors} \\
(0,0,1) &\longleftarrow 1\ \text{vector} \\
\end{align}
$$

What about the number of cards?  This time we need to count pairs of basis vectors, and to prevent double counting we must assume that

  1. both basis vectors have $1$ in their leading non-zero component
  2. the second basis vector has more leading zeros than the first
  3. each basis vector has a $0$ where it's partner has it's leading non-zero value

This tells us that the number of cards is also $p^2 + p + 1$ since

$$
\begin{align}
\{(1,0,x),(0,1,y)\} &\longleftarrow p^2\ \text{pairs} \\
\{(1,x,0),(0,0,1)\} &\longleftarrow p\ \text{pairs} \\
\{(0,1,0),(0,0,1)\} &\longleftarrow 1\ \text{pair} \\
\end{align}
$$

And the number of symbols per card is $p + 1$ since if $u,v$ are the basis vectors for a 2D subspace the 1D subspaces of this are defined by the basis vectors

$$
\begin{align}
u + \lambda v &\longleftarrow p\ \text{vectors} \\
v &\longleftarrow 1\ \text{vector} \\
\end{align}
$$

If we set $p=7$ we get 57 cards, 57 symbols, and 8 symbols per card, just like in the original Dobble.  However, this doesn't work for Dobble Connect which has 91 cards, 91 symbols, and 10 symbols per card.  For that we need a field with 9 elements.  Luckily there are fields of size $p^n$ for every non-negative $n$ (in fact these are the only permitted sizes for finite fields.)  So all that remains is to construct $\mathbb{F}_{3^2}$.

First we need to construct the polynomials over $\mathbb{F}_p$, namely $\mathbb{F}_3[X]$.  These are polynomials like $a(X) = X^3 + 2X + 1$ where all of the coefficients are from $\mathbb{F}_3$.  These add and multiply just like ordinary polynomials.  Next we need to choose an irreducible polynomial $p(X)$ of degree $n$.  This is a polynomial which is not divisible by any polynomial of degree $m$ with $0 < m < n$.  (It's a simple exercise in combinatorics to show that there will always be at least one.)  In the case $p=3, n=2$  we can choose $p(X) = X^2 + 1$.  The final step is to construct $\mathbb{F}_{p^n}$ by doing arithmetic modulo $p(X)$ in $\mathbb{F}_p[X]$.  How do we know that $a(X)$ has an inverse in $\mathbb{F}_{p^n}$?  The answer is that Euclid's algorithm works equally well with polynomials as with integers, so the argument we used to show that $\mathbb{F}_p$ has inverses can also be used to show $\mathbb{F}_{p^n}$ has inverses.

We now have a mechanism for constructing a Dobble pack for any prime $p$ and any integer $n \ge 1$, with $p^{2n} + p^n + 1$ cards, $p^{2n} + p^n + 1$ symbols, and $p^n + 1$ symbols per card.  And the following 500 lines of code does just that.

#!/usr/bin/env python3
"""
Functions for creating dobble packs - no error or bounds checking yet
"""

from functools import reduce


def gcd(a, b) -> tuple:
    """
    find gcd n of a and b (where a > b)
    return (n, x, y) where x and y are the coefficients in the equation

        n = x*a + y*b

    NB: should work with either ints or Polynomials (over any field)
    """
    q = a // b
    r = a - q * b

    # we need to do this to support case where a and b are polynomials
    zero = a * 0
    one = zero + 1

    if r == zero:
        return b, zero, one
    n, x, y = gcd(b, r)

    return n, y, x - y * q


class FiniteField_p:
    """
    Finite field of p elements (where p is prime)
    Usage:

          field = FiniteField_p(7)
          element1 = field(3)
          element2 = field(4)
          print(element1 * element2) # prints 5
    """

    class Element:
        def __init__(self, k: int, field):
            self._field = field
            self._p = field._p
            self._k = k % self._p

        def field(self):
            return self._field

        def inv(self):
            return type(self)(gcd(self._p, self._k)[2], self._field)

        def __truediv__(self, other):
            if type(other) is not type(self):
                return self / type(self)(other, self._field)
            return type(self).__mul__(self, other.inv())

        def __mul__(self, other):
            if type(other) is type(self):
                return type(self)(self._k * other._k, self._field)
            if type(other) is Polynomial:
                return Polynomial([self]) * other
            if type(other) is Vector:
                return Vector([self * other[i] for i in range(other.dim)])
            return type(self)(self._k * other, self._field)

        def __rmul__(self, other):
            return self * other

        def __add__(self, other):
            if type(other) is type(self):
                return type(self)(self._k + other._k, self._field)
            return type(self)(self._k + other, self._field)

        def __radd__(self, other):
            return self + other

        def __neg__(self):
            return type(self)(-self._k, self._field)

        def __sub__(self, other):
            return type(self)(self._k - other._k, self._field)

        def __eq__(self, other):
            if type(other) is type(self):
                return self._k == other._k
            return self._k == other

        def __ne__(self, other):
            return not self == other

        def __repr__(self):
            return f"{self._k}"

    def __init__(self, p):
        self._p = p

    def __call__(self, k: int):
        if type(k) == self.Element:
            return k
        return self.Element(k, self)

    def get_elements(self):
        return [self.Element(i, self) for i in range(self._p)]


class FiniteField_p_n:
    """
    Finite field of p^n elements (where p is prime, n > 1)
    Usage:

          subfield = FiniteField_p(3)
          poly = Polynomial([subfield(1), subfield(0), subfield(1)]) # irreducible 1 + X^2
          field = FiniteField_p_n(subfield, 2, poly)                 # field of 3^2 elements
          element1 = field(4)                                        # 1 + X
          element2 = field(5)                                        # 2 + X
          print(element1 * element2) # prints 1                      # 1 = (1 + X)(2 + X) - (1 + X^2)
    """

    class Element:
        def __init__(self, k: int, field):
            self._field = field
            if type(k) is Polynomial:
                _k = 0
                power = 1
                for v in k._a:
                    _k += v._k * power
                    power *= self._field._p
                self._k = _k
                self._poly = k
            else:
                self._k = k % field._N

                f = self._field._field
                p = self._field._p
                r = self._k
                q = r % p
                l = [f(q)]
                r -= q

                # invariant: Polynomial(l)(p) + r does not change, r multiple of p
                while r:
                    r //= p
                    q = r % p
                    l.append(f(q))
                    r -= q

                self._poly = Polynomial(l)

        def field(self):
            return self._field

        def inv(self):
            r, x, y = gcd(self._field._poly, self._poly)
            y = y // r
            return type(self)(y, self._field)

        def __truediv__(self, other):
            if type(other) is not type(self):
                return self / type(self)(other, self._field)
            return type(self).__mul__(self, other.inv())

        def __mul__(self, other):
            if type(other) is type(self):
                p = self._poly * other._poly
                q = p // self._field._poly
                r = p - q * self._field._poly
                return type(self)(r, self._field)
            if type(other) is Vector:
                return Vector([self * other[i] for i in range(other.dim)])
            return self * type(self)(other, self._field)

        def __rmul__(self, other):
            return self * other

        def __add__(self, other):
            if type(other) is type(self):
                return type(self)(self._poly + other._poly, self._field)
            return self + type(self)(other, self._field)

        def __radd__(self, other):
            return self + other

        def __neg__(self):
            return type(self)(-self._poly, self._field)

        def __sub__(self, other):
            return type(self)(self._poly - other._poly, self._field)

        def __eq__(self, other):
            if type(other) is type(self):
                return self._poly == other._poly
            return self._poly == other

        def __ne__(self, other):
            return not self == other

        def __repr__(self):
            return f"{self._k}"

    def __init__(self, field, n, poly):
        self._p = field._p
        self._n = n
        self._N = field._p**n
        assert type(n) is int
        assert n > 1
        assert poly[n] == 1
        self._poly = poly
        self._field = field

    def __call__(self, k: int):
        if type(k) == self.Element:
            return k
        return self.Element(k, self)

    def get_elements(self):
        return [self.Element(i, self) for i in range(self._N)]


class Polynomial:
    """
    Polynomial over a field.  a correponds to

       a[0] + a[1]*X + ... + a[n]*X^n

    where a[i] all belong to the same field and a[n] != 0 unless n = 0

    Usage:
          field = FiniteField_p(7)
          p = Polynomial([field(1), field(3)])
          q = Polynomial([field(0), field(1)])
          print(p)            # prints 3X + 1
          print(q)            # prints X
          print(p * q)        # prints 3X^2 + X
          print((p * q) // q) # prints 3X + 1
    """

    def __init__(self, a: list):
        self._a = a

        # prune trailing zeros
        while len(self._a) > 1 and self._a[-1] == 0:
            self._a = self._a[:-1]

        self.deg = len(self._a) - 1
        if self.deg == 0 and self._a[0] == 0:
            self.deg = -1

    def __floordiv__(self, other):
        zero = self._a[0] * 0
        one = zero + 1
        X = Polynomial([zero, one])
        if type(other) is not type(self):
            # other is field element
            inv_other = one / other
            return self * inv_other
        if self.deg < other.deg:
            return Polynomial([zero])

        a = self
        b = other
        q = Polynomial([zero])
        r = a

        # invariant is that a = bq + r
        while b.deg <= r.deg:
            n = r.deg - b.deg
            l = [zero] * (n + 1)
            l[n] = r._a[r.deg] / b._a[b.deg]
            p = Polynomial(l)
            q = q + p
            r = r - b * p
        return q

    def __call__(self, x):
        v = self._a[0]
        n = 1
        xn = x
        while n < len(self._a):
            v += self._a[n] * xn
            xn *= x
            n += 1
        return v

    def __add__(self, other):
        if type(other) is not type(self):
            zero = self._a[0] * 0
            one = zero + 1
            other = Polynomial([one * other])
        n = max(len(self._a), len(other._a))
        a = []
        for i in range(n):
            if i < len(self._a) and i < len(other._a):
                a.append(self._a[i] + other._a[i])
            elif i < len(self._a):
                a.append(self._a[i])
            else:
                a.append(other._a[i])
        return type(self)(a)

    def __mul__(self, other):
        if type(other) is type(self):
            a = self._a
            b = other._a
            m = len(a) - 1
            n = len(b) - 1
            c = []
            for k in range(m + n + 1):
                l = []
                for i in range(max(k - n, 0), min(m, k) + 1):
                    j = k - i
                    l.append(a[i] * b[j])
                c.append(reduce(lambda x, y: x + y, l))
            return type(self)(c)
        return type(self)([a * other for a in self._a])

    def __rmul__(self, other):
        return self * other

    def __neg__(self):
        return type(self)([-ai for ai in self._a])

    def __sub__(self, other):
        return self + (-other)

    def __getitem__(self, i):
        if i >= len(self._a):
            return self._a[0] * 0
        return self._a[i]

    def __eq__(self, other):
        if type(other) is type(self):
            if other.deg != self.deg:
                return False
            for i in range(self.deg):
                if other._a[i] != self._a[i]:
                    return False
            return True
        zero = self._a[0] * 0
        one = zero + 1
        other = one * other  # cast!!
        return self.deg == -1 and self._a[0] == other

    def __ne__(self, other):
        return not self == other

    def __repr__(self):
        arr = []
        n = len(self._a) - 1
        if n == 0:
            return repr(self._a[0])
        while n >= 0:
            a = self._a[n]
            power = f"X^{n}" if n > 1 else ("X" if n == 1 else "")
            if a != 0:
                if n > 0:
                    arr.append(f"{a}{power}" if a != 1 else f"{power}")
                else:
                    arr.append(f"{a}")
            n -= 1
        return " + ".join(arr)


class Vector:
    """
    Vector over a field.

    v[i] all belong to the same field

    Usage:
          field = FiniteField_p(7) # or FiniteField_p_n
          u = Vector([field(0), field(1), field(2)])
          print(4*u) # prints [0, 4, 1]
    """

    def __init__(self, v: list):
        self._v = v
        self.dim = len(self._v)

    def __getitem__(self, i):
        return self._v[i]

    def __repr__(self):
        return repr(self._v)

    def __mul__(self, other):
        return Vector([a * other for a in self._v])

    def __add__(self, other):
        return Vector([self._v[i] + other._v[i] for i in range(self.dim)])

    def __rmul__(self, other):
        return self * other


def get_2D_subspaces_of_3D_vector_space(field):
    # These correspond to all unique 2D subspaces of the 3D space over finite field

    subspaces = []

    # we can assume first basis vector starts with some leading 0s then a 1
    # we can also assume 2nd basis vector starts with more leading zeros then a 1
    # we can assume a 0 in the first basis vector where of the leading occurs in the 2nd basis vector

    # options where first basis vector starts [1,] and second starts [0,1,]
    for x in field.get_elements():
        for y in field.get_elements():
            u = Vector([field(1), field(0), field(x)])
            v = Vector([field(0), field(1), field(y)])
            basis = [u, v]
            subspaces.append(basis)

    # options where first basis vector starts [1,] and second starts [0,0,1]
    for x in field.get_elements():
        u = Vector([field(1), field(x), field(0)])
        v = Vector([field(0), field(0), field(1)])
        basis = [u, v]
        subspaces.append(basis)

    # options where first basis vector starts [0,1,]
    u = Vector([field(0), field(1), field(0)])
    v = Vector([field(0), field(0), field(1)])
    basis = [u, v]
    subspaces.append(basis)

    return subspaces


def get_1D_subspaces_of_2D_subspace_of_3D_vector_space(u, v, field):
    """
    u and v are 3D Vectors spanning the 2D subspace
    """

    # Every 1D subspace can be represented by a vector whose leading value is 1

    # These are all scalar multiples of u and v: e.g. a*u + b*v

    # To avoid double counting we should assume a is 1, or a,b is 0,1

    subspaces = []
    a = 1
    for b in field.get_elements():
        subspaces.append(a * u + b * v)

    a = 0
    b = 1
    subspaces.append(a * u + b * v)

    return subspaces


def get_dobble_cards(p, n):
    assert type(p) is int
    assert type(n) is int
    assert n > 0

    if n > 1:
        subfield = FiniteField_p(p)

        # find an irreducible polynomial - there must always be one due to combinatorics
        def find_irreducible_poly():
            a2 = subfield(1)
            for a1 in subfield.get_elements():
                for a0 in subfield.get_elements():
                    poly = Polynomial([a0, a1, a2])
                    values = [poly(x) for x in subfield.get_elements()]
                    is_irreducible = 0 not in values
                    if is_irreducible:
                        return poly

        irreducible_poly = find_irreducible_poly()
        print(f"Using polynomial {irreducible_poly} over F{p}")
        field = FiniteField_p_n(subfield, 2, irreducible_poly)
    else:
        field = FiniteField_p(p)

    _2D_subspaces = get_2D_subspaces_of_3D_vector_space(field)
    next_symbol = 0
    symbol_lookup = {}

    card_number = 0
    cards = []
    for card_number, _2D_subspace in enumerate(_2D_subspaces):
        u, v = _2D_subspace
        _1D_subspaces = get_1D_subspaces_of_2D_subspace_of_3D_vector_space(u, v, field)

        card = []
        for w in _1D_subspaces:
            s = repr(w)
            if s in symbol_lookup:
                card.append(symbol_lookup[s])
            else:
                card.append(next_symbol)
                symbol_lookup[s] = next_symbol
                next_symbol += 1
        card.sort()
        cards.append(card)

    return cards


if __name__ == "__main__":
    cards = get_dobble_cards(3, 2)
    for card in cards:
        print(card)

Comments