"""Root systems: port of rootSystem.pm.

Types A, B, C, D, E (rank 6, 7, 8), F, F*, G, G* (F* and G* are F4 and
G2 with coordinates dual to the Bourbaki ones).  All matrices are lists
of rows.  Entries are Perl-like values (see :mod:`perlval`): integers,
fractions, or strings produced by :func:`fix`, exactly as in the Perl
original, so that printed output matches it.

``labels`` is '' (Bourbaki labelling), 'gap' (the Gap labelling: the
reverse order in types B, C, D, G) or a list giving a permutation of
1..n.
"""

import re

from .perlval import (PFrac, frac, pstr, pnum, padd, pmul, pneg, pcmp,
                      fix_str, n_pow)

test = 0


def f(x):
    """Perl f(): Math::Fraction::frac of a Perl number."""
    return frac(x)


def _is_gap(labels):
    return isinstance(labels, str) and re.search("gap", labels, re.I)


def _is_list(labels):
    return isinstance(labels, (list, tuple))


def zeroMatrix(numberOfRows, numberOfColumns):
    """An all-zero matrix."""
    return [[0] * numberOfColumns for _ in range(numberOfRows)]


def cartanMatrix(type, rank, labels=""):
    """Cartan matrix of the given type and rank (permuted by labels)."""
    if _is_gap(labels):
        # The Perl code sets $labels=($rank) here (a bug which is never
        # reached: callers resolve 'gap' before calling).
        labels = getGapLabels(type, rank)
    matrix = zeroMatrix(rank, rank)
    for i in range(rank):
        matrix[i][i] = 2
    start = 2 if type == "E" else 0
    end = rank - 3 if type == "D" else rank - 2
    for i in range(start, end + 1):
        matrix[i][i + 1] = -1
        matrix[i + 1][i] = -1
    if type == "B" and rank >= 2:
        matrix[rank - 2][rank - 1] = -2
    elif type == "C" and rank >= 2:
        matrix[rank - 1][rank - 2] = -2
    elif type == "D" and rank >= 3:
        matrix[rank - 3][rank - 1] = -1
        matrix[rank - 1][rank - 3] = -1
    elif type == "F":
        matrix[1][2] = -2
    elif type == "F*":
        matrix[2][1] = -2
    elif type == "G":
        matrix[1][0] = -3
    elif type == "G*":
        matrix[0][1] = -3
    elif type == "E":
        matrix[0][2] = -1
        matrix[2][0] = -1
        matrix[1][3] = -1
        matrix[3][1] = -1
    if labels:
        P = permutationMatrix(labels)
        matrix = matrixMult(matrixMult(P, matrix), transpose(P))
    return matrix


def showMatrix(m, out):
    """Print a matrix, one comma separated row per line."""
    for row in m:
        out.append(",".join(pstr(x) for x in fix(row)) + "\n")


def fix(row):
    """Strip a trailing "/1" from every entry (Perl s/\\/1$//)."""
    return [fix_str(x) for x in row]


