pjt33

#!/usr/bin/pypy3

from functools import lru_cache

# Limit memory usage to 8GB
import resource
resource.setrlimit(resource.RLIMIT_AS, (8 * 1024**3, 8 * 1024**3))


"""

pjt33

#!/usr/bin/pypy3

"""
Ethan Dlugie asked [1] "How do we know there are no more Deligne-Mostow/Thurston lattices?"
than the 94 enumerated by Thurston [3]. This program combines case analysis with exhaustive search
to show that the list is in fact complete.

There is a subtle difference between Thurston's claim (Theorem 0.2) and Mostow's ΣINT condition [2]:
Mostow only allows one set of equal angles to give half-integer reciprocals (condition (c) below),
whereas Thurston allows multiple sets. This code uses Thurston's condition, as more straightforward,

pjt33

One solution per line. Format
n ((a_1, b_1), (a_2, b_2), ..., (a_n, b_n))

1 ((0, 0))

2 ((0, 1), (1, 0))

3 ((0, 1), (1, 1), (2, 0))
3 ((0, 2), (1, 0), (1, 1))

pjt33

R.<n> = PolynomialRing(QQ)

def p(x):
	return x^5 + n^2*x^4 - (2*n^3 + 6*n^2 + 10*n + 10)*x^3 + (n^4 + 5*n^3 + 11*n^2 + 15*n + 5)*x^2 + (n^3 + 4*n^2 + 10*n + 10)*x + 1

def s(x):
	return (n + 2 + n*x - x*x) / (1 + 2*x + n*x)

R2.<xtmp> = PolynomialRing(FractionField(R))
R3.<x1> = R2.quotient(p(xtmp))

pjt33

"""
Identify and count (by size) the subsets with odd beta_n using MacMahon's inclusion-exclusion formula for beta_n
as presented in theorem 2.3 of https://arxiv.org/pdf/1710.11033.pdf
"""

def bit_subsets_of(bitmask):
	# Nice bit hack from Thomas Beuman
	current = bitmask
	yield current
	while current:

pjt33

from collections import Counter
from math import gcd


def phi(n):
	# Not exactly over-optimised...
	rv = 1
	for p in range(2, n + 1):
		if n == 1:
			break

pjt33

def extend_closure(omega, included, extension, excluded):
	# We start with a closure which we want to extend
	closure = set(included)

	# Use a queue
	q = set([extension])
	while q:
		x = next(iter(q))
		q.remove(x)
		closure.add(x)

pjt33

from collections import defaultdict
from itertools import combinations
from sage.rings.polynomial.real_roots import real_roots

R.<alpha, x, y> = PolynomialRing(QQ, order='lex')
S.<z> = PolynomialRing(QQ)

def P(n, p, q):
    return ((1 - x^p) + alpha * (x^n - x^q)) // (1 - x^gcd(p, abs(n-q)))

pjt33

#!/usr/bin/pypy3

from heapq import heappush, heappop

def mod_inv(a: int, m: int) -> int:
	lastremainder, remainder = a % m, m
	x, lastx, y, lasty = 0, 1, 1, 0
	while remainder:
		lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
		x, lastx = lastx - quotient*x, x

pjt33

using System;
using System.Collections.Generic;
using System.Linq;
using BITMASK = System.UInt64;

namespace Sandbox
{
    using STATE = System.ValueTuple<BITMASK, System.Byte, System.Byte>;

    // One way to formalise the snake constraint is to treat internal vertices as deleted.