designs_pyx.pyx

r"""
Cython functions for combinatorial designs

This module implements the design methods that need to be somewhat efficient.

Functions
---------
"""

include "sage/data_structures/bitset.pxi"

from libc.string cimport memset

from cysignals.memory cimport sig_malloc, sig_calloc, sig_realloc, sig_free

from sage.misc.unknown import Unknown

def is_orthogonal_array(OA, int k, int n, int t=2, verbose=False, terminology="OA"):
    r"""
    Check that the integer matrix `OA` is an `OA(k,n,t)`.

    See :func:`~sage.combinat.designs.orthogonal_arrays.orthogonal_array`
    for a definition.

    INPUT:

    - ``OA`` -- the Orthogonal Array to be tested

    - ``k,n,t`` (integers) -- only implemented for `t=2`.

    - ``verbose`` (boolean) -- whether to display some information when ``OA``
      is not an orthogonal array `OA(k,n)`.

    - ``terminology`` (string) -- how to phrase the information when ``verbose =
      True``. Possible values are `"OA"`, `"MOLS"`.

    EXAMPLES::

        sage: from sage.combinat.designs.designs_pyx import is_orthogonal_array
        sage: OA = designs.orthogonal_arrays.build(8,9)
        sage: is_orthogonal_array(OA,8,9)
        True
        sage: is_orthogonal_array(OA,8,10)
        False
        sage: OA[4][3] = 1
        sage: is_orthogonal_array(OA,8,9)
        False
        sage: is_orthogonal_array(OA,8,9,verbose=True)
        Columns 0 and 3 are not orthogonal
        False
        sage: is_orthogonal_array(OA,8,9,verbose=True,terminology="MOLS")
        Squares 0 and 3 are not orthogonal
        False

    TESTS::

        sage: is_orthogonal_array(OA,8,9,t=3)
        Traceback (most recent call last):
        ...
        NotImplementedError: only implemented for t=2
        sage: is_orthogonal_array([[3]*8],8,9,verbose=True)
        The number of rows is 1 instead of 9^2=81
        False
        sage: is_orthogonal_array([[3]*8],8,9,verbose=True,terminology="MOLS")
        All squares do not have dimension n^2=9^2
        False
        sage: is_orthogonal_array([[3]*7],8,9,verbose=True)
        Some row does not have length 8
        False
        sage: is_orthogonal_array([[3]*7],8,9,verbose=True,terminology="MOLS")
        The number of squares is not 6
        False

    Up to relabelling, there is a unique `OA(3,2)`. So their number is just the
    cardinality of the relabeling group which is `S_2^3 \times S_3` and has
    cardinality `48`::

        sage: from itertools import product
        sage: n = 0
        sage: for a in product(product((0,1), repeat=3), repeat=4):
        ....:     if is_orthogonal_array(a,3,2):
        ....:          n += 1
        sage: n
        48
    """
    cdef int n2 = n*n
    cdef int x

    if t != 2:
        raise NotImplementedError("only implemented for t=2")

    for R in OA:
        if len(R) != k:
            if verbose:
                print({"OA"   : "Some row does not have length "+str(k),
                       "MOLS" : "The number of squares is not "+str(k-2)}[terminology])
            return False

    if len(OA) != n2:
        if verbose:
            print({"OA"   : "The number of rows is {} instead of {}^2={}".format(len(OA),n,n2),
                   "MOLS" : "All squares do not have dimension n^2={}^2".format(n)}[terminology])
        return False

    if n == 0:
        return True

    cdef int i,j,l

    # A copy of OA
    cdef unsigned short * OAc = <unsigned short *> sig_malloc(k*n2*sizeof(unsigned short))

    cdef unsigned short * C1
    cdef unsigned short * C2

    # failed malloc ?
    if OAc is NULL:
        raise MemoryError

    # Filling OAc
    for i,R in enumerate(OA):
        for j,x in enumerate(R):
            if x < 0 or x >= n:
                if verbose:
                    print({"OA"   : "{} is not in the interval [0..{}]".format(x,n-1),
                           "MOLS" : "Entry {} was expected to be in the interval [0..{}]".format(x,n-1)}[terminology])
                sig_free(OAc)
                return False
            OAc[j*n2+i] = x

    # A bitset to keep track of pairs of values
    cdef bitset_t seen
    bitset_init(seen, n2)

    for i in range(k): # For any column C1
        C1 = OAc+i*n2
        for j in range(i+1,k): # For any column C2 > C1
            C2 = OAc+j*n2
            bitset_set_first_n(seen, 0) # No pair has been seen yet
            for l in range(n2):
                bitset_add(seen,n*C1[l]+C2[l])

            if bitset_len(seen) != n2: # Have we seen all pairs ?
                sig_free(OAc)
                bitset_free(seen)
                if verbose:
                    print({"OA"   : "Columns {} and {} are not orthogonal".format(i,j),
                           "MOLS" : "Squares {} and {} are not orthogonal".format(i,j)}[terminology])
                return False

    sig_free(OAc)
    bitset_free(seen)
    return True