def fundamentalWeights(type, rank, labels=""):
    """Fundamental weights in standard coordinates, one per row."""
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    fw = []
    if type == "A":
        n = rank + 1
        for i in range(1, rank + 1):
            fw.append([frac(n - i, n)] * i + [frac(-i, n)] * (n - i))
    elif type == "B":
        for i in range(1, rank):
            fw.append([1] * i + [0] * (rank - i))
        fw.append([f(1 / 2)] * rank)
    elif type == "C":
        for i in range(1, rank + 1):
            fw.append([1] * i + [0] * (rank - i))
    elif type == "D":
        for i in range(1, rank - 1):
            fw.append([1] * i + [0] * (rank - i))
        fw.append([f(1 / 2)] * (rank - 1) + [f(-1 / 2)])
        fw.append([f(1 / 2)] * rank)
    elif type == "F":
        fw = [[1, 1, 0, 0], [2, 1, 1, 0],
              [f(3 / 2), f(1 / 2), f(1 / 2), f(1 / 2)], [1, 0, 0, 0]]
    elif type == "F*":
        fw = [[1, 1, 0, 0], [2, 1, 1, 0], [3, 1, 1, 1], [2, 0, 0, 0]]
    elif type == "G":
        fw = [[0, -1, 1], [-1, -1, 2]]
    elif type == "G*":
        fw = [[0, -1, 1], [f(-1 / 3), f(-1 / 3), f(2 / 3)]]
    elif type == "E":
        h = f(1 / 2)
        mh = f(-1 / 2)
        if rank == 6:
            fw = [[0, 0, 0, 0, 0, f(-2 / 3), f(-2 / 3), f(2 / 3)],
                  [h, h, h, h, h, mh, mh, h],
                  [mh, h, h, h, h, f(-5 / 6), f(-5 / 6), f(5 / 6)],
                  [0, 0, 1, 1, 1, -1, -1, 1],
                  [0, 0, 0, 1, 1, f(-2 / 3), f(-2 / 3), f(2 / 3)],
                  [0, 0, 0, 0, 1, f(-1 / 3), f(-1 / 3), f(1 / 3)]]
        elif rank == 7:
            fw = [[0, 0, 0, 0, 0, 0, -1, 1],
                  [h, h, h, h, h, h, -1, 1],
                  [mh, h, h, h, h, h, f(-3 / 2), f(3 / 2)],
                  [0, 0, 1, 1, 1, 1, -2, 2],
                  [0, 0, 0, 1, 1, 1, f(-3 / 2), f(3 / 2)],
                  [0, 0, 0, 0, 1, 1, -1, 1],
                  [0, 0, 0, 0, 0, 1, mh, h]]
        elif rank == 8:
            fw = [[0, 0, 0, 0, 0, 0, 0, 2],
                  [h, h, h, h, h, h, h, f(5 / 2)],
                  [mh, h, h, h, h, h, h, f(7 / 2)],
                  [0, 0, 1, 1, 1, 1, 1, 5],
                  [0, 0, 0, 1, 1, 1, 1, 4],
                  [0, 0, 0, 0, 1, 1, 1, 3],
                  [0, 0, 0, 0, 0, 1, 1, 2],
                  [0, 0, 0, 0, 0, 0, 1, 1]]
    if labels:
        P = permutationMatrix(labels)
        fw = matrixMult(P, fw)
    return fw


def rho(type, rank, labels=""):
    """One half the sum of the positive roots (sum of fundamental weights)."""
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    fws = fundamentalWeights(type, rank, labels)
    r = [0] * len(fws[0])
    for row in fws:
        for i in range(len(row)):
            r[i] = padd(r[i], row[i])
    return r


def permutationMatrix(permutation):
    """Permutation matrix P with P[j-1][i] = 1 for j = permutation[i]."""
    rank = len(permutation)
    P = zeroMatrix(rank, rank)
    for i, j in enumerate(permutation):
        P[int(pnum(j)) - 1][i] = 1
    return P


def fundamentalCoWeights(type, rank, labels=""):
    """Fundamental coweights = fundamental weights of the dual type."""
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    return fundamentalWeights(dual(type), rank, labels)


def simpleRoots(type, rank, labels=""):
    """Simple roots = Cartan matrix times fundamental weights."""
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    sr = matrixMult(cartanMatrix(type, rank), fundamentalWeights(type, rank))
    if labels:
        P = permutationMatrix(labels)
        sr = matrixMult(P, sr)
    return sr


def simpleCoRoots(type, rank, labels=""):
    """Simple coroots = simple roots of the dual type."""
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    return simpleRoots(dual(type), rank, labels)


def transpose(matrix):
    rows = len(matrix)
    cols = len(matrix[0])
    return [[matrix[j][i] for j in range(rows)] for i in range(cols)]


def matrixMult(A, B):
    """Matrix product; entries are passed through fix() like the Perl."""
    Arows, Acols = len(A), len(A[0])
    Brows, Bcols = len(B), len(B[0])
    if Acols != Brows:
        raise ValueError("IndexError: matrices don't match: %d != %d"
                         % (Acols, Brows))
    result = []
    for i in range(Arows):
        row = []
        for j in range(Bcols):
            acc = None
            for k in range(Acols):
                acc = padd(acc, pmul(A[i][k], B[k][j]))
            row.append(acc)
        result.append(fix(row))
    return result


