# S1: check Macdonald's product formula for the Poincare polynomial with parameters,
#   W(t) = sum_w t_w = prod_{a>0} (1 - t_a t^{a}) / (1 - t^{a}),  t^{a} = prod_i t_i^{c_i(a)} (simple-(co)root coords),
# on the (co)root systems A3, B3, C3, D4, G2, B2, F4, with independent parameters per W-class,
# at random rational points (exact Fractions), including negative parameters.  Standard library only.
from fractions import Fraction as Fr
import random

def cartan(t, n):
    A = [[2 if i == j else 0 for j in range(n)] for i in range(n)]
    if t in "ABCD":
        for i in range(n - 1): A[i][i + 1] = A[i + 1][i] = -1
        if t == "B": A[n - 1][n - 2] = -2   # convention A[i][j] = <a_i^vee, a_j>
        if t == "C": A[n - 2][n - 1] = -2
        if t == "D":
            A[n - 2][n - 1] = A[n - 1][n - 2] = 0; A[n - 3][n - 1] = A[n - 1][n - 3] = -1
    if t == "G": A = [[2, -1], [-3, 2]]
    if t == "F": A = [[2, -1, 0, 0], [-1, 2, -2, 0], [0, -1, 2, -1], [0, 0, -1, 2]]
    return A

def system(A):
    n = len(A)
    # vectors in simple-root coordinates; s_i(v) = v - <v, a_i^vee> a_i, <a_j, a_i^vee> = A[i][j]
    def pair(v, i): return sum(v[j] * A[i][j] for j in range(n))
    def s(i, v):
        p = pair(v, i); w = list(v); w[i] -= p; return tuple(w)
    simple = [tuple(1 if j == i else 0 for j in range(n)) for i in range(n)]
    pos = set(simple); fr = list(simple)
    while fr:
        new = []
        for v in fr:
            for i in range(n):
                w = s(i, v)
                if all(x >= 0 for x in w) and w not in pos: pos.add(w); new.append(w)
        fr = new
    # W elements by BFS on images of a regular vector
    reg = tuple(Fr(1, 1 + 3 * i) for i in range(n))
    def act(word, v):
        for i in reversed(word): v = s(i, v)
        return v
    elems = {act((), reg): ()}; fr = [()]
    while fr:
        new = []
        for w in fr:
            for i in range(n):
                w2 = (i,) + w; k = act(w2, reg)
                if k not in elems: elems[k] = w2; new.append(w2)
        fr = new
    # W-classes of simple roots
    cls = list(range(n))
    for i in range(n):
        for j in range(i):
            if any(act(w, simple[i]) == simple[j] for w in elems.values()):
                cls[i] = cls[j]; break
    # class of each positive root
    pcls = {}
    for p in pos:
        for i in range(n):
            if any(act(w, simple[i]) == p for w in elems.values()):
                pcls[p] = cls[i]; break
    return n, sorted(pos), list(elems.values()), cls, pcls

random.seed(1)
for (t, n) in [("A", 3), ("B", 2), ("B", 3), ("C", 3), ("D", 4), ("G", 2), ("F", 4)]:
    n, pos, W, cls, pcls = system(cartan(t, n))
    ok = True
    for trial in range(4):
        par = {c: Fr(random.choice([-1, 1]) * random.randint(1, 9), random.randint(10, 40)) for c in set(cls)}
        tp = [par[cls[i]] for i in range(n)]
        lhs = sum((lambda w: __import__("functools").reduce(lambda a, b: a * b, [tp[i] for i in w], Fr(1)))(w) for w in W)
        rhs = Fr(1)
        for p in pos:
            tpow = Fr(1)
            for i in range(n): tpow *= tp[i] ** p[i]
            rhs *= (1 - par[pcls[p]] * tpow) / (1 - tpow)
        ok = ok and lhs == rhs
    print("%s%d: |W| = %d, classes of simple reflections %s, product formula exact on 4 random parameter points (some negative): %s" % (t, n, len(W), cls, ok))