def is_group_divisible_design(groups,blocks,v,G=None,K=None,lambd=1,verbose=False):
    r"""
    Checks that input is a Group Divisible Design on `\{0,...,v-1\}`

    For more information on Group Divisible Designs, see
    :class:`~sage.combinat.designs.group_divisible_designs.GroupDivisibleDesign`.

    INPUT:

    - ``groups`` -- a partition of `X`. If set to ``None`` the groups are
      guessed automatically, and the function returns ``(True, guessed_groups)``
      instead of ``True``

    - ``blocks`` -- collection of blocks

    - ``v`` (integers) -- size of the ground set assumed to be `X=\{0,...,v-1\}`.

    - ``G`` -- list of integers of which the sizes of the groups must be
      elements. Set to ``None`` (automatic guess) by default.

    - ``K`` -- list of integers of which the sizes of the blocks must be
      elements. Set to ``None`` (automatic guess) by default.

    - ``lambd`` -- value of `\lambda`. Set to `1` by default.

    - ``verbose`` (boolean) -- whether to display some information when the
      design is not a GDD.

    EXAMPLES::

        sage: from sage.combinat.designs.designs_pyx import is_group_divisible_design
        sage: TD = designs.transversal_design(4,10)
        sage: groups = [list(range(i*10,(i+1)*10)) for i in range(4)]
        sage: is_group_divisible_design(groups,TD,40,lambd=1)
        True

    TESTS::

        sage: TD = designs.transversal_design(4,10)
        sage: groups = [list(range(i*10,(i+1)*10)) for i in range(4)]
        sage: is_group_divisible_design(groups,TD,40,lambd=2,verbose=True)
        the pair (0,10) has been seen 1 times but lambda=2
        False
        sage: is_group_divisible_design([[1,2],[3,4]],[[1,2]],40,lambd=1,verbose=True)
        groups is not a partition of [0,...,39]
        False
        sage: is_group_divisible_design([list(range(40))],[[1,2]],40,lambd=1,verbose=True)
        the pair (1,2) belongs to a group but appears in some block
        False
        sage: is_group_divisible_design([list(range(40))],[[2,2]],40,lambd=1,verbose=True)
        The following block has repeated elements: [2, 2]
        False
        sage: is_group_divisible_design([list(range(40))],[["e",2]],40,lambd=1,verbose=True)
        e does not belong to [0,...,39]
        False
        sage: is_group_divisible_design([list(range(40))],[list(range(40))],40,G=[5],lambd=1,verbose=True)
        a group has size 40 while G=[5]
        False
        sage: is_group_divisible_design([list(range(40))],[["e",2]],40,K=[1],lambd=1,verbose=True)
        a block has size 2 while K=[1]
        False

        sage: p = designs.projective_plane(3)
        sage: is_group_divisible_design(None, p.blocks(), 13)
        (True, [[0], [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12]])
        sage: is_group_divisible_design(None, p.blocks()*2, 13, verbose=True)
        the pair (0,1) has been seen 2 times but lambda=1
        False
        sage: is_group_divisible_design(None, p.blocks()*2, 13, lambd=2)
        (True, [[0], [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12]])
    """
    cdef int n = v
    cdef int i,ii,j,jj,s,isok
    cdef int l = lambd
    cdef bint guess_groups = groups is None

    if v < 0 or lambd < 0:
        if verbose:
            print("v={} and lambda={} must be non-negative integers".format(v,l))
        return False

    # Block sizes are element of K
    if K is not None:
        K = set(K)
        for b in blocks:
            if not len(b) in K:
                if verbose:
                    print("a block has size {} while K={}".format(len(b),list(K)))
                return False

    # Check that "groups" consists of disjoint sets whose union has length n
    if (groups is not None and
        (sum(len(g) for g in groups) != n or
         len(set().union(*groups)) != n)):
        if verbose:
            print("groups is not a partition of [0,...,{}]".format(n-1))
        return False

    # Checks that the blocks are indeed sets and do not repeat elements
    for b in blocks:
        if len(b) != len(set(b)):
            if verbose:
                print("The following block has repeated elements: {}".format(b))
            return False

    # Check that the groups/blocks belong to [0,...,n-1]
    from itertools import chain
    for b in chain(groups if groups is not None else [],blocks):
        for x in b:
            try:
                i = x
            except TypeError:
                i = -1
            if i < 0 or i >= n:
                if verbose:
                    print("{} does not belong to [0,...,{}]".format(x, n-1))
                return False

    cdef unsigned short * matrix = <unsigned short *> sig_calloc(n*n,sizeof(unsigned short))
    if matrix is NULL:
        raise MemoryError

    # Counts the number of occurrences of each pair of points
    for b in blocks:
        s = len(b)
        for i in range(s):
            ii = b[i]
            for j in range(i+1,s):
                jj = b[j]
                matrix[ii*n+jj] += 1
                matrix[jj*n+ii] += 1

    # Guess the groups (if necessary)
    if groups is None:
        from sage.sets.disjoint_set import DisjointSet_of_integers
        groups = DisjointSet_of_integers(n)
        for i in range(n):
            for j in range(i + 1, n):
                if matrix[i * n + j] == 0:
                    groups.union(i, j)
        groups = list(groups.root_to_elements_dict().values())

    # Group sizes are element of G
    if G is not None:
        G = set(G)
        for g in groups:
            if not len(g) in G:
                if verbose:
                    print("a group has size {} while G={}".format(len(g),list(G)))
                sig_free(matrix)
                return False

    # Checks that two points of the same group were never covered
    for g in groups:
        s = len(g)
        for i in range(s):
            ii = g[i]
            for j in range(i+1,s):
                jj = g[j]
                if matrix[ii*n+jj] != 0:
                    if verbose:
                        print("the pair ({},{}) belongs to a group but appears in some block".format(ii, jj))
                    sig_free(matrix)
                    return False

                # We fill the entries with what is expected by the next loop
                matrix[ii*n+jj] = l
                matrix[jj*n+ii] = l

    # Checking that what should be equal to lambda IS equal to lambda
    for i in range(n):
        for j in range(i+1,n):
            if matrix[i*n+j] != l:
                if verbose:
                    print("the pair ({},{}) has been seen {} times but lambda={}".format(i,j,matrix[i*n+j],l))
                sig_free(matrix)
                return False

    sig_free(matrix)

    return True if not guess_groups else (True, groups)


