Coverage for local/lib/python2.7/site-packages/sage/combinat/designs/bibd.py : 82%

r""" 

Balanced Incomplete Block Designs (BIBD) 

This module gathers everything related to Balanced Incomplete Block Designs. One can build a 

BIBD (or check that it can be built) with :func:`balanced_incomplete_block_design`:: 

sage: BIBD = designs.balanced_incomplete_block_design(7,3) 

In particular, Sage can build a `(v,k,1)`-BIBD when one exists for all `k\leq 

5`. The following functions are available: 

.. csv-table:: 

:class: contentstable 

:widths: 30, 70 

:delim: | 

:func:`balanced_incomplete_block_design` | Return a BIBD of parameters `v,k`. 

:func:`BIBD_from_TD` | Return a BIBD through TD-based constructions. 

:func:`BIBD_from_difference_family` | Return the BIBD associated to the difference family ``D`` on the group ``G``. 

:func:`BIBD_from_PBD` | Return a `(v,k,1)`-BIBD from a `(r,K)`-PBD where `r=(v-1)/(k-1)`. 

:func:`PBD_from_TD` | Return a `(kt,\{k,t\})`-PBD if `u=0` and a `(kt+u,\{k,k+1,t,u\})`-PBD otherwise. 

:func:`steiner_triple_system` | Return a Steiner Triple System. 

:func:`v_5_1_BIBD` | Return a `(v,5,1)`-BIBD. 

:func:`v_4_1_BIBD` | Return a `(v,4,1)`-BIBD. 

:func:`PBD_4_5_8_9_12` | Return a `(v,\{4,5,8,9,12\})`-PBD on `v` elements. 

:func:`BIBD_5q_5_for_q_prime_power` | Return a `(5q,5,1)`-BIBD with `q\equiv 1\pmod 4` a prime power. 

**Construction of BIBD when** `k=4` 

Decompositions of `K_v` into `K_4` (i.e. `(v,4,1)`-BIBD) are built following 

Douglas Stinson's construction as presented in [Stinson2004]_ page 167. It is 

based upon the construction of `(v\{4,5,8,9,12\})`-PBD (see the doc of 

:func:`PBD_4_5_8_9_12`), knowing that a `(v\{4,5,8,9,12\})`-PBD on `v` points 

can always be transformed into a `((k-1)v+1,4,1)`-BIBD, which covers all 

possible cases of `(v,4,1)`-BIBD. 

**Construction of BIBD when** `k=5` 

Decompositions of `K_v` into `K_4` (i.e. `(v,4,1)`-BIBD) are built following 

Clayton Smith's construction [ClaytonSmith]_. 

.. [ClaytonSmith] On the existence of `(v,5,1)`-BIBD. 

http://www.argilo.net/files/bibd.pdf 

Clayton Smith 

Functions 

--------- 

""" 

# python3 

from __future__ import division, print_function 

from __future__ import absolute_import 

from six.moves import range 

from sage.categories.sets_cat import EmptySetError 

from sage.misc.unknown import Unknown 

from .design_catalog import transversal_design 

from .block_design import BlockDesign 

from sage.arith.all import binomial, is_prime_power 

from .group_divisible_designs import GroupDivisibleDesign 

from .designs_pyx import is_pairwise_balanced_design 

def balanced_incomplete_block_design(v, k, existence=False, use_LJCR=False): 