def cleanUp(values):
    return [fix_str(x) for x in values]


def showAll(type, rank, labels="", out=None):
    """All the root data for a type and rank, as text (rootSystem.cgi)."""
    if out is None:
        out = []
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    p = out.append
    if "*" in type:
        typeNoStar = type.replace("*", "", 1)
        p("Coordinates for %s%s are dual to those of %s%s\n"
          % (type, rank, typeNoStar, rank))
    p("\nDynkin Diagram:\n\n")
    diagram(type, rank, labels, out)
    p("\n\nCartan Matrix:%s\n" % type)
    showMatrix(cartanMatrix(type, rank, labels), out)
    p("\nSimple Roots:\n")
    showMatrix(simpleRoots(type, rank, labels), out)
    p("\nFundamental Weights:\n")
    showMatrix(fundamentalWeights(type, rank, labels), out)
    p("\nOne half sum of the positive roots: ")
    # (The Perl printed "$type $rank $labels" here, a debugging leftover.)
    r = rho(type, rank, labels)
    p(",".join(pstr(x) for x in cleanUp(r)))
    p("\n\nOne half sum of the positive co-roots: ")
    p(",".join(pstr(x) for x in cleanUp(rho(dual(type), rank, labels))))
    p("\n\nSimple CoRoots:\n")
    showMatrix(simpleCoRoots(type, rank, labels), out)
    p("\nFundamental CoWeights:\n")
    fcw = fundamentalCoWeights(type, rank, labels)
    showMatrix(fundamentalCoWeights(type, rank, labels), out)
    order = commify(orderWeylGroup(type, rank))
    p("\nOrder of Weyl Group: %s\n" % order)
    pos = positiveRoots(type, rank, labels)
    p("\nNumber of Positive roots: %d" % len(pos))
    m = matrixMult(pos, transpose(fcw))
    rhoCheck = transpose([rho(dual(type), rank, labels)])
    p("\n\nPositive roots:\n")
    p("Standard Coordinates:Simple Root Coordinates:Height\n")
    for i in range(len(m)):
        alpha = ",".join(pstr(x) for x in pos[i])
        beta = ",".join(pstr(x) for x in m[i])
        height = matrixMult([pos[i]], rhoCheck)[0][0]
        p("%s\t %s\t   %s\n" % (alpha, beta, pstr(height)))
    return out


def orderWeylGroup(type, rank):
    """Order of the Weyl group."""
    if type == "A":
        return factorial(rank + 1)
    if re.search("B|C", type):
        return n_pow(2, rank) * factorial(rank)
    if type == "D":
        return n_pow(2, rank - 1) * factorial(rank)
    if type == "E":
        return {6: 51840, 7: 2903040, 8: 696729600}.get(rank)
    if type == "F":
        return 1152
    if type == "G":
        return 12
    return None


def factorial(n):
    return n * factorial(n - 1) if n > 1 else 1


def commify(n):
    s = pstr(n)
    while True:
        new = re.sub(r"^(-?\d+)(\d{3})", r"\1,\2", s)
        if new == s:
            return s
        s = new


def diagram(type, rank, labels="", out=None):
    """ASCII Dynkin diagram with the given node labels."""
    if out is None:
        out = []
    p = out.append
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    if not labels:
        labels = list(range(1, rank + 2))
    labels = [pstr(x) for x in labels]

    class _L(object):
        """Perl array indexing: negative indices count from the end,
        missing elements are empty."""
        def __getitem__(self, i):
            if i < 0:
                i += len(labels)
            return labels[i] if 0 <= i < len(labels) else ""
    L = _L()
    tail = {"A": "---", "B": "=>=", "C": "=<=", "D": ""}
    if type == "D":
        spaces = (rank - 3) * 4
        p(" " * spaces + L[rank - 2] + "\n")
        p(" " * spaces + "|\n")
        p(" " * spaces + "|\n")
    if re.search("A|B|C|D", type):
        for i in range(1, rank - 1):
            p(L[i - 1] + "---")
        if type != "D":
            p(L[rank - 2])
        p(tail.get(type, ""))
        p(L[rank - 1])
    elif type == "E":
        p("\n        %s\n\t|\n\t|\n%s---%s---%s---%s---%s"
          % (L[1], L[0], L[2], L[3], L[4], L[5]))
        for i in range(7, rank + 1):
            p("---" + L[i - 1])
    elif type == "F":
        p("%s---%s=>=%s---%s" % (L[0], L[1], L[2], L[3]))
    elif type == "F*":
        p("%s---%s=<=%s---%s" % (L[0], L[1], L[2], L[3]))
    elif type == "G":
        p("1=<=2")
    elif type == "G*":
        p("1=>=2")
    return out