def is_pairwise_balanced_design(blocks,v,K=None,lambd=1,verbose=False):
    r"""
    Checks that input is a Pairwise Balanced Design (PBD) on `\{0,...,v-1\}`

    For more information on Pairwise Balanced Designs (PBD), see
    :class:`~sage.combinat.designs.bibd.PairwiseBalancedDesign`.

    INPUT:

    - ``blocks`` -- collection of blocks

    - ``v`` (integers) -- size of the ground set assumed to be `X=\{0,...,v-1\}`.

    - ``K`` -- list of integers of which the sizes of the blocks must be
      elements. Set to ``None`` (automatic guess) by default.

    - ``lambd`` -- value of `\lambda`. Set to `1` by default.

    - ``verbose`` (boolean) -- whether to display some information when the
      design is not a PBD.

    EXAMPLES::

        sage: from sage.combinat.designs.designs_pyx import is_pairwise_balanced_design
        sage: sts = designs.steiner_triple_system(9)
        sage: is_pairwise_balanced_design(sts,9,[3],1)
        True
        sage: TD = designs.transversal_design(4,10).blocks()
        sage: groups = [list(range(i*10,(i+1)*10)) for i in range(4)]
        sage: is_pairwise_balanced_design(TD+groups,40,[4,10],1,verbose=True)
        True

    TESTS::

        sage: from sage.combinat.designs.designs_pyx import is_pairwise_balanced_design
        sage: is_pairwise_balanced_design(TD+groups,40,[4,10],2,verbose=True)
        the pair (0,1) has been seen 1 times but lambda=2
        False
        sage: is_pairwise_balanced_design(TD+groups,40,[10],1,verbose=True)
        a block has size 4 while K=[10]
        False
        sage: is_pairwise_balanced_design([[2,2]],40,[2],1,verbose=True)
        The following block has repeated elements: [2, 2]
        False
        sage: is_pairwise_balanced_design([["e",2]],40,[2],1,verbose=True)
        e does not belong to [0,...,39]
        False
    """
    return is_group_divisible_design([[i] for i in range(v)],
                                     blocks,
                                     v,
                                     K=K,
                                     lambd=lambd,
                                     verbose=verbose)