r""" 

Return a BIBD of parameters `v,k`. 

A Balanced Incomplete Block Design of parameters `v,k` is a collection 

`\mathcal C` of `k`-subsets of `V=\{0,\dots,v-1\}` such that for any two 

distinct elements `x,y\in V` there is a unique element `S\in \mathcal C` 

such that `x,y\in S`. 

More general definitions sometimes involve a `\lambda` parameter, and we 

assume here that `\lambda=1`. 

For more information on BIBD, see the 

:wikipedia:`corresponding Wikipedia entry <Block_design#Definition_of_a_BIBD_.28or_2-design.29>`. 

INPUT: 

- ``v,k`` (integers) 

- ``existence`` (boolean) -- instead of building the design, return: 

- ``True`` -- meaning that Sage knows how to build the design 

- ``Unknown`` -- meaning that Sage does not know how to build the 

design, but that the design may exist (see :mod:`sage.misc.unknown`). 

- ``False`` -- meaning that the design does not exist. 

- ``use_LJCR`` (boolean) -- whether to query the La Jolla Covering 

Repository for the design when Sage does not know how to build it (see 

:func:`~sage.combinat.designs.covering_design.best_known_covering_design_www`). This 

requires internet. 

.. SEEALSO:: 

* :func:`steiner_triple_system` 

* :func:`v_4_1_BIBD` 

* :func:`v_5_1_BIBD` 

.. TODO:: 

Implement other constructions from the Handbook of Combinatorial 

Designs. 

EXAMPLES:: 

sage: designs.balanced_incomplete_block_design(7, 3).blocks() 

[[0, 1, 3], [0, 2, 4], [0, 5, 6], [1, 2, 6], [1, 4, 5], [2, 3, 5], [3, 4, 6]] 

sage: B = designs.balanced_incomplete_block_design(66, 6, use_LJCR=True) # optional - internet 

sage: B # optional - internet 

Incidence structure with 66 points and 143 blocks 

sage: B.blocks() # optional - internet 

[[0, 1, 2, 3, 4, 65], [0, 5, 24, 25, 39, 57], [0, 6, 27, 38, 44, 55], ... 

sage: designs.balanced_incomplete_block_design(66, 6, use_LJCR=True) # optional - internet 

Incidence structure with 66 points and 143 blocks 

sage: designs.balanced_incomplete_block_design(216, 6) 

Traceback (most recent call last): 

... 

NotImplementedError: I don't know how to build a (216,6,1)-BIBD! 

TESTS:: 

sage: designs.balanced_incomplete_block_design(85,5,existence=True) 

True 

sage: _ = designs.balanced_incomplete_block_design(85,5) 

A BIBD from a Finite Projective Plane:: 

sage: _ = designs.balanced_incomplete_block_design(21,5) 

Some trivial BIBD:: 

sage: designs.balanced_incomplete_block_design(10,10) 

(10,10,1)-Balanced Incomplete Block Design 

sage: designs.balanced_incomplete_block_design(1,10) 

(1,0,1)-Balanced Incomplete Block Design 

Existence of BIBD with `k=3,4,5`:: 

sage: [v for v in range(50) if designs.balanced_incomplete_block_design(v,3,existence=True)] 

[1, 3, 7, 9, 13, 15, 19, 21, 25, 27, 31, 33, 37, 39, 43, 45, 49] 

sage: [v for v in range(100) if designs.balanced_incomplete_block_design(v,4,existence=True)] 

[1, 4, 13, 16, 25, 28, 37, 40, 49, 52, 61, 64, 73, 76, 85, 88, 97] 

sage: [v for v in range(150) if designs.balanced_incomplete_block_design(v,5,existence=True)] 

[1, 5, 21, 25, 41, 45, 61, 65, 81, 85, 101, 105, 121, 125, 141, 145] 

For `k > 5` there are currently very few constructions:: 

sage: [v for v in range(300) if designs.balanced_incomplete_block_design(v,6,existence=True) is True] 

[1, 6, 31, 66, 76, 91, 96, 106, 111, 121, 126, 136, 141, 151, 156, 171, 181, 186, 196, 201, 211, 241, 271] 

sage: [v for v in range(300) if designs.balanced_incomplete_block_design(v,6,existence=True) is Unknown] 

[51, 61, 81, 166, 216, 226, 231, 246, 256, 261, 276, 286, 291] 

Here are some constructions with `k \geq 7` and `v` a prime power:: 

sage: designs.balanced_incomplete_block_design(169,7) 

(169,7,1)-Balanced Incomplete Block Design 

sage: designs.balanced_incomplete_block_design(617,8) 

(617,8,1)-Balanced Incomplete Block Design 

sage: designs.balanced_incomplete_block_design(433,9) 

(433,9,1)-Balanced Incomplete Block Design 

sage: designs.balanced_incomplete_block_design(1171,10) 

(1171,10,1)-Balanced Incomplete Block Design 

And we know some inexistence results:: 

sage: designs.balanced_incomplete_block_design(21,6,existence=True) 

False 

""" 

lmbd = 1 