def positiveRoots(type, rank, labels=""):
    """Positive roots in standard coordinates, sorted by height."""
    if _is_gap(labels):
        labels = getGapLabels(type, rank)
    sr = simpleRoots(type, rank, labels)
    size = len(sr[0]) - 1
    P = []
    if re.search("A|B|C|D|F", type):
        for i in range(0, size):
            for j in range(i + 1, size + 1):
                P.append(e(1, -1, i, j, size))
    if re.search("B|C|D|F", type):
        for i in range(0, size):
            for j in range(i + 1, size + 1):
                P.append(e(1, 1, i, j, size))
    if re.search("B|C|F", type):
        length = 2 if re.search("C|f", type) else 1
        for i in range(0, size + 1):
            v = [0] * (size + 1)
            v[i] = length
            P.append(v)
    if re.search("E", type):
        n = rank - 2
        m = 7 if rank == 8 else n
        tails = {6: [frac(-1 / 2), frac(-1 / 2), frac(1 / 2)],
                 7: [frac(-1 / 2), frac(1 / 2)],
                 8: [frac(1 / 2)]}
        for i in range(0, m + 1):
            for j in range(i + 1, m + 1):
                P.append(e(1, 1, i, j, size))
                P.append(e(-1, 1, i, j, size))
        if rank == 7:
            P.append([0, 0, 0, 0, 0, 0, -1, 1])
        for i in range(1, 2 ** n + 1):
            parity = 1
            v = [frac(1 / 2)] * n
            for j in range(0, n):
                if i & (2 ** j):
                    v[j] = frac(-1 / 2)
                    parity = -parity
            v.append(pmul(parity * (-1) ** rank, frac(1 / 2)))
            v.extend(tails[rank])
            P.append(v)
    if type == "G":
        P = [[1, -1, 0], [-2, 1, 1], [-1, 0, 1], [0, -1, 1], [1, -2, 1],
             [-1, -1, 2]]
    if type == "G*":
        P = [[1, -1, 0], [f(-2 / 3), f(1 / 3), f(1 / 3)], [-1, 0, 1],
             [0, -1, 1], [f(1 / 3), f(-2 / 3), f(1 / 3)],
             [f(-1 / 3), f(-1 / 3), f(2 / 3)]]
    if re.search("F", type):
        c = 1 if re.search("f", type) else frac(1 / 2)
        for i in range(1, 2 ** size + 1):
            v = [c] * (size + 1)
            for j in range(1, size + 1):
                v[j] = pneg(c) if (i & 2 ** (j - 1)) else c
            P.append(v)
    rhoCheck = transpose([rho(dual(type), rank, labels)])
    heights = [matrixMult([v], rhoCheck)[0][0] for v in P]
    order = sorted(range(len(P)),
                   key=lambda k: _Key(heights[k]))
    return [P[k] for k in order]


class _Key(object):
    """Sort key using Perl's numeric comparison (<=>)."""
    __slots__ = ("v",)

    def __init__(self, v):
        self.v = v

    def __lt__(self, other):
        return pcmp(self.v, other.v) < 0


def e(sign_i, sign_j, i, j, size):
    v = [0] * (size + 1)
    v[i] = sign_i
    v[j] = sign_j
    return v


def getGapLabels(type, rank):
    """The Gap labelling: n, n-1, ..., 1 in types B, C, D, G; else ''."""
    if not re.search("B|C|D|G", type):
        return ""
    return [rank - i for i in range(rank)]


def dual(type):
    return {"B": "C", "C": "B", "F": "F*", "F*": "F", "G": "G*",
            "G*": "G"}.get(type, type)