def is_projective_plane(blocks, verbose=False):
    r"""
    Test whether the blocks form a projective plane on `\{0,...,v-1\}`

    A *projective plane* is an incidence structure that has the following properties:

    1. Given any two distinct points, there is exactly one line incident with both of them.
    2. Given any two distinct lines, there is exactly one point incident with both of them.
    3. There are four points such that no line is incident with more than two of them.

    For more informations, see :wikipedia:`Projective_plane`.

    :meth:`~IncidenceStructure.is_t_design` can also check if an incidence structure is a projective plane
    with the parameters `v=k^2+k+1`, `t=2` and `l=1`.

    INPUT:

    - ``blocks`` -- collection of blocks

    - ``verbose`` -- whether to print additional information


    EXAMPLES::

        sage: from sage.combinat.designs.designs_pyx import is_projective_plane
        sage: p = designs.projective_plane(4)
        sage: b = p.blocks()
        sage: is_projective_plane(b, verbose=True)
        True

        sage: p = designs.projective_plane(2)
        sage: b = p.blocks()
        sage: is_projective_plane(b)
        True
        sage: b[0][2] = 5
        sage: is_projective_plane(b, verbose=True)
        the pair (0,5) has been seen 2 times but lambda=1
        False

        sage: is_projective_plane([[0,1,2],[1,2,4]], verbose=True)
        the pair (0,3) has been seen 0 times but lambda=1
        False

        sage: is_projective_plane([[1]], verbose=True)
        First block has less than 3 points.
        False

        sage: p = designs.projective_plane(2)
        sage: b = p.blocks()
        sage: b[2].append(4)
        sage: is_projective_plane(b, verbose=True)
        a block has size 4 while K=[3]
        False
    """
    if not blocks:
        if verbose:
            print('There is no block.')
        return False
    k = len(blocks[0])-1
    if k < 2:
        if verbose:
            print('First block has less than 3 points.')
        return False
    v = k**2 + k + 1
    return is_group_divisible_design([[i] for i in range(v)],
                                     blocks,
                                     v,
                                     K=[k+1],
                                     lambd=1,
                                     verbose=verbose)