# Trivial BIBD 

if v == 1: 

if existence: 

return True 

return BalancedIncompleteBlockDesign(v, [], check=False) 

if k == v: 

if existence: 

return True 

return BalancedIncompleteBlockDesign(v, [list(range(v))], check=False, copy=False) 

# Non-existence of BIBD 

if (v < k or 

k < 2 or 

(v-1) % (k-1) != 0 or 

(v*(v-1)) % (k*(k-1)) != 0 or 

# From the Handbook of combinatorial designs: 

# 

# With lambda>1 other exceptions are 

# (15,5,2),(21,6,2),(22,7,2),(22,8,4). 

(k==6 and v in [36,46]) or 

(k==7 and v == 43) or 

# Fisher's inequality 

(v*(v-1))/(k*(k-1)) < v): 

if existence: 

return False 

raise EmptySetError("There exists no ({},{},{})-BIBD".format(v,k,lmbd)) 

if k == 2: 

if existence: 

return True 

from itertools import combinations 

return BalancedIncompleteBlockDesign(v, combinations(list(range(v)),2), check=False, copy=True) 

if k == 3: 

if existence: 

return v%6 == 1 or v%6 == 3 

return steiner_triple_system(v) 

if k == 4: 

if existence: 

return v%12 == 1 or v%12 == 4 

return BalancedIncompleteBlockDesign(v, v_4_1_BIBD(v), copy=False) 

if k == 5: 

if existence: 

return v%20 == 1 or v%20 == 5 

return BalancedIncompleteBlockDesign(v, v_5_1_BIBD(v), copy=False) 

from .difference_family import difference_family 

from .database import BIBD_constructions 

if (v,k,1) in BIBD_constructions: 

if existence: 

return True 

return BlockDesign(v,BIBD_constructions[(v,k,1)](), copy=False) 

if BIBD_from_arc_in_desarguesian_projective_plane(v,k,existence=True): 

if existence: 

return True 

B = BIBD_from_arc_in_desarguesian_projective_plane(v,k) 

return BalancedIncompleteBlockDesign(v, B, copy=False) 

if BIBD_from_TD(v,k,existence=True): 

if existence: 

return True 

return BalancedIncompleteBlockDesign(v, BIBD_from_TD(v,k), copy=False) 

if v == (k-1)**2+k and is_prime_power(k-1): 

if existence: 

return True 

from .block_design import projective_plane 

return BalancedIncompleteBlockDesign(v, projective_plane(k-1),copy=False) 

if difference_family(v,k,existence=True): 

if existence: 

return True 

G,D = difference_family(v,k) 

return BalancedIncompleteBlockDesign(v, BIBD_from_difference_family(G,D,check=False), copy=False) 

if use_LJCR: 

from .covering_design import best_known_covering_design_www 

B = best_known_covering_design_www(v,k,2) 

# Is it a BIBD or just a good covering ? 

expected_n_of_blocks = binomial(v,2)//binomial(k,2) 

if B.low_bd() > expected_n_of_blocks: 

if existence: 

return False 

raise EmptySetError("There exists no ({},{},{})-BIBD".format(v,k,lmbd)) 

B = B.incidence_structure() 

if B.num_blocks() == expected_n_of_blocks: 

if existence: 

return True 

else: 

return B 

if existence: 

return Unknown 

else: 

raise NotImplementedError("I don't know how to build a ({},{},1)-BIBD!".format(v,k)) 

def steiner_triple_system(n): 

