A script using the Z3 constraint solver to find projective planes.

projective_planes.py

#!/usr/bin/env python3

"""
Uses the Z3 constraint solver to find projective planes.

Usage:
- Create a venv. Install Z3 with Python extensions
- Ipython:

      %run projective_planes.py
      run(4)
"""

import math
import pprint
from z3 import *

def generate_pp_solver(size):
    """Generates a solver to find a projective plane of size `size`.
    """

    s = Solver()

    n_points = size * size + size + 1
    incidence = [
            BitVec('incidence_{}'.format(i), n_points) for i in range(n_points)
    ]

    # Rule 1: weight of each row and column is `size + 1`.
    n_bits = math.ceil(math.log2(n_points))
    for i in range(n_points):
        s.add(sum(ZeroExt(n_bits - 1, Extract(j, j, incidence[i])) for j in range(n_points)) == size + 1)
        s.add(sum(ZeroExt(n_bits - 1, Extract(i, i, incidence[j])) for j in range(n_points)) == size + 1)

    # Rule 2: each pair of rows has exactly one symbol in common.
    # In other words, their logical AND is a power of two.
    for i in range(n_points):
        for j in range(i+1, n_points):
            common = incidence[i] & incidence[j]
            is_power_of_two = And(common != 0, common & (common - 1) == 0)
            s.add(is_power_of_two)

    return s, incidence

def run(size):
    s, incidence = generate_pp_solver(size)
    if s.check() == sat:
        print("sat")
        m = s.model()
        result = [ m[i] for i in incidence ]
        print("result =")
        for row in result:
            print("  {row:0{width}b}".format(row=row.as_long(), width=len(incidence)))
    else:
        print("unsat")