def is_difference_matrix(M,G,k,lmbda=1,verbose=False):
    r"""
    Test if `M` is a `(G,k,\lambda)`-difference matrix.

    A matrix `M` is a `(G,k,\lambda)`-difference matrix if its entries are
    element of `G`, and if for any two rows `R,R'` of `M` and `x\in G` there
    are exactly `\lambda` values `i` such that `R_i-R'_i=x`.

    INPUT:

    - ``M`` -- a matrix with entries from ``G``

    - ``G`` -- a group

    - ``k`` -- integer

    - ``lmbda`` (integer) -- set to `1` by default.

    - ``verbose`` (boolean) -- whether to print some information when the answer
      is ``False``.

    EXAMPLES::

        sage: from sage.combinat.designs.designs_pyx import is_difference_matrix
        sage: q = 3**3
        sage: F = GF(q,'x')
        sage: M = [[x*y for y in F] for x in F]
        sage: is_difference_matrix(M,F,q,verbose=1)
        True

        sage: B = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        ....:      [0, 1, 2, 3, 4, 2, 3, 4, 0, 1],
        ....:      [0, 2, 4, 1, 3, 3, 0, 2, 4, 1]]
        sage: G = GF(5)
        sage: B = [[G(b) for b in R] for R in B]
        sage: is_difference_matrix(list(zip(*B)),G,3,2)
        True

    Bad input::

        sage: for R in M: R.append(None)
        sage: is_difference_matrix(M,F,q,verbose=1)
        The matrix has 28 columns but k=27
        False
        sage: for R in M: _=R.pop(-1)
        sage: M.append([None]*3**3)
        sage: is_difference_matrix(M,F,q,verbose=1)
        The matrix has 28 rows instead of lambda(|G|-1+2u)+mu=1(27-1+2.0)+1=27
        False
        sage: _= M.pop(-1)
        sage: for R in M: R[-1] = 0
        sage: is_difference_matrix(M,F,q,verbose=1)
        Columns 0 and 26 generate 0 exactly 27 times instead of the expected mu(=1)
        False
        sage: for R in M: R[-1] = 1
        sage: M[-1][-1] = 0
        sage: is_difference_matrix(M,F,q,verbose=1)
        Columns 0 and 26 do not generate all elements of G exactly lambda(=1) times. The element x appeared 0 times as a difference.
        False
    """
    return is_quasi_difference_matrix(M,G,k,lmbda=lmbda,mu=lmbda,u=0,verbose=verbose)