r""" 

Return a Steiner Triple System 

A Steiner Triple System (STS) of a set `\{0,...,n-1\}` 

is a family `S` of 3-sets such that for any `i \not = j` 

there exists exactly one set of `S` in which they are 

both contained. 

It can alternatively be thought of as a factorization of 

the complete graph `K_n` with triangles. 

A Steiner Triple System of a `n`-set exists if and only if 

`n \equiv 1 \pmod 6` or `n \equiv 3 \pmod 6`, in which case 

one can be found through Bose's and Skolem's constructions, 

respectively [AndHonk97]_. 

INPUT: 

- ``n`` return a Steiner Triple System of `\{0,...,n-1\}` 

EXAMPLES: 

A Steiner Triple System on `9` elements :: 

sage: sts = designs.steiner_triple_system(9) 

sage: sts 

(9,3,1)-Balanced Incomplete Block Design 

sage: list(sts) 

[[0, 1, 5], [0, 2, 4], [0, 3, 6], [0, 7, 8], [1, 2, 3], 

[1, 4, 7], [1, 6, 8], [2, 5, 8], [2, 6, 7], [3, 4, 8], 

[3, 5, 7], [4, 5, 6]] 

As any pair of vertices is covered once, its parameters are :: 

sage: sts.is_t_design(return_parameters=True) 

(True, (2, 9, 3, 1)) 

An exception is raised for invalid values of ``n`` :: 

sage: designs.steiner_triple_system(10) 

Traceback (most recent call last): 

... 

EmptySetError: Steiner triple systems only exist for n = 1 mod 6 or n = 3 mod 6 

REFERENCE: 

.. [AndHonk97] A short course in Combinatorial Designs, 

Ian Anderson, Iiro Honkala, 

Internet Editions, Spring 1997, 

http://www.utu.fi/~honkala/designs.ps 

""" 

name = "Steiner Triple System on "+str(n)+" elements" 

if n%6 == 3: 

t = (n-3) // 6 

Z = list(range(2 * t + 1)) 

T = lambda x_y : x_y[0] + (2*t+1)*x_y[1] 

sts = [[(i,0),(i,1),(i,2)] for i in Z] + \ 

[[(i,k),(j,k),(((t+1)*(i+j)) % (2*t+1),(k+1)%3)] for k in range(3) for i in Z for j in Z if i != j] 

elif n%6 == 1: 

t = (n-1) // 6 

N = list(range(2 * t)) 

T = lambda x_y : x_y[0]+x_y[1]*t*2 if x_y != (-1,-1) else n-1 