def is_quasi_difference_matrix(M,G,int k,int lmbda,int mu,int u,verbose=False):
    r"""
    Test if the matrix is a `(G,k;\lambda,\mu;u)`-quasi-difference matrix

    Let `G` be an abelian group of order `n`. A
    `(n,k;\lambda,\mu;u)`-quasi-difference matrix (QDM) is a matrix `Q_{ij}`
    with `\lambda(n-1+2u)+\mu` rows and `k` columns, with each entry either
    equal to ``None`` (i.e. the 'missing entries') or to an element of `G`. Each
    column contains exactly `\lambda u` empty entries, and each row contains at
    most one ``None``. Furthermore, for each `1\leq i<j\leq k`, the multiset

    .. MATH::

        \{q_{li}-q_{lj}:1\leq l\leq \lambda (n-1+2u)+\mu, \text{ with } q_{li}\text{ and }q_{lj}\text{ not empty}\}

    contains `\lambda` times every nonzero element of `G` and contains `\mu`
    times `0`.

    INPUT:

    - ``M`` -- a matrix with entries from ``G`` (or equal to ``None`` for
      missing entries)

    - ``G`` -- a group

    - ``k,lmbda,mu,u`` -- integers

    - ``verbose`` (boolean) -- whether to print some information when the answer
      is ``False``.

    EXAMPLES:

    Differences matrices::

        sage: from sage.combinat.designs.designs_pyx import is_quasi_difference_matrix
        sage: q = 3**3
        sage: F = GF(q,'x')
        sage: M = [[x*y for y in F] for x in F]
        sage: is_quasi_difference_matrix(M,F,q,1,1,0,verbose=1)
        True

        sage: B = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        ....:      [0, 1, 2, 3, 4, 2, 3, 4, 0, 1],
        ....:      [0, 2, 4, 1, 3, 3, 0, 2, 4, 1]]
        sage: G = GF(5)
        sage: B = [[G(b) for b in R] for R in B]
        sage: is_quasi_difference_matrix(list(zip(*B)),G,3,2,2,0)
        True

    A quasi-difference matrix from the database::

        sage: from sage.combinat.designs.database import QDM
        sage: G,M = QDM[38,1][37,1,1,1][1]()
        sage: is_quasi_difference_matrix(M,G,k=6,lmbda=1,mu=1,u=1)
        True

    Bad input::

        sage: is_quasi_difference_matrix(M,G,k=6,lmbda=1,mu=1,u=3,verbose=1)
        The matrix has 39 rows instead of lambda(|G|-1+2u)+mu=1(37-1+2.3)+1=43
        False
        sage: is_quasi_difference_matrix(M,G,k=6,lmbda=1,mu=2,u=1,verbose=1)
        The matrix has 39 rows instead of lambda(|G|-1+2u)+mu=1(37-1+2.1)+2=40
        False
        sage: M[3][1] = None
        sage: is_quasi_difference_matrix(M,G,k=6,lmbda=1,mu=1,u=1,verbose=1)
        Row 3 contains more than one empty entry
        False
        sage: M[3][1] = 1
        sage: M[6][1] = None
        sage: is_quasi_difference_matrix(M,G,k=6,lmbda=1,mu=1,u=1,verbose=1)
        Column 1 contains 2 empty entries instead of the expected lambda.u=1.1=1
        False
    """
    from .difference_family import group_law

    assert k>=2
    assert lmbda >=1
    assert mu>=0
    assert u>=0

    cdef int n = G.cardinality()
    cdef int M_nrows = len(M)
    cdef int i,j,ii
    cdef bint bit

    # Height of the matrix
    if lmbda*(n-1+2*u)+mu != M_nrows:
        if verbose:
            print("The matrix has {} rows instead of lambda(|G|-1+2u)+mu={}({}-1+2.{})+{}={}".format(M_nrows,lmbda,n,u,mu,lmbda*(n-1+2*u)+mu))
        return False

    # Width of the matrix
    for R in M:
        if len(R)!=k:
            if verbose:
                print("The matrix has {} columns but k={}".format(len(R),k))
            return False

    # When |G|=0
    if M_nrows == 0:
        return True

    # Map group element with integers
    cdef list int_to_group = list(G)
    cdef dict group_to_int = {v:i for i,v in enumerate(int_to_group)}

    # Allocations
    cdef int ** x_minus_y     = <int **> sig_malloc((n+1)*sizeof(int *))
    cdef int * x_minus_y_data = <int *>  sig_malloc((n+1)*(n+1)*sizeof(int))
    cdef int * M_c            = <int *>  sig_malloc(k*M_nrows*sizeof(int))
    cdef int * G_seen         = <int *>  sig_malloc((n+1)*sizeof(int))
    if (x_minus_y == NULL or x_minus_y_data == NULL or M_c == NULL or G_seen == NULL):
        sig_free(x_minus_y)
        sig_free(x_minus_y_data)
        sig_free(G_seen)
        sig_free(M_c)
        raise MemoryError

    # The "x-y" table. If g_i, g_j \in G, then x_minus_y[i][j] is equal to
    # group_to_int[g_i-g_j].
    #
    # In order to handle empty values represented by n, we have
    # x_minus_y[?][n]=x_minus_y[n][?]=n
    zero, op, inv = group_law(G)
    x_minus_y[0] = x_minus_y_data
    for i in range(1,n+1):
        x_minus_y[i] = x_minus_y[i-1] + n+1

    # Elements of G
    for j,Gj in enumerate(int_to_group):
        minus_Gj = inv(Gj)
        assert op(Gj, minus_Gj) == zero
        for i,Gi in enumerate(int_to_group):
            x_minus_y[i][j] = group_to_int[op(Gi,minus_Gj)]

    # Empty values
    for i in range(n+1):
        x_minus_y[n][i]=n
        x_minus_y[i][n]=n

    # A copy of the matrix
    for i,R in enumerate(M):
        for j,x in enumerate(R):
            M_c[i*k+j] = group_to_int[G(x)] if x is not None else n

    # Each row contains at most one empty entry
    if u:
        for i in range(M_nrows):
            bit = False
            for j in range(k):
                if M_c[i*k+j] == n:
                    if bit:
                        if verbose:
                            print("Row {} contains more than one empty entry".format(i))
                        sig_free(x_minus_y_data)
                        sig_free(x_minus_y)
                        sig_free(G_seen)
                        sig_free(M_c)
                        return False
                    bit = True

    # Each column contains lmbda*u empty entries
        for j in range(k):
            ii = 0
            for i in range(M_nrows):
                if M_c[i*k+j] == n:
                    ii += 1
            if ii!=lmbda*u:
                if verbose:
                    print("Column {} contains {} empty entries instead of the expected "
                          "lambda.u={}.{}={}".format(j, ii, lmbda, u, lmbda*u))
                sig_free(x_minus_y_data)
                sig_free(x_minus_y)
                sig_free(G_seen)
                sig_free(M_c)
                return False

    # We are now ready to test every pair of columns
    for i in range(k):
        for j in range(i+1,k):
            memset(G_seen, 0, (n+1)*sizeof(int))
            for ii in range(M_nrows):
                G_seen[x_minus_y[M_c[ii*k+i]][M_c[ii*k+j]]] += 1

            if G_seen[0] != mu: # Bad number of 0
                if verbose:
                    print("Columns {} and {} generate 0 exactly {} times "
                          "instead of the expected mu(={})".format(i,j,G_seen[0],mu))
                sig_free(x_minus_y_data)
                sig_free(x_minus_y)
                sig_free(G_seen)
                sig_free(M_c)
                return False

            for ii in range(1,n): # bad number of g_ii\in G
                if G_seen[ii] != lmbda:
                    if verbose:
                        print("Columns {} and {} do not generate all elements of G "
                         "exactly lambda(={}) times. The element {} appeared {} "
                         "times as a difference.".format(i,j,lmbda,int_to_group[ii],G_seen[ii]))
                    sig_free(x_minus_y_data)
                    sig_free(x_minus_y)
                    sig_free(G_seen)
                    sig_free(M_c)
                    return False

    sig_free(x_minus_y_data)
    sig_free(x_minus_y)
    sig_free(G_seen)
    sig_free(M_c)
    return True

# Cached information for OA constructions (see .pxd file for more info)

_OA_cache = <cache_entry *> sig_malloc(2*sizeof(cache_entry))
if (_OA_cache == NULL):
    sig_free(_OA_cache)
    raise MemoryError
_OA_cache[0].max_true = -1
_OA_cache[1].max_true = -1
_OA_cache_size = 2

cpdef _OA_cache_set(int k,int n,truth_value):
    r"""
    Sets a value in the OA cache of existence results

    INPUT:

    - ``k,n`` (integers)

    - ``truth_value`` -- one of ``True,False,Unknown``
    """
    global _OA_cache, _OA_cache_size
    cdef int i
    if _OA_cache_size <= n:
        new_cache_size = n+100
        _OA_cache = <cache_entry *> sig_realloc(_OA_cache,new_cache_size*sizeof(cache_entry))
        if _OA_cache == NULL:
            sig_free(_OA_cache)
            raise MemoryError

        for i in range(_OA_cache_size,new_cache_size):
            _OA_cache[i].max_true = 0
            _OA_cache[i].min_unknown = -1
            _OA_cache[i].max_unknown = 0
            _OA_cache[i].min_false = -1

        _OA_cache_size = new_cache_size

    if truth_value is True:
        _OA_cache[n].max_true    = k if k>_OA_cache[n].max_true    else _OA_cache[n].max_true
    elif truth_value is Unknown:
        _OA_cache[n].min_unknown = k if k<_OA_cache[n].min_unknown else _OA_cache[n].min_unknown
        _OA_cache[n].max_unknown = k if k>_OA_cache[n].max_unknown else _OA_cache[n].max_unknown
    else:
        _OA_cache[n].min_false   = k if k<_OA_cache[n].min_false   else _OA_cache[n].min_false

cpdef _OA_cache_get(int k,int n):
    r"""
    Gets a value from the OA cache of existence results

    INPUT:

    ``k,n`` (integers)
    """
    if n>=_OA_cache_size:
        return None
    if k <= _OA_cache[n].max_true:
        return True
    elif (k >= _OA_cache[n].min_unknown and k <= _OA_cache[n].max_unknown):
        return Unknown
    elif k >= _OA_cache[n].min_false:
        return False

    return None

cpdef _OA_cache_construction_available(int k,int n):
    r"""
    Tests if a construction is implemented using the cache's information

    INPUT:

    - ``k,n`` (integers)
    """
    if n>=_OA_cache_size:
        return Unknown
    if k <= _OA_cache[n].max_true:
        return True
    if k >= _OA_cache[n].min_unknown:
        return False
    else:
        return Unknown