L1 = lambda i,j : (i+j) % ((n-1)//3) 

L = lambda i,j : L1(i,j)//2 if L1(i,j)%2 == 0 else t+(L1(i,j)-1)//2 

sts = [[(i,0),(i,1),(i,2)] for i in range(t)] + \ 

[[(-1,-1),(i,k),(i-t,(k+1) % 3)] for i in range(t,2*t) for k in [0,1,2]] + \ 

[[(i,k),(j,k),(L(i,j),(k+1) % 3)] for k in [0,1,2] for i in N for j in N if i < j] 

else: 

raise EmptySetError("Steiner triple systems only exist for n = 1 mod 6 or n = 3 mod 6") 

# apply T and remove duplicates 

sts = set(frozenset(T(xx) for xx in x) for x in sts) 

return BalancedIncompleteBlockDesign(n, sts, name=name,check=False) 

def BIBD_from_TD(v,k,existence=False): 

r""" 

Return a BIBD through TD-based constructions. 

INPUT: 

- ``v,k`` (integers) -- computes a `(v,k,1)`-BIBD. 

- ``existence`` (boolean) -- instead of building the design, return: 

- ``True`` -- meaning that Sage knows how to build the design 

- ``Unknown`` -- meaning that Sage does not know how to build the 

design, but that the design may exist (see :mod:`sage.misc.unknown`). 

- ``False`` -- meaning that the design does not exist. 

This method implements three constructions: 

- If there exists a `TD(k,v)` and a `(v,k,1)`-BIBD then there exists a 

`(kv,k,1)`-BIBD. 

The BIBD is obtained from all blocks of the `TD`, and from the blocks of 

the `(v,k,1)`-BIBDs defined over the `k` groups of the `TD`. 

- If there exists a `TD(k,v)` and a `(v+1,k,1)`-BIBD then there exists a 

`(kv+1,k,1)`-BIBD. 

The BIBD is obtained from all blocks of the `TD`, and from the blocks of 

the `(v+1,k,1)`-BIBDs defined over the sets `V_1\cup \infty,\dots,V_k\cup 

\infty` where the `V_1,\dots,V_k` are the groups of the TD. 

- If there exists a `TD(k,v)` and a `(v+k,k,1)`-BIBD then there exists a 

`(kv+k,k,1)`-BIBD. 

The BIBD is obtained from all blocks of the `TD`, and from the blocks of 

the `(v+k,k,1)`-BIBDs defined over the sets `V_1\cup 

\{\infty_1,\dots,\infty_k\},\dots,V_k\cup \{\infty_1,\dots,\infty_k\}` 

where the `V_1,\dots,V_k` are the groups of the TD. By making sure that 

all copies of the `(v+k,k,1)`-BIBD contain the block 

`\{\infty_1,\dots,\infty_k\}`, the result is also a BIBD. 

These constructions can be found in 

`<http://www.argilo.net/files/bibd.pdf>`_. 

EXAMPLES: 

First construction:: 

sage: from sage.combinat.designs.bibd import BIBD_from_TD 

sage: BIBD_from_TD(25,5,existence=True) 

True 

sage: _ = BlockDesign(25,BIBD_from_TD(25,5)) 

Second construction:: 

sage: from sage.combinat.designs.bibd import BIBD_from_TD 

sage: BIBD_from_TD(21,5,existence=True) 

True 

sage: _ = BlockDesign(21,BIBD_from_TD(21,5)) 

Third construction:: 

sage: from sage.combinat.designs.bibd import BIBD_from_TD 

sage: BIBD_from_TD(85,5,existence=True) 

True 

sage: _ = BlockDesign(85,BIBD_from_TD(85,5)) 

No idea:: 

sage: from sage.combinat.designs.bibd import BIBD_from_TD 

sage: BIBD_from_TD(20,5,existence=True) 

Unknown 

sage: BIBD_from_TD(20,5) 

Traceback (most recent call last): 

... 

NotImplementedError: I do not know how to build a (20,5,1)-BIBD! 

""" 

# First construction 

if (v%k == 0 and 

balanced_incomplete_block_design(v//k,k,existence=True) and 

transversal_design(k,v//k,existence=True)): 

if existence: 

return True 

v = v//k 

BIBDvk = balanced_incomplete_block_design(v,k)._blocks 

TDkv = transversal_design(k,v,check=False) 

BIBD = TDkv._blocks 

for i in range(k): 

BIBD.extend([[x+i*v for x in B] for B in BIBDvk]) 

# Second construction 

elif ((v-1)%k == 0 and 

balanced_incomplete_block_design((v-1)//k+1,k,existence=True) and 

transversal_design(k,(v-1)//k,existence=True)): 

if existence: 

return True 

v = (v-1)//k 

BIBDv1k = balanced_incomplete_block_design(v+1,k)._blocks 

TDkv = transversal_design(k,v,check=False)._blocks 

inf = v*k 

BIBD = TDkv 

for i in range(k): 

BIBD.extend([[inf if x == v else x+i*v for x in B] for B in BIBDv1k]) 

# Third construction 

elif ((v-k)%k == 0 and 

balanced_incomplete_block_design((v-k)//k+k,k,existence=True) and 

transversal_design(k,(v-k)//k,existence=True)): 

if existence: 

return True 

v = (v-k)//k 

BIBDvpkk = balanced_incomplete_block_design(v+k,k) 

TDkv = transversal_design(k,v,check=False)._blocks 

inf = v*k 

BIBD = TDkv 

# makes sure that [v,...,v+k-1] is a block of BIBDvpkk. Then, we remove it. 

BIBDvpkk = _relabel_bibd(BIBDvpkk,v+k) 

BIBDvpkk = [B for B in BIBDvpkk if min(B) < v] 

for i in range(k): 

BIBD.extend([[(x-v)+inf if x >= v else x+i*v for x in B] for B in BIBDvpkk]) 

BIBD.append(list(range(k * v, v * k + k))) 

# No idea ... 

else: 

if existence: 

return Unknown 

else: 

raise NotImplementedError("I do not know how to build a ({},{},1)-BIBD!".format(v,k)) 

return BIBD 

def BIBD_from_difference_family(G, D, lambd=None, check=True): 

r""" 

Return the BIBD associated to the difference family ``D`` on the group ``G``. 

Let `G` be a group. A `(G,k,\lambda)`-*difference family* is a family `B = 

\{B_1,B_2,\ldots,B_b\}` of `k`-subsets of `G` such that for each element of 

`G \backslash \{0\}` there exists exactly `\lambda` pairs of elements 

`(x,y)`, `x` and `y` belonging to the same block, such that `x - y = g` (or 

x y^{-1} = g` in multiplicative notation). 

If `\{B_1, B_2, \ldots, B_b\}` is a `(G,k,\lambda)`-difference family then 

its set of translates `\{B_i \cdot g; i \in \{1,\ldots,b\}, g \in G\}` is a 

`(v,k,\lambda)`-BIBD where `v` is the cardinality of `G`. 

INPUT: 

- ``G`` - a finite additive Abelian group 

- ``D`` - a difference family on ``G`` (short blocks are allowed). 

- ``lambd`` - the `\lambda` parameter (optional, only used if ``check`` is 

``True``) 

- ``check`` - whether or not we check the output (default: ``True``) 

EXAMPLES:: 

sage: G = Zmod(21) 

sage: D = [[0,1,4,14,16]] 

sage: sorted(G(x-y) for x in D[0] for y in D[0] if x != y) 

[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] 

sage: from sage.combinat.designs.bibd import BIBD_from_difference_family 

sage: BIBD_from_difference_family(G, D) 

[[0, 1, 4, 14, 16], 

[1, 2, 5, 15, 17], 

[2, 3, 6, 16, 18], 

[3, 4, 7, 17, 19], 

[4, 5, 8, 18, 20], 

[5, 6, 9, 19, 0], 

[6, 7, 10, 20, 1], 

[7, 8, 11, 0, 2], 

[8, 9, 12, 1, 3], 

[9, 10, 13, 2, 4], 

[10, 11, 14, 3, 5], 

[11, 12, 15, 4, 6], 

[12, 13, 16, 5, 7], 

[13, 14, 17, 6, 8], 

[14, 15, 18, 7, 9], 

[15, 16, 19, 8, 10], 

[16, 17, 20, 9, 11], 

[17, 18, 0, 10, 12], 

[18, 19, 1, 11, 13], 

[19, 20, 2, 12, 14], 

[20, 0, 3, 13, 15]] 

""" 

from .difference_family import group_law, block_stabilizer 

identity, mul, inv = group_law(G) 

bibd = [] 

Gset = set(G) 

p_to_i = {g:i for i,g in enumerate(Gset)} 

for b in D: 

b = [G(_) for _ in b] 

S = block_stabilizer(G,b) 

GG = Gset.copy() 

while GG: 

g = GG.pop() 

if S: GG.difference_update(mul(s,g) for s in S) 

bibd.append([p_to_i[mul(i,g)] for i in b]) 

if check: 

if lambd is None: 

k = len(bibd[0]) 

v = G.cardinality() 

lambd = (len(bibd) * k * (k-1)) // (v * (v-1)) 

assert is_pairwise_balanced_design(bibd, G.cardinality(), [len(D[0])], lambd=lambd) 

return bibd 

################ 

# (v,4,1)-BIBD # 

################ 

def v_4_1_BIBD(v, check=True): 

r""" 

Return a `(v,4,1)`-BIBD. 

A `(v,4,1)`-BIBD is an edge-decomposition of the complete graph `K_v` into 

copies of `K_4`. For more information, see 

:func:`balanced_incomplete_block_design`. It exists if and only if `v\equiv 1,4 

\pmod {12}`. 

See page 167 of [Stinson2004]_ for the construction details. 

.. SEEALSO:: 

* :func:`balanced_incomplete_block_design` 

INPUT: 

- ``v`` (integer) -- number of points. 

- ``check`` (boolean) -- whether to check that output is correct before 

returning it. As this is expected to be useless (but we are cautious 

guys), you may want to disable it whenever you want speed. Set to ``True`` 

by default. 

EXAMPLES:: 

sage: from sage.combinat.designs.bibd import v_4_1_BIBD # long time 

sage: for n in range(13,100): # long time 

....: if n%12 in [1,4]: # long time 

....: _ = v_4_1_BIBD(n, check = True) # long time 

TESTS: 

Check that the `(25,4)` and `(37,4)`-difference family are available:: 

sage: assert designs.difference_family(25,4,existence=True) 

sage: _ = designs.difference_family(25,4) 

sage: assert designs.difference_family(37,4,existence=True) 

sage: _ = designs.difference_family(37,4) 

Check some larger `(v,4,1)`-BIBD (see :trac:`17557`):: 

sage: for v in range(400): # long time 

....: if v%12 in [1,4]: # long time 

....: _ = designs.balanced_incomplete_block_design(v,4) # long time 

""" 

k = 4 

if v == 0: 

return [] 

if v <= 12 or v%12 not in [1,4]: 

raise EmptySetError("A K_4-decomposition of K_v exists iif v=2,4 mod 12, v>12 or v==0") 

# Step 1. Base cases. 

if v == 13: 

# note: this construction can also be obtained from difference_family 

from .block_design import projective_plane 

return projective_plane(3)._blocks 

if v == 16: 

from .block_design import AffineGeometryDesign 

from sage.rings.finite_rings.finite_field_constructor import FiniteField 

return AffineGeometryDesign(2,1,FiniteField(4,'x'))._blocks 

if v == 25 or v == 37: 

from .difference_family import difference_family 

G,D = difference_family(v,4) 

return BIBD_from_difference_family(G,D,check=False) 

if v == 28: 

return [[0, 1, 23, 26], [0, 2, 10, 11], [0, 3, 16, 18], [0, 4, 15, 20], 

[0, 5, 8, 9], [0, 6, 22, 25], [0, 7, 14, 21], [0, 12, 17, 27], 

[0, 13, 19, 24], [1, 2, 24, 27], [1, 3, 11, 12], [1, 4, 17, 19], 

[1, 5, 14, 16], [1, 6, 9, 10], [1, 7, 20, 25], [1, 8, 15, 22], 

[1, 13, 18, 21], [2, 3, 21, 25], [2, 4, 12, 13], [2, 5, 18, 20], 

[2, 6, 15, 17], [2, 7, 19, 22], [2, 8, 14, 26], [2, 9, 16, 23], 

[3, 4, 22, 26], [3, 5, 7, 13], [3, 6, 14, 19], [3, 8, 20, 23], 

[3, 9, 15, 27], [3, 10, 17, 24], [4, 5, 23, 27], [4, 6, 7, 8], 

[4, 9, 14, 24], [4, 10, 16, 21], [4, 11, 18, 25], [5, 6, 21, 24], 

[5, 10, 15, 25], [5, 11, 17, 22], [5, 12, 19, 26], [6, 11, 16, 26], 

[6, 12, 18, 23], [6, 13, 20, 27], [7, 9, 17, 18], [7, 10, 26, 27], 

[7, 11, 23, 24], [7, 12, 15, 16], [8, 10, 18, 19], [8, 11, 21, 27], 

[8, 12, 24, 25], [8, 13, 16, 17], [9, 11, 19, 20], [9, 12, 21, 22], 

[9, 13, 25, 26], [10, 12, 14, 20], [10, 13, 22, 23], [11, 13, 14, 15], 

[14, 17, 23, 25], [14, 18, 22, 27], [15, 18, 24, 26], [15, 19, 21, 23], 

[16, 19, 25, 27], [16, 20, 22, 24], [17, 20, 21, 26]] 

# Step 2 : this is function PBD_4_5_8_9_12 

PBD = PBD_4_5_8_9_12((v-1)//(k-1),check=False) 

# Step 3 : Theorem 7.20 

bibd = BIBD_from_PBD(PBD,v,k,check=False) 

if check: 

assert is_pairwise_balanced_design(bibd,v,[k]) 

return bibd 

def BIBD_from_PBD(PBD,v,k,check=True,base_cases={}): 

r""" 

Return a `(v,k,1)`-BIBD from a `(r,K)`-PBD where `r=(v-1)/(k-1)`. 

This is Theorem 7.20 from [Stinson2004]_. 

INPUT: 

- ``v,k`` -- integers. 

- ``PBD`` -- A PBD on `r=(v-1)/(k-1)` points, such that for any block of 

``PBD`` of size `s` there must exist a `((k-1)s+1,k,1)`-BIBD. 

- ``check`` (boolean) -- whether to check that output is correct before 

returning it. As this is expected to be useless (but we are cautious 

guys), you may want to disable it whenever you want speed. Set to ``True`` 

by default. 

- ``base_cases`` -- caching system, for internal use. 

EXAMPLES:: 

sage: from sage.combinat.designs.bibd import PBD_4_5_8_9_12 

sage: from sage.combinat.designs.bibd import BIBD_from_PBD 

sage: from sage.combinat.designs.bibd import is_pairwise_balanced_design 

sage: PBD = PBD_4_5_8_9_12(17) 

sage: bibd = is_pairwise_balanced_design(BIBD_from_PBD(PBD,52,4),52,[4]) 

""" 

r = (v-1) // (k-1) 

bibd = [] 

for X in PBD: 

n = len(X) 

N = (k-1)*n+1 

if not (n,k) in base_cases: 

base_cases[n,k] = _relabel_bibd(balanced_incomplete_block_design(N,k), N) 

for XX in base_cases[n,k]: 

if N-1 in XX: 

continue 

bibd.append([X[x//(k-1)] + (x%(k-1))*r for x in XX]) 

for x in range(r): 

bibd.append([x+i*r for i in range(k-1)]+[v-1]) 

if check: 

assert is_pairwise_balanced_design(bibd,v,[k]) 

return bibd 

def _relabel_bibd(B,n,p=None): 

r""" 

Relabels the BIBD on `n` points and blocks of size k such that 

`\{0,...,k-2,n-1\},\{k-1,...,2k-3,n-1\},...,\{n-k,...,n-2,n-1\}` are blocks 

of the BIBD. 

INPUT: 

- ``B`` -- a list of blocks. 

- ``n`` (integer) -- number of points. 

- ``p`` (optional) -- the point that will be labeled with n-1. 

EXAMPLES:: 

sage: designs.balanced_incomplete_block_design(40,4).blocks() # indirect doctest 

[[0, 1, 2, 12], [0, 3, 6, 9], [0, 4, 8, 10], 

[0, 5, 7, 11], [0, 13, 26, 39], [0, 14, 25, 28], 

[0, 15, 27, 38], [0, 16, 22, 32], [0, 17, 23, 34], 

... 

""" 

if p is None: 

p = n-1 

found = 0 

last = n-1 

d = {} 

for X in B: 

if last in X: 

for x in X: 

if x == last: 

continue 

d[x] = found 

found += 1 

if found == n-1: 

break 

d[p] = n-1 

return [[d[x] for x in X] for X in B] 

def PBD_4_5_8_9_12(v, check=True): 

""" 

Return a `(v,\{4,5,8,9,12\})`-PBD on `v` elements. 

A `(v,\{4,5,8,9,12\})`-PBD exists if and only if `v\equiv 0,1 \pmod 4`. The 

construction implemented here appears page 168 in [Stinson2004]_. 

INPUT: 

- ``v`` -- an integer congruent to `0` or `1` modulo `4`. 

- ``check`` (boolean) -- whether to check that output is correct before 

returning it. As this is expected to be useless (but we are cautious 

guys), you may want to disable it whenever you want speed. Set to ``True`` 

by default. 

EXAMPLES:: 

sage: designs.balanced_incomplete_block_design(40,4).blocks() # indirect doctest 

[[0, 1, 2, 12], [0, 3, 6, 9], [0, 4, 8, 10], 

[0, 5, 7, 11], [0, 13, 26, 39], [0, 14, 25, 28], 

[0, 15, 27, 38], [0, 16, 22, 32], [0, 17, 23, 34], 

... 

Check that :trac:`16476` is fixed:: 

sage: from sage.combinat.designs.bibd import PBD_4_5_8_9_12 

sage: for v in (0,1,4,5,8,9,12,13,16,17,20,21,24,25): 

....: _ = PBD_4_5_8_9_12(v) 

""" 

if not v%4 in [0,1]: 

raise ValueError 

if v <= 1: 

PBD = [] 

elif v <= 12: 

PBD = [list(range(v))] 

elif v == 13 or v == 28: 

PBD = v_4_1_BIBD(v, check=False) 

elif v == 29: 

TD47 = transversal_design(4,7)._blocks 

four_more_sets = [[28]+[i*7+j for j in range(7)] for i in range(4)] 

PBD = TD47 + four_more_sets 

elif v == 41: 

TD59 = transversal_design(5,9) 

PBD = ([[x for x in X if x<41] for X in TD59]