#!/usr/bin/env sage
"""\
Usage:
    fuchsia [options] <command> <args>...

Commands:
    reduce [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>
        find an epsilon form of the given matrix

    fuchsify [-x <name>] [-m <path>] [-t <path>] <matrix>
        find a transformation that will transform a given matrix
        into Fuchsian form

    normalize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>
        find a transformation that will transform a given Fuchsian
        matrix into normalized form

    factorize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>
        find a transformation that will make a given normalized
        matrix proportional to the infinitesimal parameter

    info [-x <name>] <matrix>
        show a description of a given matrix, listing singular
        points and residue eigenvalues (in verbose more)

    sort [-m <path>] [-t <path>] <matrix>
        find a block-triangular form of the given matrix

    transform [-x <name>] [-m <path>] <matrix> <transform>
        transform a given matrix using a given transformation

    changevar [-x <name>] [-y <name>] [-m <path>] <matrix> <expr>
        transform a given matrix by susbtituting free variable
        by a given expression

Options:
    -h          show this help message
    -f <fmt>    matrix file format: mtx or m (default: m)
    -l <path>   write log to this file
    -v          produce a more verbose log
    -P <path>   save profile report into this file
    -x <name>   use this name for the free variable (default: x)
    -y <name>   use this name for the new free variable (default: y)
    -e <name>   use this name for the infinitesimal parameter (default: eps)
    -m <path>   save the resulting matrix into this file
    -t <path>   save the resulting transformation into this file
    -s <number> use this random seed when applicable (default: 0)
    --use-maple speed up calculations by using Maple when possible

Arguments:
    <matrix>    read the input matrix from this file
    <transform> read the transformation matrix from this file
    <expr>      arbitrary expression
"""

__author__ = "<NAME>, <NAME>"
__author_email__ = "<EMAIL>"
__version__ = "17.06.28"

__all__ = [
    "FuchsianSystem",
    "NormalizeAssistant",

    "balance",
    "balance_transform",
    "block_triangular_form",
    "epsilon_form",
    "export_matrix_to_file",
    "factorize",
    "fuchsify",
    "fuchsify_off_diagonal_blocks",
    "import_matrix_from_file",
    "is_fuchsian",
    "matrix_c0",
    "matrix_c1",
    "matrix_complexity",
    "matrix_residue",
    "normalize",
    "reduce_diagonal_blocks",
    "setup_fuchsia",
    "simple",
    "simplify_by_factorization",
    "simplify_by_jordanification",
    "singularities",
    "transform"
]

from   collections import defaultdict
from   functools import wraps
from   itertools import combinations, permutations
from   random import Random
import inspect
import logging
import time

from   sage.all import *
from   sage.misc.parser import Parser
from   sage.libs.ecl import ecl_eval

def logcall(f):
    @wraps(f)
    def wrapper(*args, **kwargs):
        logger.enter(f.__name__)
        r = f(*args, **kwargs)
        logger.exit(f.__name__)
        return r
    return wrapper

class ElapsedTimeFormatter(logging.Formatter):
    def __init__(self, fmt=None, datefmt=None):
        super(ElapsedTimeFormatter, self).__init__(fmt, datefmt)
        self.start_time = time.time()
    def formatTime(self, record, datefmt=None):
        if datefmt is None:
            return "%7.1fs" % (record.created - self.start_time)
        return time.strftime(datefmt, time.localtime(record.created))

class FuchsiaLogger(object):
    def __init__(self):
        log_handler = logging.StreamHandler()
        log_handler.setFormatter(ElapsedTimeFormatter(
            "\033[32m[%(asctime)s]\033[0m %(message)s"
        ))
        logger = logging.getLogger('fuchsia')
        logger.addHandler(log_handler)
        logger.setLevel(logging.WARNING)

        self.addHandler = logger.addHandler
        self.error = logger.error
        self.isEnabledFor = logger.isEnabledFor
        self.logger = logger
        self.setLevel = logger.setLevel

        self.depth = 0

    def indent(self):
        return "   " * self.depth

    def is_verbose(self):
        return self.isEnabledFor(logging.INFO)

    def enter(self, name):
        self.info("-> %s" % name)
        self.depth += 1

    def exit(self, name):
        self.depth -= 1
        self.info("<- %s" % name)

    def debug(self, msg):
        self.logger.debug(self.indent() + msg)

    def info(self, msg):
        self.logger.info(self.indent() + msg)

if True:
    ecl_eval("(ext:set-limit 'ext:heap-size 0)")
    logger = FuchsiaLogger()
    USE_MAPLE = False

def setup_fuchsia(verbosity=0, use_maple=False):
    global USE_MAPLE
    USE_MAPLE = bool(use_maple)
    logger.setLevel(
        logging.WARNING if verbosity <= 0 else \
        logging.INFO if verbosity == 1 else
        logging.DEBUG
    )

class FuchsiaError(Exception):
    pass

def cross_product(v1, v2):
    m1, m2 = matrix(v1), matrix(v2)
    return m1.transpose() * m2

def dot_product(v1, v2):
    m1, m2 = matrix(v1), matrix(v2)
    sp = m1 * m2.transpose()
    return sp[0,0]

def partial_fraction(M, var):
    return M.apply_map(lambda ex: ex.partial_fraction(var))

#@logcall
def fuchsia_simplify(obj, x=None):
    if USE_MAPLE:
        def maple_simplify(ex):
            if hasattr(ex, "variables") and ((x is None) or (x in ex.variables())):
                res = maple.factor(maple.radnormal(ex))
                return parse(str(res))
            else:
                return ex
        if hasattr(obj, "apply_map"):
            return obj.apply_map(maple_simplify)
        else:
            return maple_simplify(obj)
    else:
        return obj.simplify_rational()

def fuchsia_solve(eqs, var):
    if USE_MAPLE:
        s = maple.solve(eqs, var)
        solutions = s.parent().get(s._name).strip('[]').split('],[')
        solutions = [s.split(',') for s in solutions if s != '']
        result = []
        for solution in solutions:
            r = []
            for s in solution:
                try:
                    expr = fuchsia_simplify(parse(s))
                    r.append(expr)
                except SyntaxError as error:
                    print("ERROR:  \n%s\n  %s\n" % (s, error))
                    continue
            result.append(r)
        return result
    else:
        return solve(eqs, var, solution_dict=True)

def change_variable(m, x, y, fy):
    mm = m.subs({x: fy}) * derivative(fy, y)
    return mm

def transform(M, x, T):
    """Given a system of differential equations dF/dx = M*F,
    and a transformation of base functions F = T*F', compute
    and return M', such that dF'/dx = M'*F'.

    Note: M' = inverse(T)*(M*T - dT/dx)
    """
    mm = T.inverse()*(M*T - derivative(T, x))
    return mm

def balance(P, x1, x2, x):
    assert P.is_square()
    assert (P*P - P).is_zero()
    assert not bool(x1 == x2)
    coP = identity_matrix(P.nrows()) - P
    if x1 == oo:
        b = coP - (x - x2)*P
    elif x2 == oo:
        b = coP - 1/(x - x1)*P
    else:
        b = coP + (x - x2)/(x - x1)*P
    return b

def balance_transform(M, P, x1, x2, x):
    """Same thing as transform(M, x, balance(P, x1, x2, x)), but faster."""
    assert P.is_square()
    #assert (P*P - P).is_zero()
    assert not bool(x1 == x2)
    coP = identity_matrix(P.nrows()) - P
    if x1 == oo:
        k = -(x - x2)
        d = -1
    elif x2 == oo:
        k = -1/(x - x1)
        d = 1/(x - x1)**2
    else:
        k = (x - x2)/(x - x1)
        d = (x2 - x1)/(x - x1)**2
    mm = (coP + 1/k*P)*M*(coP + k*P) - d/k*P
    return mm

#@logcall
def limit_fixed(expr, x, x0):
    #logger.debug("[%s -> %s] %s" % (x, x0, expr))
    res = None
    if USE_MAPLE:
        res = limit_fixed_maple(expr, x, x0)
    else:
        res = limit_fixed_maxima(expr, x, x0)
    #logger.debug("res = %s" % res)
    return res

def limit_fixed_maple(expr, x, x0):
    res = maple.limit(expr, **{str(x): x0})
    try:
        res = parse(str(res))
    except:
        res = NaN
    return res

def limit_fixed_maxima(expr, x, x0):
    """Return a limit of expr when x->lim.

    The standard 'limit()' function of SageMath does not allow
    you to specify the variable, only it's name as a keyword
    argument. If you have a variable, and not it's name, use
    this function instead.
    """
    l = maxima_calculus.sr_limit(expr, x, x0)
    return expr.parent()(l)

@logcall
def singularities(m, x):
    """Find values of x around which rational matrix M has
    a singularity; return a dictionary with {val: p} entries,
    where p is the Poincare rank of M at x=val.

    Example:
    >>> x = var("x")
    >>> s = singularities(matrix([[1/x, 0], [1/(x+1)**3, 1/(x+1)]]), x)
    >>> sorted(s.items())
    [(-1, 2), (0, 0), (+Infinity, 0)]

    >>> s = singularities(matrix([[x, 1, 1/x]]), x)
    >>> sorted(s.items())
    [(0, 0), (+Infinity, 2)]
    """
    m = fuchsia_simplify(m)
    result = {}
    for expr in m.list():
        points_expr = singularities_expr(expr, x)
        for x0, p in points_expr.items():
            if x0 in result:
                result[x0] = max(result[x0], p)
            else:
                result[x0] = p
    return result

def singularities_expr(expr, x):
    if USE_MAPLE:
        return singularities_expr_maple(expr, x)
    else:
        return singularities_expr_maxima(expr, x)

def singularities_expr_maple(expr, x):
    if bool(expr == 0):
        return {}

    result = {}
    sols = fuchsia_solve(1/expr, x)
    points = [x0 for x0 in sols[0]] if len(sols) > 0 else []
    for x0 in points:
        if x0 not in result:
            result[x0] = 0
        else:
            result[x0] += 1

    sols = fuchsia_solve((1/(expr.subs({x: 1/x})/x**2)).simplify_rational(), x)
    points = [x0 for x0 in sols[0]] if len(sols) > 0 else []
    for x0 in points:
        if x0 == 0:
            if oo not in result:
                result[oo] = 0
            else:
                result[oo] += 1
    return result

def singularities_expr_maxima(expr, x):
    if expr.is_zero():
        return {}

    result = {}
    eq = factor(1/expr)
    points, ps = solve(eq, x, solution_dict=True, multiplicities=True)
    for sol, p in zip(points, ps):
        if x not in sol:
            raise Exception("Maxima can't solve this equation: %s" % eq)
        x0 = expand(sol[x])
        if p >= 1:
            result[x0] = p-1

    eq = factor(1/(expr.subs({x: 1/x})/x**2))
    points, ps = solve(eq, x, solution_dict=True, multiplicities=True)
    for sol, p in zip(points, ps):
        if x not in sol:
            raise Exception("Maxima can't solve this equation: %s" % eq)
        if sol[x] == 0 and p >= 1:
            result[oo] = p-1
    return result

def matrix_taylor0(m, x, x0, exp):
    """Return the 0-th coefficient of Taylor expansion of
    a matrix M around a finite point x=point, assuming that
    M(x->point)~1/(x-point)**exp.

    Example:
    >>> x = var('x')
    >>> matrix_taylor0(matrix([[x/(x-1), 1/x, x, 1]]), x, 1, 1)
    [1 0 0 0]
    """
    return m.apply_map(lambda ex: limit_fixed(ex*(x-x0)**exp, x, x0))

def matrix_taylor1(M, x, point, exp):
    """Return the 1-th coefficient of Taylor expansion of
    a matrix M around a finite point x=point, assuming that
    M(x->point)~1/(x-point)**exp.

    Example:
    >>> x = var('x')
    >>> matrix_taylor1(matrix([[x/(x-1), 1/x, x, 1]]), x, 0, 1)
    [0 0 0 1]
    """
    def taylor1(e, x):
        l0 = limit_fixed(e, x, 0)
        l1 = limit_fixed((e - l0)/x, x, 0)
        return l1
    return matrix([
        [taylor1(e, x) for e in row]
        for row in M.subs({x: x+point})*x**exp
    ])

def matrix_c0(M, x, point, p):
    """Return the 0-th coefficient of M's expansion at x=point,
    assuming Poincare rank of M at that point is p. If point is
    +Infinity, return minus the coefficient at the highest power of x.

    Examples:
    >>> x = var("x")
    >>> m = matrix([[1/x, 1/x**2], [1, 1/(x-1)]])
    >>> matrix_c0(m, x, 0, 1)
    [0 1]
    [0 0]
    >>> matrix_c0(m, x, 1, 0)
    [0 0]
    [0 1]
    >>> matrix_c0(m, x, oo, 1)
    [ 0  0]
    [-1  0]
    >>> matrix_c0(m*x, x, oo, 2)
    [ 0  0]
    [-1  0]
    """
    if point == oo:
        return -matrix_taylor0(M.subs({x: 1/x}), x, 0, p-1)
    else:
        return matrix_taylor0(M, x, point, p+1)

def matrix_c1(M, x, point, p):
    """Return the 1-st coefficient of M's expansion at x=point,
    assuming Poincare rank of M at that point is p. If point is
    +Infinity, return minus the coefficient at the second-to-highest
    power of x.

    Examples:
    >>> x = var("x")
    >>> m = matrix([[1/x, 1/x**2], [1, 1/(x-1)]])
    >>> matrix_c1(m, x, 0, 1)
    [1 0]
    [0 0]
    >>> matrix_c1(m, x, oo, 1)
    [-1  0]
    [ 0 -1]
    """
    if point == oo:
        return -matrix_taylor1(M.subs({x: 1/x}), x, 0, p-1)
    else:
        return matrix_taylor1(M, x, point, p+1)

def matrix_residue(M, x, x0):
    """Return matrix residue of M at x=x0, assuming that M's
    Poincare rank at x=x0 is zero.

    Example:
    >>> x = var("x")
    >>> m = matrix([[1/x, 2/x], [3/x, 4/x]])
    >>> matrix_residue(m, x, 0)
    [1 2]
    [3 4]
    >>> matrix_residue(m, x, oo)
    [-1 -2]
    [-3 -4]
    """
    if M._cache is None:
        M._cache = {}
    key = "matrix_residue_%s_%s" % (x, x0)
    if key in M._cache:
        return M._cache[key]

    m0 = matrix_c0(M, x, x0, 0)
    M._cache[key] = m0
    return m0

def matrix_is_nilpotent(M):
    """Return True if M is always nilpotent, False otherwise.

    Examples:
    >>> a, b, c = var('a b c')
    >>> matrix_is_nilpotent(matrix([[0,a,b],[0,0,c],[0,0,0]]))
    True
    >>> matrix_is_nilpotent(matrix([[a,0,0],[0,b,0],[0,0,c]]))
    False
    """
    for v in M.eigenvalues():
        if not v.is_zero():
            return False
    return True

def jordan_cell_sizes(J):
    """Return a tuple of Jordan cell sizes from a matrix J in Jordan
    normal form.

    Examples:
    >>> jordan_cell_sizes(matrix([[1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,0,1]]))
    (2, 2)
    >>> jordan_cell_sizes(matrix([[1,1,0,0],[0,1,1,0],[0,0,1,0],[0,0,0,1]]))
    (3, 1)
    >>> jordan_cell_sizes(zero_matrix(5,5))
    (1, 1, 1, 1, 1)
    """
    assert J.is_square()
    sizes = []
    n = 1
    for i in range(J.nrows() - 1):
        if J[i, i+1].is_zero():
            sizes.append(n)
            n = 1
        else:
            n += 1
    sizes.append(n)
    assert sum(sizes) == J.nrows()
    return tuple(sizes)

def solve_right_fixed(A, B):
    """As of SageMath 6.10, 'Matrix.solve_right' method uses a broken
    check for solution correctness; this function corrects that check.
    """
    C = A.solve_right(B, check=False)
    if C.nrows() and C.ncols():
        # Sometimes 'solve_right' returns a giant expression,
        # which can be simplified into a tiny one; without this
        # simplification 'is_zero' call below may take forever
        # to finish.
        #
        # Note that the condition above is there to make sure
        # that an empty matrix is not passed into 'simplify_rational',
        # otherwise you'll get this error:
        #   TypeError: unable to make sense of Maxima expression
        #   'matrix()' in Sage
        C = C.simplify_rational()
    if not (A*C - B).simplify_rational().is_zero():
        raise ValueError("matrix equation has no solutions")
    return C

def solve_left_fixed(A, B):
    """As of SageMath 6.10, 'Matrix.solve_left' method uses a broken
    check for solution correctness; this function corrects that check.
    """
    return solve_right_fixed(A.transpose(), B.transpose()).transpose()

def any_integer(rng, ring, excluded):
    r = 2
    while True:
        p = ring(rng.randint(-r, r))
        if p not in excluded:
            return p
        r *= 2

#==================================================================================================
# Transformation routines
#==================================================================================================

def block_triangular_form(m):
    """Find a lower block-triangular form of a given matrix.

    Return a tuple `(M, T, B)` where:
      * `M` is a new matrix;
      * `T` is a transformation matrix;
      * `B` is a list of tuples (ki, ni), such that M's i-th
        diagonal block is given by `M.submatrix(ki, ki, ni, ni)`.
    """
    logger.enter("block_triangular_form")
    if logger.is_verbose():
        logger.debug("matrix before transformation:\n%s\n" % matrix_mask_str(m))

    n = m.nrows()
    deps_1_1 = {}
    for i, row in enumerate(m.rows()):
        deps_1_1[i] = set(j for j,ex in enumerate(row) if not bool(ex==0))

    deps_1_all = {}
    def find_deps_1_all(i):
        if i in deps_1_all:
            return deps_1_all[i]
        deps_1_all[i] = deps_1_1[i]
        for j in deps_1_1[i]:
            if i == j:
                continue
            find_deps_1_all(j)
            deps_1_all[i] = deps_1_all[i].union(deps_1_all[j])
    [find_deps_1_all(j) for j in range(n)]

    deps_coup = dict((i, set([])) for i in range(n))
    for i in range(n):
        for j in range(n):
            if (i in deps_1_all[j]) and (j in deps_1_all[i]):
                deps_coup[i].update([i,j])
                deps_coup[j].update([i,j])

    shuffle = []
    error = False
    while not error:
        if not deps_coup:
            break
        error = True
        for i in range(n):
            if i not in deps_coup:
                continue
            b = deps_coup[i].copy()
            if b != deps_1_all[i]:
                continue
            if b == set([]):
                b = set([i])
            shuffle += [b]
            for j in deps_1_all[i].copy():
                if j in deps_1_all:
                    del deps_1_all[j]
                    del deps_coup[j]
            for j in range(n):
                if j not in deps_coup:
                    continue
                deps_coup[j].difference_update(b)
                deps_1_all[j].difference_update(b)
            if i in deps_coup:
                del deps_coup[i]
            if i in deps_1_all:
                del deps_1_all[i]
            error = False
            break
    if error:
        raise FuchsiaError("Infinite loop")

    t = zero_matrix(SR,n)
    i = 0
    blocks = []
    for block in shuffle:
        blocks.append((i, len(block)))
        for j in list(block):
            t[j,i] = 1
            i += 1
    logger.info("found %d blocks" % len(blocks))
    mt = transform(m, None, t)
    if logger.is_verbose():
        logger.debug("matrix after transformation:\n%s\n" % matrix_mask_str(mt))
    logger.exit("block_triangular_form")
    return mt, t, blocks

def epsilon_form(m, x, eps, seed=0):
    logger.enter("epsilon_form")
    m, t1, b = block_triangular_form(m)
    m, t2 = reduce_diagonal_blocks(m, x, eps, b=b, seed=seed)
    m, t3 = fuchsify_off_diagonal_blocks(m, x, eps, b=b)
    m, t4 = factorize(m, x, eps, b=b, seed=seed)
    t = t1*t2*t3*t4
    logger.exit("epsilon_form")
    return m, t

def matrix_mask(m):
    n = m.nrows()
    return matrix(SR, n, n, [int(not bool(ex==0)) for ex in m.list()])

def matrix_mask_str(m):
    s = ''
    for row in matrix_mask(m).rows():
        s += ' '.join([str(ex) for ex in row]).replace('0', '.').replace('1','x') + "\n"
    return s

def matrix_str(m, n=2):
    buf = ""
    ind = " "*n
    for col in m.columns():
        for ex in col:
            buf += ind + str(ex) + "\n"
        buf += "\n"
    buf = buf[:-1]
    return buf

#==================================================================================================
# Step I: Fuchsify
#==================================================================================================

def is_fuchsian(m, x):
    for i, expr in enumerate(m.list()):
        if expr.is_zero():
            continue
        points = singularities_expr(expr, x)
        for x0, p in points.items():
            if p > 0:
                return False
    return True

def fuchsify(M, x, seed=0):
    """Given a system of differential equations of the form dF/dx=M*F,
    try to find a transformation T, which will reduce M to Fuchsian
    form. Return the transformed M and T. Raise FuchsiaError if
    M can not be transformed into Fuchsian form.

    Note that such transformations are not unique; you can obtain
    different ones by supplying different seeds.
    """
    logger.enter("fuchsify")
    assert M.is_square()
    rng = Random(seed)
    poincare_map = singularities(M, x)
    def iter_reductions(p1, U):
        for p2, prank2 in poincare_map.items():
            if bool(p2 == p1): continue
            while prank2 >= 0:
                B0 = matrix_c0(M, x, p2, prank2)
                if not B0.is_zero(): break
                poincare_map[p2] = prank2 = prank2 - 1
            if prank2 < 0: continue
            v = find_dual_basis_spanning_left_invariant_subspace(B0, U, rng)
            if v is None: continue
            P = fuchsia_simplify(U*v, x)
            M2 = fuchsia_simplify(balance_transform(M, P, p1, p2, x), x)
            yield p2, P, M2
    combinedT = identity_matrix(M.base_ring(), M.nrows())
    reduction_points = [pt for pt,p in poincare_map.items() if p >= 1]
    reduction_points.sort()
    if reduction_points == []:
        logger.info("already fuchsian")
    else:
        for pt in reduction_points:
            logger.info("rank = %d, x = %s" % (poincare_map[pt], pt))
    while reduction_points:
        pointidx = rng.randint(0, len(reduction_points) - 1)
        point = reduction_points[pointidx]
        prank = poincare_map[point]
        if prank < 1:
            del reduction_points[pointidx]
            continue
        while True:
            A0 = matrix_c0(M, x, point, prank)
            if A0.is_zero(): break
            A1 = matrix_c1(M, x, point, prank)
            try:
                U, V = alg1x(A0, A1, x)
            except FuchsiaError as e:
                logger.debug("Managed to fuchsify matrix to this state:\n"
                        "%s\nfurther reduction is pointless:\n%s" % (M, e))
                raise FuchsiaError("matrix cannot be reduced to Fuchsian form")
            try:
                point2, P, M = min(iter_reductions(point, U), \
                        key=lambda point2_P_M2: matrix_complexity(point2_P_M2[2]))
            except ValueError as e:
                point2 = any_integer(rng, M.base_ring(), poincare_map)
                P = fuchsia_simplify(U*V, x)
                M = balance_transform(M, P, point, point2, x)
                M = fuchsia_simplify(M, x)
                logger.info("Will introduce an apparent singularity at %s." % point2)
            logger.debug(
                "Applying balance between %s and %s with projector:\n%s" % (point, point2, P))
            combinedT = combinedT * balance(P, point, point2, x)
            if point2 not in poincare_map:
                poincare_map[point2] = 1
        poincare_map[point] = prank - 1
    combinedT = fuchsia_simplify(combinedT, x)
    logger.exit("fuchsify")
    return M, combinedT

def fuchsify_off_diagonal_blocks(m, x, eps, b=None):
    logger.enter("fuchsify_off_diagonal_blocks")
    n = m.nrows()
    if b is None:
        m, t, b = block_triangular_form(m)
    else:
        t = identity_matrix(SR, n)
    for i, (ki, ni) in enumerate(b):
        for j, (kj, nj) in enumerate(reversed(b[:i])):
            pts = singularities(m.submatrix(ki, kj, ni, nj), x)
            printed = False
            while any(pts.values()):
                bj = m.submatrix(ki, kj, ni, nj)
                if bj.is_zero():
                    break
                for x0, p in pts.items():
                    if p < 1:
                        continue
                    if not printed:
                        printed = True
                        logger.info("Fuchsifying block (%d, %d) (%d, %d)" % (ki, ni, kj, nj))
                        logger.debug("  singular points = %s" % (pts,))
                    a0 = matrix_residue(m.submatrix(ki, ki, ni, ni)/eps, x, x0)
                    b0 = matrix_c0(bj, x, x0, p)
                    c0 = matrix_residue(m.submatrix(kj, kj, nj, nj)/eps, x, x0)

                    d_vars = [gensym() for i in range(ni*nj)]
                    d = matrix(SR, ni, nj, d_vars)
                    eq = d + eps/p*(a0*d - d*c0) + b0/p
                    sol = solve(eq.list(), *d_vars, solution_dict=True)
                    d = d.subs(sol[0])

                    t0 = identity_matrix(SR, n)
                    t0[ki:ki+ni, kj:kj+nj] = \
                            d/(x-x0)**p if not (x0 == oo) else d*(x**p)
                    m = fuchsia_simplify(transform(m, x, t0), x)

                    t = fuchsia_simplify(t*t0, x)
                    pts[x0] -= 1
    logger.exit("fuchsify_off_diagonal_blocks")
    return m, t

def reduce_at_one_point(M, x, v, p, v2=oo):
    """Given a system of differential equations of the form dF/dx=M*F,
    with M having a singularity around x=v with Poincare rank p,
    try to find a transformation T, which will reduce p by 1 (but
    possibly introduce another singularity at x=v2). Return the
    transformed M and T.
    """
    assert M.is_square()
    assert p > 0
    n = M.nrows()
    combinedT = identity_matrix(M.base_ring(), n)
    while True:
        A0 = matrix_c0(M, x, v, p)
        if A0.is_zero(): break
        A1 = matrix_c1(M, x, v, p)
        U, V = alg1x(A0, A1, x)
        P = U*V
        M = balance_transform(M, P, v, v2, x)
        M = fuchsia_simplify(M, x)
        combinedT = combinedT * balance(P, v, v2, x)
    combinedT = fuchsia_simplify(combinedT, x)
    return M, combinedT

def find_dual_basis_spanning_left_invariant_subspace(A, U, rng):
    """Find matrix V, such that it's rows span a left invariant
    subspace of A and form a dual basis with columns of U (that
    is, V*U=I).
    """
    evlist = []
    for eigenval, eigenvects, evmult in A.eigenvectors_left():
        evlist.extend(eigenvects)
    rng.shuffle(evlist)
    W = matrix(evlist)
    try:
        M = solve_left_fixed(W*U, identity_matrix(U.ncols()))
        return M*W
    except ValueError:
        return None

def alg1x(A0, A1, x):
    #if not matrix_is_nilpotent(A0):
    #    raise FuchsiaError("matrix is irreducible (non-nilpotent residue)")
    # We will rely on undocumented behavior of jordan_form() call:
    # we will take it upon faith that it sorts Jordan cells by
    # their size, and puts the largest ones closer to (0, 0).
    A0J, U = A0.jordan_form(transformation=True)
    invU = U.inverse()
    A0J_cs = jordan_cell_sizes(A0J)
    assert all(A0J_cs[i] >= A0J_cs[i+1] for i in range(len(A0J_cs) - 1))
    ncells = len(A0J_cs)
    A0J_css = [sum(A0J_cs[:i]) for i in range(ncells + 1)]
    nsimplecells = sum(1 if s == 1 else 0 for s in A0J_cs)
    u0 = [U[:,A0J_css[i]] for i in range(ncells)]
    v0t = [invU[A0J_css[i+1]-1,:] for i in range(ncells)]
    L0 = matrix([
        [(v0t[k]*(A1)*u0[l])[0,0] for l in range(ncells)]
        for k in range(ncells)
    ])
    L0 = fuchsia_simplify(L0, x)
    L1 = matrix([
        [(v0t[k]*u0[l])[0,0] for l in range(ncells)]
        for k in range(ncells)
    ])
    assert (L1 - diagonal_matrix(
            [0]*(ncells-nsimplecells) + [1]*nsimplecells)).is_zero()
    #zero_rows = [i for i in range(A0.nrows()) if A0J[i,:].is_zero()]
    #zero_cols = [j for j in range(A0.nrows()) if A0J[:,j].is_zero()]
    #assert len(zero_rows) == len(zero_cols) == ncells
    #_L0 = (invU*A1*U)[zero_rows, zero_cols]
    #assert (L0 - _L0).is_zero()
    lam = SR.symbol()
    if not (L0 - lam*L1).determinant().is_zero():
        raise FuchsiaError("matrix is Moser-irreducible")
    S, D = alg1(L0, A0J_cs)
    I_E = identity_matrix(D.base_ring(), A0.nrows())
    for i in range(ncells):
        for j in range(ncells):
            if not D[i,j].is_zero():
                ni = A0J_css[i]
                nj = A0J_css[j]
                for k in range(min(A0J_cs[i], A0J_cs[j])):
                    I_E[ni+k,nj+k] += D[i,j]
    U_t = U*I_E
    invU_t = U_t.inverse()
    return \
        U_t[:, [A0J_css[i] for i in S]], \
        invU_t[[A0J_css[i] for i in S], :]

def alg1(L0, jordan_cellsizes):
    assert L0.is_square()
    ring = L0.base_ring()
    N = L0.nrows()
    I_N = identity_matrix(N)
    S = set()
    D = zero_matrix(ring, N)
    K = sum(1 if s > 1 else 0 for s in jordan_cellsizes)
    while True:
        Lx = copy(L0)
        for i in S:
            Lx[:, i] = 0
            Lx[i, :] = 0
        fi = None
        c = None
        for i in range(0, N):
            if i not in S:
                try:
                    c = solve_right_fixed(Lx[:,0:i], Lx[:,i])
                except ValueError:
                    # No solution found; vectors are independent.
                    continue
                fi = i
                break
        assert fi is not None
        assert c is not None
        D0 = matrix(ring, N)
        invD0 = matrix(ring, N)
        for j in range(fi):
            D0[j, fi] = -c[j, 0]
            invD0[j, fi] = -c[j, 0] \
                    if jordan_cellsizes[j] == jordan_cellsizes[fi] else 0
        L0 = (I_N - invD0)*L0*(I_N + D0)
        D = D + D0 + D*D0
        if fi < K:
            break
        S.add(fi)
    # Check Alg.1 promise
    for j in range(N):
        for k in range(N):
            if (j not in S) and ((k in S) or k == fi):
                assert L0[j,k] == 0
    S.add(fi)
    return S, D

#==================================================================================================
# Step II: Normalize
#==================================================================================================

def is_normalized(M, x, eps):
    """Return True if (a Fuchsian) matrix M is normalized, that
    is all the eigenvalues of it's residues in x lie in [-1/2, 1/2)
    range (in limit eps->0). Return False otherwise.

    Examples:
    >>> x, e = var("x epsilon")
    >>> is_normalized(matrix([[(1+e)/3/x, 0], [0, e/x]]), x, e)
    True
    """
    points = singularities(M, x)
    for x0, p in points.items():
        M0 = matrix_residue(M, x, x0)
        for ev in M0.eigenvalues():
            ev = limit_fixed(ev, eps, 0)
            if not (Rational((-1, 2)) <= ev and ev < Rational((1, 2))):
                return False
    return True

def are_diagonal_blocks_reduced(m, b, x, eps):
    """Return True if diagonal blocks of `m` are normalized, that is all the eigenvalues of theirs
    residues in `x` lie in the range [-1/2, 1/2) in eps->0 limit; return False otherwise. Diagonal
    blocks are defined by the list `b` which corresponds to the equivalent value returned by the
    `block_triangular_form` routine.

    Examples:
    >>> x, e = var("x epsilon")
    >>> are_diagonal_blocks_reduced(matrix([[(1+e)/3/x, 0], [0, e/x]]), [(0,1),(1,1)], x, e)
    True
    """
    for ki, ni in b:
        mi = fuchsia_simplify(m.submatrix(ki, ki, ni, ni), x)
        if not is_normalized(mi, x, eps):
            return False
    return True

def reduce_diagonal_blocks(m, x, eps, b=None, seed=0):
    """Given a lower block-triangular system of differential equations of the form dF/dx=m*F,
    find a transformation that will shift all eigenvalues of all residues of all its diagonal
    blocks into the range [-1/2, 1/2) in eps->0 limit. Return the transformed matrix m and the
    transformation; raise FuchsiaError if such transformation is not found. Diagonal blocks
     are defined by the list `b` which corresponds to the equivalent value returned by the
    `block_triangular_form` routine.
    """
    logger.enter("reduce_diagonal_blocks")
    n = m.nrows()
    if b is None:
        m, t, b = block_triangular_form(m)
    else:
        t = identity_matrix(SR, n)
    for i, (ki, ni) in enumerate(b):
        mi = fuchsia_simplify(m.submatrix(ki, ki, ni, ni), x)
        ti = identity_matrix(SR, ni)
        logger.info("reducing block #%d (%d,%d)" % (i,ki,ni))
        if logger.is_verbose():
            logger.debug("\n%s" % matrix_str(mi, 2))

        mi_fuchs, ti_fuchs = fuchsify(mi, x, seed=seed)
        ti = ti*ti_fuchs

        mi_norm, ti_norm = normalize(mi_fuchs, x, eps, seed=seed)
        ti = ti*ti_norm

        mi_eps, ti_eps = factorize(mi_norm, x, eps, seed=seed)
        ti = ti*ti_eps

        t[ki:ki+ni, ki:ki+ni] = ti
    mt = fuchsia_simplify(transform(m, x, t), x)
    logger.exit("reduce_diagonal_blocks")
    return mt, t


def cmap_add_div(Cmap, C, pi, ki, x0):
    # M += C*(x-pi)^ki/(x-x0)
    if pi == x0:
        Cmap[pi, ki-1] += C
    elif ki >= 0:
        assert pi == 0
        for k in range(0, ki):
            Cmap[pi, k] += C * x0**(ki-1-k)
        Cmap[x0, -1] += C * x0**ki
    else:
        for k in range(ki, 0):
            Cmap[pi, k] -= C * (x0-pi)**(ki-1-k)
        Cmap[x0, -1] += C * (x0-pi)**ki

def cmap_add_mul(Cmap, C, pi, ki, x0):
    # M += C*(x-pi)^ki*(x-x0) = C*(x-pi)^ki*{(x-pi) + (pi-x0)}
    Cmap[pi, ki] += C*(pi-x0)
    if ki == -1:
        Cmap[SR(0), ki+1] += C
    else:
        Cmap[pi, ki+1] += C

def memoize(f):
    @wraps(f)
    def wrapper(self):
        key = f.__name__+"_cache"
        if not hasattr(self, key):
#            print("cache miss: %s" % key)
            setattr(self, key, f(self))
#        else:
#            print("cache hit: %s" % key)
        return getattr(self, key)
    return wrapper

class FuchsianSystem(object):
    @staticmethod
    def from_M(m, x, eps):
        # singular points
        xs = sorted(singularities(m, x).keys())
        # residues in singular points
        rs = {}
        for xi in xs:
            rs[xi] = SystemResidue(matrix_residue(m, x, xi), eps)
        # TODO: check if system's residues have half-integer eigenvalues.
        # If yes, then inform user about that and let us remove resonances only.
        # Though, such a form does not allow to construct an eps-form, it
        # however can be used to find solutions as a power series using
        # semi-analytical methods.
        return FuchsianSystem(rs, x, eps)

    def __init__(self, rs, x, eps):
        self.rs = rs
        self.xs = sorted(rs.keys())
        self.x = x
        self.eps = eps

    def __str__(self):
        #s = "zeroes:\n  %s\n" % (self.xs,)
        s = ""
        a = 0
        for xi, ri in self.rs.items():
            s += "a = %s (x = %s)\n%s\n" % (a, xi, ri)
            a += 1
        return s

    def write(self, filename):
        m = SR(0)
        for xi, ri in self.rs.items():
            if xi == oo:
                continue
            m += ri.m/(self.x-xi)
        export_matrix_to_file(filename, m)

    def _balance(self, i1, j1, i2, j2):
        x1, x2 = self.xs[i1], self.xs[i2]
        r1, r2 = self.rs[x1], self.rs[x2]
        v1, v2 = r1.eigenvectors_right()[j1][1], r2.eigenvectors_left()[j2][1]
        print("x1 = %s" % x1)
        print("x2 = %s" % x2)
        print("v1 = %s" % v1)
        print("v2 = %s" % v2)

        sp = simple(dot_product(v1, v2))
        print("sp = %s" % sp)
        P = simple(cross_product(v1, v2) / sp)
        print("P  =\n%s" % P)
        m = balance_transform(self.m, P, x1, x2, self.x)
        return FuchsianSystem.from_M(m, self.x, self.eps)

    #@logcall
    def balance(self, i1, j1, i2, j2):
        x1, x2 = self.xs[i1], self.xs[i2]
        r1, r2 = self.rs[x1], self.rs[x2]
        v1, v2 = r1.eigenvectors_right()[j1][1], r2.eigenvectors_left()[j2][1]

        sp = simple(dot_product(v1, v2))
        print("sp = %s" % sp)
        P = simple(cross_product(v1, v2) / sp)
        print("P  =\n%s" % P)

        coP = 1-P
        if x1 == oo:
            Cmap = defaultdict(SR.zero)
            for pi, Ci in self.rs.items():
                if pi == oo:
                    continue
                # coP Ci coP (x-pi)^-1 + P Ci P (x-pi)^-1
                Cmap[pi, -1] += coP*Ci.m*coP + P*Ci.m*P
                # coP Ci P -(x-x2) (x-pi)^-1
                cmap_add_mul(Cmap, -coP*Ci.m*P, pi, -1, x2)
                # P Ci coP -1/(x-x2) (x-pi)^-1
                cmap_add_div(Cmap, -P*Ci.m*coP, pi, -1, x2)
            # -P/(x-x2)
            Cmap[x2, -1] -= P
        elif x2 == oo:
            Cmap = defaultdict(SR.zero)
            for pi, Ci in self.rs.items():
                if pi == oo:
                    continue
                # coP Ci coP (x-pi)^-1 + P Ci P (x-pi)^-1
                Cmap[pi, -1] += coP*Ci.m*coP + P*Ci.m*P
                # P Ci coP -(x-x1) (x-pi)^-1
                cmap_add_mul(Cmap, -P*Ci.m*coP, pi, -1, x1)
                # coP Ci P -1/(x-x1) (x-pi)^-1
                cmap_add_div(Cmap, -coP*Ci.m*P, pi, -1, x1)
            # P/(x-x1)
            Cmap[x1, -1] += P
        else:
            Cmap = defaultdict(SR.zero, [((pi,-1),ri.m) for pi,ri in self.rs.items()])
            for pi, Ci in self.rs.items():
                if pi == oo:
                    continue
                # coP Ci P (x1-x2)/(x-x1) (x-pi)^-1
                cmap_add_div(Cmap, coP*Ci.m*P*(x1-x2), pi, -1, x1)
                # P Ci coP (x2-x1)/(x-x2) (x-pi)^-1
                cmap_add_div(Cmap, P*Ci.m*coP*(x2-x1), pi, -1, x2)
            # P/(x-x1) - P/(x-x2)
            Cmap[x1, -1] += P
            Cmap[x2, -1] -= P
        rs = {}
        Coo = SR.zero()
        for key, C in Cmap.items():
            pi, ki = key
            if pi == oo:
                continue
            Cmap[key] = C = simple(C)
            if C.is_zero():
                del Cmap[key]
            else:
                assert(ki == -1)
                rs[pi] = SystemResidue(C)
                Coo -= C
        Coo = simple(Coo)
        rs[oo] = SystemResidue(Coo)
        #T = fuchsia_simplify(self.T * balance(P, x1, x2, self.x))
        #return RationalSystem(dict(Cmap), self.x, T, self.trace + [("balance", P, x1, x2)])
        return FuchsianSystem(rs, self.x, self.eps) 

    #@logcall
    @memoize
    def balances(self):
        res = []
        for i1,x1 in enumerate(self.xs):
            for i2,x2 in enumerate(self.xs):
                if i1 == i2:
                    continue
                for j1, j2, aux in self.balances2(x1, x2):
                    res.append((i1,j1, i2,j2, aux))
        return res

    #@logcall
    def balances2(self, x1, x2):
        #assert(not (x1-x2).is_zero())
        res = []
        r1, r2 = self.rs[x1], self.rs[x2]
        for j1, (a1,v1) in enumerate(r1.eigenvectors_right()):
            for j2, (a2,v2) in enumerate(r2.eigenvectors_left()):
                sp = simple(dot_product(v1, v2))
                if sp.is_zero():
                    continue
                res.append((j1, j2, (a1,a2, sp)))
        return res

    #@logcall
    @memoize
    def eigenvalues(self):
        res = []
        for xi in self.xs:
            res.append((xi, self.rs[xi].eigenvalues()))
        return res

#    @memoize
#    @logcall
#    def eigenvectors_left(self):
#        return self._eigenvectors(left=True)
#
#    @memoize
#    @logcall
#    def eigenvectors_right(self):
#        return self._eigenvectors(left=False)
#
#    def _eigenvectors(self, left=True):
#        res = []
#        for xi, ri in self.rs.items():
#            eigenvectors = ri.eigenvectors_left
#            if left is not True:
#                eigenvectors = ri.eigenvectors_right
#            res.append((xi, eigenvectors()))
#        return res

    def normalize(self):
        pass

class SystemResidue(object):
    def __init__(self, m, eps=None):
        self.m = m
        self.eps = eps

    def __str__(self):
        return str(self.m)

    def __repr__(self):
        return self.m.__repr__()

    @logcall
    @memoize
    def eigenvalues(self):
        return simple(mathematica.Eigenvalues(self.m).sage())

    @logcall
    @memoize
    def eigenvectors_left(self):
        return self._eigenvectors(left=True)

    @logcall
    @memoize
    def eigenvectors_right(self):
        return self._eigenvectors(left=False)

    def _eigenvectors(self, left=True):
        m = self.m.transpose() if left is True else self.m
        s = simple(mathematica.Eigensystem(m).sage())
        res = zip(s[0], s[1])
        return res

#    def _eigenvectors(self, left=True):
#        res = []
#        eigenvectors = self.m.eigenvectors_left
#        if left is not True:
#            eigenvectors = self.m.eigenvectors_right
#        res = [vvm[0:2] for vvm in eigenvectors()]
#        return res

# helpers
@logcall
def simple(obj, x=None):
    def _simplify(ex):
            return mathematica.Simplify(ex).sage()
    if hasattr(obj, "apply_map"):
        return obj.apply_map(_simplify)
    else:
        return _simplify(obj)

def ex_in(ex, lst):
    #if any((ex-t).is_zero() for t in printed):
    for t in lst:
        if (ex-t).is_zero():
            return True
    return False
        

class NormalizeAssistant(object):
    def __init__(self, m):
        self.m = m           # current matrix to balance
        self.history = []    # matrices balanced so far

    # Entry point
    def start(self):
#        cmd = self.cmd_main
#        while cmd != None:
#            cmd = cmd()
#        return
        try:
            cmd = self.cmd_main
            while cmd != None:
                cmd = cmd()
        except Exception as error:
            print("UNEXPECTED ERROR:\n%s" % (error,))
            return self.start()

    # Helpers
    def choose_balance(self, i1, i2):
        x1, x2 = self.m.xs[i1], self.m.xs[i2]
        bs = self.m.balances2(x1, x2)
        for i, (j1, j2, aux) in enumerate(bs):
            print("[%s] %s" % (i, aux))
        n = raw_input("(choose balance) > ")
        try:
            n = int(n)
        except ValueError:
            return None
        if n < 0 or n > len(bs):
            return None
        b = bs[n]
        return (i1, b[0], i2, b[1])

    # Commands
    def cmd_balance(self, b):
        print("balancing with %s" % (b,))

        mm = apply(self.m.balance, b)
        self.history.append(mm)
        self.m = mm
        return self.cmd_main

    def cmd_balance_res(self, args):
        if len(args) != 2:
            return self.cmd_error("needs exactly 2 arguments")
        try:
            i1, i2 = [int(n) for n in args]
        except ValueError:
            return self.cmd_error("arguments out of range (type e for valid values)")
        b = self.choose_balance(i1, i2)
        if b is None:
            return self.cmd_error("unknown balance")
        return self.cmd_balance(b[0:4])

    def cmd_error(self, msg):
        print("ERROR: %s" % msg)
        return self.cmd_main

    def cmd_main(self):
        cmd = raw_input("> ").split(' ')
        cmd, args = cmd[0], cmd[1:]
        if cmd == 's':
            return self.cmd_save(args)
        elif cmd == 'b':
            return self.cmd_balance_res(args)
        elif cmd == 'e':
            return self.cmd_print_eigenvalues
        elif cmd == 'h':
            return self.cmd_print_help
        elif cmd == 'q':
            return None
        elif cmd == '':
            return self.cmd_print_help
        return self.cmd_unknown_command(cmd)

    def cmd_print_eigenvalues(self):
        for i, (xi, a) in enumerate(self.m.eigenvalues()):
            print("[%s] x = %s\n    %s" % (i, xi, a))
        return self.cmd_main

    def cmd_print_help(self):
        """
        Available commands:
            b <i> <j>    balance xi, xj points
            e            show eigenvectors
            h            print this help
            q            quit
            s <file>     save matrix to file
        """
        print(inspect.cleandoc(self.cmd_print_help.__doc__))
        return self.cmd_main

    def cmd_save(self, args):
        if len(args) != 1:
            return self.cmd_error("needs exactly 1 arguments")
        name = args[0]
        self.m.write(name)
        return self.cmd_main

    def cmd_unknown_balance(self, n):
        print("Unknown balance: '%s'" % n)
        return self.cmd_main

    def cmd_unknown_command(self, cmd):
        print("Unknown command: '%s'" % cmd)
        return self.cmd_print_help

def normalize(m, x, eps, seed=0):
    """Given a Fuchsian system of differential equations of the
    form dF/dx=m*F, find a transformation that will shift all
    the eigenvalues of m's residues into [-1/2, 1/2) range (in
    the limit eps->0). Return the transformed matrix m and the
    transformation. Raise FuchsiaError if such transformation
    is not found.
    """
    class State(object):
        def __init__(self, m, x, eps, seed):
            # random
            self.random = Random(seed)
            # points
            self.points = singularities(m, x).keys()
            # ev_cum
            self.ev_cum = {}
            for x0 in self.points:
                pos, neg = 0, 0
                a0 = matrix_residue(m, x, x0)
                for ev in a0.eigenvalues():
                    ev0 = limit_fixed(ev, eps, 0)
                    if ev0 > 0:
                        pos += ev0
                    if ev0 < 0:
                        neg += ev0
                    logger.debug("x = %s, ev = %s" % (x0, ev0))
                self.ev_cum[x0] = [pos,neg]
            # x0
            logger.debug("state.ev_cum = %s" % self.ev_cum)
            self.x0 = None

        def is_normalized(self):
            for ev in self.ev_cum.values():
                if ev != [0,0]:
                    return False
            return True

        def pairs(self):
            points = [x0 for x0 in self.points if self.ev_cum[x0] != [0,0]]
            if len(points) == 1:
                self.select_x0(points[0])
                points.append(self.x0)
            return permutations(points, 2)

        def select_x0(self, x1):
            if self.x0 is not None:
                return
            for x0 in self.points:
                if x0 != x1:
                    self.x0 = x0
                    return

        def update_ev_cum(self, x2, x1):
            ev_cum_x1, ev_cum_x2 = self.ev_cum[x1], self.ev_cum[x2]
            if ev_cum_x1[0] > 0:
                ev_cum_x1[0] -= 1
            else:
                ev_cum_x1[1] -= 1
            if ev_cum_x2[1] < 0:
                ev_cum_x2[1] += 1
            else:
                ev_cum_x2[0] += 1

    logger.enter("normalize")
    state = State(m, x, eps, seed)
    T = identity_matrix(m.base_ring(), m.nrows())
    if state.is_normalized():
        logger.info("already normalized")
    i = 0
    while not state.is_normalized():
        i += 1
        m = fuchsia_simplify(m, x)
        logger.info("step %s" % i)
        balances = find_balances(m, x, eps, state)
        b = select_balance(balances, eps, state)
        if b is None:
            raise FuchsiaError("can not balance matrix")
        logger.info("  balancing x = %s and x = %s" % (b[1],b[2]))
        if logger.is_verbose():
            logger.debug("\n      use the balance:\n        %s\n" % b)

        cond, x1, x2, a0_eval, b0_eval, a0_evec, b0_evec, scale = b
        if cond == 1:
            P = cross_product(a0_evec, b0_evec) / scale
            m = balance_transform(m, P, x1, x2, x)
            T0 = balance(P, x1, x2, x)
            state.update_ev_cum(x1, x2)
        else:
            P = cross_product(b0_evec, a0_evec) / scale
            m = balance_transform(m, P, x2, x1, x)
            T0 = balance(P, x2, x1, x)
            state.update_ev_cum(x2, x1)

        T = fuchsia_simplify(T*T0, x)
    logger.exit("normalize")
    return m, T

def find_balances(m, x, eps, state={}):
    residues = {}
    for x1, x2 in state.pairs():
        logger.debug("trying to balance x = %s and x = %s" % (x1,x2))
        for xi in [x1,x2]:
            if xi not in residues:
                residues[xi] = matrix_residue(m, x, xi)
        a0, b0 = residues[x1], residues[x2]

        a0_evr, b0_evl = eigenvectors_right(a0), eigenvectors_left(b0)

        if logger.is_verbose():
            msg = "\n  Eigenvalues:\n"
            msg += "    x = %s:\n" % x1
            a0_evals = [];
            for ev, evec, emult in a0_evr:
                a0_evals += [ev]*emult
            msg += "        %s\n" % str(a0_evals)
            msg += "    x = %s:\n" % x2
            b0_evals = [];
            for ev, evec, emult in b0_evl:
                b0_evals += [ev]*emult
            msg += "        %s\n" % str(b0_evals)
            logger.debug(msg)

        balances_1 = find_balances_by_cond(a0_evr, b0_evl, lambda a0_eval, b0_eval: limit_fixed(a0_eval, eps, 0) < -0.5)
        for balance in balances_1:
            balance = [1, x1, x2] + balance
            yield balance

        a0_evl, b0_evr = eigenvectors_left(a0), eigenvectors_right(b0)
        balances_2 = find_balances_by_cond(a0_evl, b0_evr, lambda a0_eval, b0_eval: limit_fixed(a0_eval, eps, 0) >= 0.5)
        for balance in balances_2:
            balance = [2, x1, x2] + balance
            yield balance

def find_balances_by_cond(a0_ev, b0_ev, cond):
    res = []
    for a0_eval, a0_evecs, a0_evmult in a0_ev:
        for b0_eval, b0_evecs, b0_evmult in b0_ev:
            if not cond(a0_eval, b0_eval):
                logger.debug("Balance rejected:\n    a0_eval = %s\n    b0_eval = %s" % (a0_eval, b0_eval))
                continue
            for a0_evec in a0_evecs:
                for b0_evec in b0_evecs:
                    scale = fuchsia_simplify(dot_product(a0_evec, b0_evec))
                    balance = [a0_eval, b0_eval, a0_evec, b0_evec, scale]
                    if scale == 0:
                        logger.debug("Balance rejected:\n    a0_eval = %s\n    b0_eval = %s\n    a0_evec = %s\n    b0_evec = %s\n    scale   = %s" % tuple(balance))
                        continue
                    logger.debug("Balance found:\n    a0_eval = %s\n    b0_eval = %s\n    a0_evec = %s\n    b0_evec = %s\n    scale   = %s" % tuple(balance))
                    res.append(balance)
    return res

def select_balance(balances, eps, state={}):
    min_degree, min_balance = None, None
    bs = []
    for b in balances:
        cond, x1, x2, a0_eval, b0_eval, a0_evec, b0_evec, scale = b
        if (cond == 1) and limit_fixed(a0_eval, eps, 0) < -0.5 and \
                limit_fixed(b0_eval, eps, 0) >= 0.5:
            degree = max(scale.numerator().degree(eps), scale.denominator().degree(eps))
            if degree < 4:
                return b
            if (min_degree is None) or (min_degree > degree):
                 min_degree = degree
                 min_balance = b
        elif (cond == 2) and limit_fixed(a0_eval, eps, 0) >= 0.5 and \
                limit_fixed(b0_eval, eps, 0) < -0.5:
            degree = max(scale.numerator().degree(eps), scale.denominator().degree(eps))
            if degree < 4:
                return b
            if (min_degree is None) or (min_degree > degree):
                 min_degree = degree
                 min_balance = b
        bs.append(b)
    if min_balance is not None:
        return min_balance

    x0 = state.x0
    if x0 is None:
        for b in bs:
            cond, x1, x2, ev1, ev2 = b[:5]
            if cond == 1:
                x0 = x2
                break
            if cond == 2:
                x0 = x1
                break
        logger.info("  select x0 = %s" % x0)
        state.x0 = x0

    balances_x0 = [b for b in bs if (b[0] == 1 and b[2] == x0) or (b[0] == 2 and b[1] == x0)]
    b = state.random.choice(balances_x0) if balances_x0 else None
    return b

def eigenvectors_left(m):
    if m._cache is None:
        m._cache = {}
    if "eigenvectors_left" not in m._cache:
        res = simplify(m.eigenvectors_left())
        m._cache["eigenvectors_left"] = res
    return m._cache["eigenvectors_left"]

def eigenvectors_right(m):
    if m._cache is None:
        m._cache = {}
    if "eigenvectors_right" not in m._cache:
        res = m.eigenvectors_right()
        m._cache["eigenvectors_right"] = res
    return m._cache["eigenvectors_right"]

#==================================================================================================
# Step III: Factorize
#==================================================================================================

def gensym():
    sym = SR.symbol()
    SR.symbols[str(sym)] = sym
    return sym

def factorize(M, x, epsilon, b=None, seed=0):
    """Given a normalized Fuchsian system of differential equations:
        dF/dx = M(x,epsilon)*F,
    try to find a transformation that will factor out an epsilon
    from M. Return a transformed M (proportional to epsilon)
    and T. Raise FuchsiaError if epsilon can not be factored.
    """
    logger.info("-> factorize")
    n = M.nrows()
    M = fuchsia_simplify(M, x)
    if epsilon not in expand(M/epsilon).variables():
        logger.info("   already in epsilon form")
        logger.info("<- factorize")
        return M, identity_matrix(SR, n)
    rng = Random(seed)
    mu = gensym()
    if b is None:
        T_symbols = [gensym() for i in range(n*n)]
        T = matrix(SR, n, n, T_symbols)
    else:
        T, T_symbols = matrix(SR, n), []
        for ki,ni in b:
            for i in range(ki,ki+ni):
                for j in range(ki+ni):
                    sym = gensym()
                    T[i,j] = sym
                    T_symbols.append(sym)
    eqs = []
    for point, prank in singularities(M, x).items():
        assert prank == 0
        logger.debug("    processing point x = %s" % point)
        R = matrix_c0(M, x, point, 0)
        R = fuchsia_simplify(R)
        eq = (R/epsilon)*T-T*(R.subs({epsilon: mu})/mu)
        eq = fuchsia_simplify(eq)
        eqs.extend(eq.list())
    logger.info("    found %d equations with %d unknowns" % (len(eqs), len(T_symbols)))
    solutions = fuchsia_solve(eqs, T_symbols)
    for solution in solutions:
        S = T.subs(solution)
        # Right now S likely has a number of free variables in
        # it; we can set them to arbibtrary values, as long as
        # it'll make S invertible.
        rndrange = 0
        while True:
            try:
                sT = S.subs([
                    e==rng.randint(-rndrange, rndrange)
                    for e in S.variables() if e != epsilon
                ])
                sT = fuchsia_simplify(sT,x)
                M = fuchsia_simplify(transform(M, x, sT), x)
                # We're leaking a bunch of temprary variables here,
                # which accumulate in SR.variables, but who cares?
                logger.info("<- factorize")
                return M, sT
            except (ZeroDivisionError, ValueError):
                rndrange += 1 + rndrange//16
                # We've tried a bunch of substituions, and they didn't
                # work. Is the matrix at all invertible? Let's check.
                if rndrange == 16 and not S.is_invertible():
                    break
    raise FuchsiaError("can not factor epsilon")

#==================================================================================================
# Helpers
#==================================================================================================

def matrix_complexity(M):
    return len(str(M.list()))

def simplify_by_jordanification(M, x):
    """Try to simplify matrix M by constant transformations that
    transform M's residues into their Jordan forms. Return the
    simplified matrix and the transformation. If none of the
    attempted transformations reduce M's complexity (as measured
    by 'matrix_complexity()'), return the original matrix and
    the identity transformation.
    """
    minM = M
    minC = matrix_complexity(M)
    minT = identity_matrix(M.base_ring(), M.nrows())
    for point, prank in singularities(M, x).items():
        R = matrix_c0(M, x, point, prank)
        J, T = R.jordan_form(transformation=True)
        MM = fuchsia_simplify(transform(M, x, T), x)
        C = matrix_complexity(MM)
        if C < minC:
            minM = MM
            minC = C
            minT = T
    return minM, minT

def common_factor(expressions, filter):
    """Factorize given expressions, select those factors for
    which 'filter(factor)' is True, and return the product of
    factors common to all the expressions.

    Examples:
    >>> x = var("x")
    >>> common_factor([x*x-1, x+1], lambda f: True)
    x + 1
    >>> common_factor([1/x**2, 2/x**3, 3/x**4], lambda f: True)
    x^(-2)

    Note that if there is a mix of positive and negative exponents
    of a given factor, this function will use (one of) the most
    frequently occurring exponent:
    >>> common_factor([x, 1/x, 2/x**2, 3/x], lambda f: True)
    1/x
    """
    factor2exp2count = defaultdict(lambda: defaultdict(lambda: 0))
    for i, expr in enumerate(expressions):
        factors = dict(expr.factor_list())
        for factor, n in factors.items():
            if not filter(factor): continue
            if factor in factor2exp2count:
                factor2exp2count[factor][n] += 1
            else:
                if i > 0: factor2exp2count[factor][0] = i
                factor2exp2count[factor][n] = 1
        for factor, exps in factor2exp2count.items():
            if factor not in factors:
                exps[0] += 1
    result = SR(1)
    for factor, exp2count in factor2exp2count.items():
        exps = exp2count.keys()
        minn = min(exps)
        maxn = max(exps)
        if minn > 0: result *= factor**minn
        if maxn < 0: result *= factor**maxn
        if minn <= 0 and maxn >= 0:
            bestn = max(exps, key=lambda exp: exp2count[exp])
            result *= factor**bestn
    return result

def simplify_by_factorization(M, x):
    """Try to simplify matrix M by a constant transformation
    that extracts common factors found in M (if any). Return
    the simplified matrix and the transformation.
    """
    n = M.nrows()
    T = identity_matrix(M.base_ring(), n)
    factors = []
    for i in range(n):
        factor = common_factor(
            [M[i,k] for k in range(n) if i != k and not M[i,k].is_zero()] +
            [1/M[k,i] for k in range(n) if k != i and not M[k,i].is_zero()],
            lambda e: x not in e.variables())
        if factor != 1:
            dT = identity_matrix(M.base_ring(), n)
            dT[i,i] = T[i,i] = factor
            M = transform(M, x, dT)
    logger.debug(
            "stripping common factors with this transform:\ndiagonal_matrix([\n  %s\n])" % ",\n  ".join(str(e) for e in T.diagonal()))
    return fuchsia_simplify(M, x), T

#==============================================================================
# Import/Export routines
#==============================================================================

_parser = Parser(make_int=SR, make_float=SR,
        make_var=lambda v: I if v in "Ii" else SR.var(v),
        make_function={"sqrt": sqrt})

def parse(s):
    """Convert a given string to a Sage expression.

    >>> parse("(1 + I*eps)*(1 - I*eps)").simplify_rational()
    eps^2 + 1
    """
    return _parser.parse(s)

def import_matrix_from_file(filename):
    """Read and return a matrix from a named file. Both Mathematica
    and MatrixMarket array formats are supported. The exact format
    will be autodetected.
    """
    with open(filename, 'r') as f:
        hd = f.read(2)
        f.seek(0)
        if hd == "%%":
            return import_matrix_matrixmarket(f)
        elif hd == "{{":
            return import_matrix_mathematica(f)
        else:
            raise ValueError("File '%s' is not in a known matrix format" % filename)

def import_matrix_mathematica(f):
    """Read and return a matrix, stored in the Mathematica format,
    from a file-like object.
    """
    data =  f.read()
    data = data.replace(' ', '').replace('\n', '')
    nr = data.count('},{')+1
    data = data.replace('{', '').replace('}', '')
    data = data.replace(',', '\n')
    data = data.replace('Sqrt[', 'sqrt(').replace(']', ')')
    nc = int((data.count('\n')+1)/nr)

    import io
    sio = io.StringIO("%d %d\n" % (nc, nr) + data)
    return import_matrix_matrixmarket(sio).T

def import_matrix_matrixmarket(f):
    """Read and return a matrix, stored in the MatrixMarket
    array format, from a file-like object.
    """
    lineno = 0
    while True:
        s = f.readline()
        lineno += 1
        if not s.startswith('%'):
            break
    try:
        nrows, ncols = map(int, s.split())
    except ValueError:
        msg = ("Can not determine size of the matrix\n"
               "Make sure the file is in the MatrixMarket array format")
        raise FuchsiaError(msg)
    data = []
    for s in f.readlines():
        lineno += 1
        if s.startswith('%'):
            continue
        try:
            expr = parse(s)
        except SyntaxError:
            msg = ("Malformad expression at line %d:\n%s\n"
                   "Make sure the file is in the MatrixMarket array format")
            raise FuchsiaError(msg % (lineno, s))
        data.append(expr)
    m = matrix(ncols, nrows, data).T
    return m

def export_matrix_to_file(filename, m, fmt="m"):
    """Save matrix into a file in a given matrix format. Set
    'fmt' to "mtx" for MatrixMarket array format, or to "m"
    for Mathematica format.
    """
    with open(filename, 'w') as f:
        if fmt == "mtx":
            export_matrix_matrixmarket(f, m)
        elif fmt == "m":
            export_matrix_mathematica(f, m)
        else:
            raise FuchsiaError("Unknown matrix format '%s'" % fmt)

def export_matrix_mathematica(f, m):
    """Write matrix to a file-like object using Mathematica format."""
    i = 0
    f.write("{")
    nc, nr = m.ncols(), m.nrows()
    for row in m.rows():
        i += 1
        f.write("{")
        j = 0
        for mij in row:
            j += 1
            f.write(str(mij).replace(' ', ''))
            if j < nc:
                f.write(',')
        f.write("}")
        if i < nr:
            f.write(",")
    f.write('}')

def export_matrix_matrixmarket(f, m):
    """Write matrix to a file-like object using MatrixMarket
    array format.
    """
    f.write("%%MatrixMarket matrix array Fuchsia[symbolic] general\n")
    f.write("%d %d\n" % (m.nrows(), m.ncols()))
    i = 0
    for col in m.columns():
        i += 1
        f.write('%% #%d\n' % i)
        for mij in col:
            f.write(str(mij).replace(' ', ''))
            f.write('\n')

#==================================================================================================

def main():
    import getopt
    from   contextlib import contextmanager

    @contextmanager
    def profile(file_stats):
        if file_stats is None:
            yield
            return

        import cProfile, pstats
        profile = cProfile.Profile()
        profile.enable()
        try:
            yield
        finally:
            profile.disable()
            with open(file_stats, 'w') as f:
                stats = pstats.Stats(profile, stream=f).sort_stats("cumulative")
                stats.print_stats(50)

    def usage():
        print(__doc__)
        exit(0)

    try:
        print('\033[35;1mFuchsia v%s\033[0m' % (__version__))
        print("Authors: %s\n" % __author__)
        mpath = tpath = profpath = fmt = None
        M = T = None
        x, y, epsilon = SR.var("x y eps")
        fmt = "m"
        seed = 0
        logger.setLevel(logging.INFO)
        kwargs, args = getopt.gnu_getopt(sys.argv[1:],
                "hvl:f:P:x:y:e:m:t:s:", ["help", "use-maple"])
        for key, value in kwargs:
            if key in ["-h", "--help"]: usage()
            if key == "-f":
                if value not in ["m", "mtx"]:
                    raise getopt.GetoptError(
                            "'%s' is not a known matrix format; "
                            "'m' and 'mtx' are supported" % value)
                fmt = value
            if key == "-l":
                fh = logging.FileHandler(value, "w")
                logger_format = '%(levelname)-8s [%(asctime)s] %(message)s'
                fh.setFormatter(ElapsedTimeFormatter(logger_format))
                logger.addHandler(fh)
            if key == "-v": logger.setLevel(logging.DEBUG)
            if key == "-P": profpath = value
            if key == "-x": x = parse(value)
            if key == "-y": y = parse(value)
            if key == "-e": epsilon = SR.var(value)
            if key == "-m": mpath = value
            if key == "-t": tpath = value
            if key == "-s": seed = int(value)
            if key == "--use-maple":
                global USE_MAPLE
                USE_MAPLE = True
        with profile(profpath):
            if len(args) == 2 and args[0] == 'fuchsify':
                M = import_matrix_from_file(args[1])
                M, t1 = simplify_by_factorization(M, x)
                M, t2 = fuchsify(M, x, seed=seed)
                T = t1*t2
            elif len(args) == 2 and args[0] == 'normalize':
                M = import_matrix_from_file(args[1])
                M, t1 = simplify_by_factorization(M, x)
                M, t2 = normalize(M, x, epsilon, seed=seed)
                T = t1*t2
            elif len(args) == 2 and args[0] == 'factorize':
                M = import_matrix_from_file(args[1])
                M, T = factorize(M, x, epsilon, seed=seed)
            elif len(args) == 2 and args[0] == 'sort':
                M = import_matrix_from_file(args[1])
                M, T, B = block_triangular_form(M)
            elif len(args) == 2 and args[0] == 'reduce':
                m = import_matrix_from_file(args[1])
                M, T = epsilon_form(m, x, epsilon, seed=seed)
            elif len(args) == 2 and args[0] == 'cat':
                M = import_matrix_from_file(args[1])
                T = None
            elif len(args) == 3 and args[0] == 'transform':
                M = import_matrix_from_file(args[1])
                t = import_matrix_from_file(args[2])
                M = transform(M, x, t)
            elif len(args) == 3 and args[0] == 'changevar':
                M = import_matrix_from_file(args[1])
                fy = parse(args[2])
                extravars = set(fy.variables()) - set(M.variables())
                if y not in extravars:
                    raise getopt.GetoptError(
                            "'%s' is missing in '%s'" % (y, fy))
                M = change_variable(M, x, y, fy)
                x = y
            elif len(args) == 2 and args[0] == 'info':
                M = import_matrix_from_file(args[1])
                print("Matrix size: %s" % M.nrows())
                print("Matrix complexity: %s" % matrix_complexity(M))
                print("Matrix expansion:")
                for point, prank in singularities(M, x).items():
                    c0 = matrix_c0(M, x, point, prank)
                    print("  C0[%s=%s; %s] = {rank=%s, complexity=%s}" % (
                        x, point, prank, c0.rank(), matrix_complexity(c0)
                    ))
                    if logger.isEnabledFor(logging.DEBUG):
                        for eigenval in sorted(set(c0.eigenvalues()), key=str):
                            print("    eigenvalue: %s" % (eigenval))
                M = None
            elif len(args) == 0:
                usage()
            else:
                raise getopt.GetoptError("unknown command: %s" % (args))
        if M is not None:
            M = partial_fraction(M, x)
            if mpath is not None:
                export_matrix_to_file(mpath, M, fmt=fmt)
            else:
                print(M)
        if tpath is not None and T is not None:
            T = partial_fraction(T, x)
            export_matrix_to_file(tpath, T, fmt=fmt)
    except getopt.GetoptError as error:
        logger.error("%s", error)
        exit(1)
    except FuchsiaError as e:
        msg, tab = str(e), '  '
        msg = tab+msg.replace('\n', '\n'+tab)
        logger.error(msg)
        exit(1)

if __name__ == '__main__':
    main()

#!/usr/bin/env sage
"""\
Usage:
    fuchsia [options] <command> <args>...

Commands:
    reduce [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>
        find an epsilon form of the given matrix

    fuchsify [-x <name>] [-m <path>] [-t <path>] <matrix>
        find a transformation that will transform a given matrix
        into Fuchsian form

    normalize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>
        find a transformation that will transform a given Fuchsian
        matrix into normalized form

    factorize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>
        find a transformation that will make a given normalized
        matrix proportional to the infinitesimal parameter

    info [-x <name>] <matrix>
        show a description of a given matrix, listing singular
        points and residue eigenvalues (in verbose more)

    sort [-m <path>] [-t <path>] <matrix>
        find a block-triangular form of the given matrix

    transform [-x <name>] [-m <path>] <matrix> <transform>
        transform a given matrix using a given transformation

    changevar [-x <name>] [-y <name>] [-m <path>] <matrix> <expr>
        transform a given matrix by susbtituting free variable
        by a given expression

Options:
    -h          show this help message
    -f <fmt>    matrix file format: mtx or m (default: m)
    -l <path>   write log to this file
    -v          produce a more verbose log
    -P <path>   save profile report into this file
    -x <name>   use this name for the free variable (default: x)
    -y <name>   use this name for the new free variable (default: y)
    -e <name>   use this name for the infinitesimal parameter (default: eps)
    -m <path>   save the resulting matrix into this file
    -t <path>   save the resulting transformation into this file
    -s <number> use this random seed when applicable (default: 0)
    --use-maple speed up calculations by using Maple when possible

Arguments:
    <matrix>    read the input matrix from this file
    <transform> read the transformation matrix from this file
    <expr>      arbitrary expression
"""

__author__ = "<NAME>, <NAME>"
__author_email__ = "<EMAIL>"
__version__ = "17.06.28"

__all__ = [
    "FuchsianSystem",
    "NormalizeAssistant",

    "balance",
    "balance_transform",
    "block_triangular_form",
    "epsilon_form",
    "export_matrix_to_file",
    "factorize",
    "fuchsify",
    "fuchsify_off_diagonal_blocks",
    "import_matrix_from_file",
    "is_fuchsian",
    "matrix_c0",
    "matrix_c1",
    "matrix_complexity",
    "matrix_residue",
    "normalize",
    "reduce_diagonal_blocks",
    "setup_fuchsia",
    "simple",
    "simplify_by_factorization",
    "simplify_by_jordanification",
    "singularities",
    "transform"
]

from   collections import defaultdict
from   functools import wraps
from   itertools import combinations, permutations
from   random import Random
import inspect
import logging
import time

from   sage.all import *
from   sage.misc.parser import Parser
from   sage.libs.ecl import ecl_eval

def logcall(f):
    @wraps(f)
    def wrapper(*args, **kwargs):
        logger.enter(f.__name__)
        r = f(*args, **kwargs)
        logger.exit(f.__name__)
        return r
    return wrapper

class ElapsedTimeFormatter(logging.Formatter):
    def __init__(self, fmt=None, datefmt=None):
        super(ElapsedTimeFormatter, self).__init__(fmt, datefmt)
        self.start_time = time.time()
    def formatTime(self, record, datefmt=None):
        if datefmt is None:
            return "%7.1fs" % (record.created - self.start_time)
        return time.strftime(datefmt, time.localtime(record.created))

class FuchsiaLogger(object):
    def __init__(self):
        log_handler = logging.StreamHandler()
        log_handler.setFormatter(ElapsedTimeFormatter(
            "\033[32m[%(asctime)s]\033[0m %(message)s"
        ))
        logger = logging.getLogger('fuchsia')
        logger.addHandler(log_handler)
        logger.setLevel(logging.WARNING)

        self.addHandler = logger.addHandler
        self.error = logger.error
        self.isEnabledFor = logger.isEnabledFor
        self.logger = logger
        self.setLevel = logger.setLevel

        self.depth = 0

    def indent(self):
        return "   " * self.depth

    def is_verbose(self):
        return self.isEnabledFor(logging.INFO)

    def enter(self, name):
        self.info("-> %s" % name)
        self.depth += 1

    def exit(self, name):
        self.depth -= 1
        self.info("<- %s" % name)

    def debug(self, msg):
        self.logger.debug(self.indent() + msg)

    def info(self, msg):
        self.logger.info(self.indent() + msg)

if True:
    ecl_eval("(ext:set-limit 'ext:heap-size 0)")
    logger = FuchsiaLogger()
    USE_MAPLE = False

def setup_fuchsia(verbosity=0, use_maple=False):
    global USE_MAPLE
    USE_MAPLE = bool(use_maple)
    logger.setLevel(
        logging.WARNING if verbosity <= 0 else \
        logging.INFO if verbosity == 1 else
        logging.DEBUG
    )

class FuchsiaError(Exception):
    pass

def cross_product(v1, v2):
    m1, m2 = matrix(v1), matrix(v2)
    return m1.transpose() * m2

def dot_product(v1, v2):
    m1, m2 = matrix(v1), matrix(v2)
    sp = m1 * m2.transpose()
    return sp[0,0]

def partial_fraction(M, var):
    return M.apply_map(lambda ex: ex.partial_fraction(var))

#@logcall
def fuchsia_simplify(obj, x=None):
    if USE_MAPLE:
        def maple_simplify(ex):
            if hasattr(ex, "variables") and ((x is None) or (x in ex.variables())):
                res = maple.factor(maple.radnormal(ex))
                return parse(str(res))
            else:
                return ex
        if hasattr(obj, "apply_map"):
            return obj.apply_map(maple_simplify)
        else:
            return maple_simplify(obj)
    else:
        return obj.simplify_rational()

def fuchsia_solve(eqs, var):
    if USE_MAPLE:
        s = maple.solve(eqs, var)
        solutions = s.parent().get(s._name).strip('[]').split('],[')
        solutions = [s.split(',') for s in solutions if s != '']
        result = []
        for solution in solutions:
            r = []
            for s in solution:
                try:
                    expr = fuchsia_simplify(parse(s))
                    r.append(expr)
                except SyntaxError as error:
                    print("ERROR:  \n%s\n  %s\n" % (s, error))
                    continue
            result.append(r)
        return result
    else:
        return solve(eqs, var, solution_dict=True)

def change_variable(m, x, y, fy):
    mm = m.subs({x: fy}) * derivative(fy, y)
    return mm

def transform(M, x, T):
    """Given a system of differential equations dF/dx = M*F,
    and a transformation of base functions F = T*F', compute
    and return M', such that dF'/dx = M'*F'.

    Note: M' = inverse(T)*(M*T - dT/dx)
    """
    mm = T.inverse()*(M*T - derivative(T, x))
    return mm

def balance(P, x1, x2, x):
    assert P.is_square()
    assert (P*P - P).is_zero()
    assert not bool(x1 == x2)
    coP = identity_matrix(P.nrows()) - P
    if x1 == oo:
        b = coP - (x - x2)*P
    elif x2 == oo:
        b = coP - 1/(x - x1)*P
    else:
        b = coP + (x - x2)/(x - x1)*P
    return b

def balance_transform(M, P, x1, x2, x):
    """Same thing as transform(M, x, balance(P, x1, x2, x)), but faster."""
    assert P.is_square()
    #assert (P*P - P).is_zero()
    assert not bool(x1 == x2)
    coP = identity_matrix(P.nrows()) - P
    if x1 == oo:
        k = -(x - x2)
        d = -1
    elif x2 == oo:
        k = -1/(x - x1)
        d = 1/(x - x1)**2
    else:
        k = (x - x2)/(x - x1)
        d = (x2 - x1)/(x - x1)**2
    mm = (coP + 1/k*P)*M*(coP + k*P) - d/k*P
    return mm

#@logcall
def limit_fixed(expr, x, x0):
    #logger.debug("[%s -> %s] %s" % (x, x0, expr))
    res = None
    if USE_MAPLE:
        res = limit_fixed_maple(expr, x, x0)
    else:
        res = limit_fixed_maxima(expr, x, x0)
    #logger.debug("res = %s" % res)
    return res

def limit_fixed_maple(expr, x, x0):
    res = maple.limit(expr, **{str(x): x0})
    try:
        res = parse(str(res))
    except:
        res = NaN
    return res

def limit_fixed_maxima(expr, x, x0):
    """Return a limit of expr when x->lim.

    The standard 'limit()' function of SageMath does not allow
    you to specify the variable, only it's name as a keyword
    argument. If you have a variable, and not it's name, use
    this function instead.
    """
    l = maxima_calculus.sr_limit(expr, x, x0)
    return expr.parent()(l)

@logcall
def singularities(m, x):
    """Find values of x around which rational matrix M has
    a singularity; return a dictionary with {val: p} entries,
    where p is the Poincare rank of M at x=val.

    Example:
    >>> x = var("x")
    >>> s = singularities(matrix([[1/x, 0], [1/(x+1)**3, 1/(x+1)]]), x)
    >>> sorted(s.items())
    [(-1, 2), (0, 0), (+Infinity, 0)]

    >>> s = singularities(matrix([[x, 1, 1/x]]), x)
    >>> sorted(s.items())
    [(0, 0), (+Infinity, 2)]
    """
    m = fuchsia_simplify(m)
    result = {}
    for expr in m.list():
        points_expr = singularities_expr(expr, x)
        for x0, p in points_expr.items():
            if x0 in result:
                result[x0] = max(result[x0], p)
            else:
                result[x0] = p
    return result

def singularities_expr(expr, x):
    if USE_MAPLE:
        return singularities_expr_maple(expr, x)
    else:
        return singularities_expr_maxima(expr, x)

def singularities_expr_maple(expr, x):
    if bool(expr == 0):
        return {}

    result = {}
    sols = fuchsia_solve(1/expr, x)
    points = [x0 for x0 in sols[0]] if len(sols) > 0 else []
    for x0 in points:
        if x0 not in result:
            result[x0] = 0
        else:
            result[x0] += 1

    sols = fuchsia_solve((1/(expr.subs({x: 1/x})/x**2)).simplify_rational(), x)
    points = [x0 for x0 in sols[0]] if len(sols) > 0 else []
    for x0 in points:
        if x0 == 0:
            if oo not in result:
                result[oo] = 0
            else:
                result[oo] += 1
    return result

def singularities_expr_maxima(expr, x):
    if expr.is_zero():
        return {}

    result = {}
    eq = factor(1/expr)
    points, ps = solve(eq, x, solution_dict=True, multiplicities=True)
    for sol, p in zip(points, ps):
        if x not in sol:
            raise Exception("Maxima can't solve this equation: %s" % eq)
        x0 = expand(sol[x])
        if p >= 1:
            result[x0] = p-1

    eq = factor(1/(expr.subs({x: 1/x})/x**2))
    points, ps = solve(eq, x, solution_dict=True, multiplicities=True)
    for sol, p in zip(points, ps):
        if x not in sol:
            raise Exception("Maxima can't solve this equation: %s" % eq)
        if sol[x] == 0 and p >= 1:
            result[oo] = p-1
    return result

def matrix_taylor0(m, x, x0, exp):
    """Return the 0-th coefficient of Taylor expansion of
    a matrix M around a finite point x=point, assuming that
    M(x->point)~1/(x-point)**exp.

    Example:
    >>> x = var('x')
    >>> matrix_taylor0(matrix([[x/(x-1), 1/x, x, 1]]), x, 1, 1)
    [1 0 0 0]
    """
    return m.apply_map(lambda ex: limit_fixed(ex*(x-x0)**exp, x, x0))

def matrix_taylor1(M, x, point, exp):
    """Return the 1-th coefficient of Taylor expansion of
    a matrix M around a finite point x=point, assuming that
    M(x->point)~1/(x-point)**exp.

    Example:
    >>> x = var('x')
    >>> matrix_taylor1(matrix([[x/(x-1), 1/x, x, 1]]), x, 0, 1)
    [0 0 0 1]
    """
    def taylor1(e, x):
        l0 = limit_fixed(e, x, 0)
        l1 = limit_fixed((e - l0)/x, x, 0)
        return l1
    return matrix([
        [taylor1(e, x) for e in row]
        for row in M.subs({x: x+point})*x**exp
    ])

def matrix_c0(M, x, point, p):
    """Return the 0-th coefficient of M's expansion at x=point,
    assuming Poincare rank of M at that point is p. If point is
    +Infinity, return minus the coefficient at the highest power of x.

    Examples:
    >>> x = var("x")
    >>> m = matrix([[1/x, 1/x**2], [1, 1/(x-1)]])
    >>> matrix_c0(m, x, 0, 1)
    [0 1]
    [0 0]
    >>> matrix_c0(m, x, 1, 0)
    [0 0]
    [0 1]
    >>> matrix_c0(m, x, oo, 1)
    [ 0  0]
    [-1  0]
    >>> matrix_c0(m*x, x, oo, 2)
    [ 0  0]
    [-1  0]
    """
    if point == oo:
        return -matrix_taylor0(M.subs({x: 1/x}), x, 0, p-1)
    else:
        return matrix_taylor0(M, x, point, p+1)

def matrix_c1(M, x, point, p):
    """Return the 1-st coefficient of M's expansion at x=point,
    assuming Poincare rank of M at that point is p. If point is
    +Infinity, return minus the coefficient at the second-to-highest
    power of x.

    Examples:
    >>> x = var("x")
    >>> m = matrix([[1/x, 1/x**2], [1, 1/(x-1)]])
    >>> matrix_c1(m, x, 0, 1)
    [1 0]
    [0 0]
    >>> matrix_c1(m, x, oo, 1)
    [-1  0]
    [ 0 -1]
    """
    if point == oo:
        return -matrix_taylor1(M.subs({x: 1/x}), x, 0, p-1)
    else:
        return matrix_taylor1(M, x, point, p+1)

def matrix_residue(M, x, x0):
    """Return matrix residue of M at x=x0, assuming that M's
    Poincare rank at x=x0 is zero.

    Example:
    >>> x = var("x")
    >>> m = matrix([[1/x, 2/x], [3/x, 4/x]])
    >>> matrix_residue(m, x, 0)
    [1 2]
    [3 4]
    >>> matrix_residue(m, x, oo)
    [-1 -2]
    [-3 -4]
    """
    if M._cache is None:
        M._cache = {}
    key = "matrix_residue_%s_%s" % (x, x0)
    if key in M._cache:
        return M._cache[key]

    m0 = matrix_c0(M, x, x0, 0)
    M._cache[key] = m0
    return m0

def matrix_is_nilpotent(M):
    """Return True if M is always nilpotent, False otherwise.

    Examples:
    >>> a, b, c = var('a b c')
    >>> matrix_is_nilpotent(matrix([[0,a,b],[0,0,c],[0,0,0]]))
    True
    >>> matrix_is_nilpotent(matrix([[a,0,0],[0,b,0],[0,0,c]]))
    False
    """
    for v in M.eigenvalues():
        if not v.is_zero():
            return False
    return True

def jordan_cell_sizes(J):
    """Return a tuple of Jordan cell sizes from a matrix J in Jordan
    normal form.

    Examples:
    >>> jordan_cell_sizes(matrix([[1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,0,1]]))
    (2, 2)
    >>> jordan_cell_sizes(matrix([[1,1,0,0],[0,1,1,0],[0,0,1,0],[0,0,0,1]]))
    (3, 1)
    >>> jordan_cell_sizes(zero_matrix(5,5))
    (1, 1, 1, 1, 1)
    """
    assert J.is_square()
    sizes = []
    n = 1
    for i in range(J.nrows() - 1):
        if J[i, i+1].is_zero():
            sizes.append(n)
            n = 1
        else:
            n += 1
    sizes.append(n)
    assert sum(sizes) == J.nrows()
    return tuple(sizes)

def solve_right_fixed(A, B):
    """As of SageMath 6.10, 'Matrix.solve_right' method uses a broken
    check for solution correctness; this function corrects that check.
    """
    C = A.solve_right(B, check=False)
    if C.nrows() and C.ncols():
        # Sometimes 'solve_right' returns a giant expression,
        # which can be simplified into a tiny one; without this
        # simplification 'is_zero' call below may take forever
        # to finish.
        #
        # Note that the condition above is there to make sure
        # that an empty matrix is not passed into 'simplify_rational',
        # otherwise you'll get this error:
        #   TypeError: unable to make sense of Maxima expression
        #   'matrix()' in Sage
        C = C.simplify_rational()
    if not (A*C - B).simplify_rational().is_zero():
        raise ValueError("matrix equation has no solutions")
    return C

def solve_left_fixed(A, B):
    """As of SageMath 6.10, 'Matrix.solve_left' method uses a broken
    check for solution correctness; this function corrects that check.
    """
    return solve_right_fixed(A.transpose(), B.transpose()).transpose()

def any_integer(rng, ring, excluded):
    r = 2
    while True:
        p = ring(rng.randint(-r, r))
        if p not in excluded:
            return p
        r *= 2

#==================================================================================================
# Transformation routines
#==================================================================================================

def block_triangular_form(m):
    """Find a lower block-triangular form of a given matrix.

    Return a tuple `(M, T, B)` where:
      * `M` is a new matrix;
      * `T` is a transformation matrix;
      * `B` is a list of tuples (ki, ni), such that M's i-th
        diagonal block is given by `M.submatrix(ki, ki, ni, ni)`.
    """
    logger.enter("block_triangular_form")
    if logger.is_verbose():
        logger.debug("matrix before transformation:\n%s\n" % matrix_mask_str(m))

    n = m.nrows()
    deps_1_1 = {}
    for i, row in enumerate(m.rows()):
        deps_1_1[i] = set(j for j,ex in enumerate(row) if not bool(ex==0))

    deps_1_all = {}
    def find_deps_1_all(i):
        if i in deps_1_all:
            return deps_1_all[i]
        deps_1_all[i] = deps_1_1[i]
        for j in deps_1_1[i]:
            if i == j:
                continue
            find_deps_1_all(j)
            deps_1_all[i] = deps_1_all[i].union(deps_1_all[j])
    [find_deps_1_all(j) for j in range(n)]

    deps_coup = dict((i, set([])) for i in range(n))
    for i in range(n):
        for j in range(n):
            if (i in deps_1_all[j]) and (j in deps_1_all[i]):
                deps_coup[i].update([i,j])
                deps_coup[j].update([i,j])

    shuffle = []
    error = False
    while not error:
        if not deps_coup:
            break
        error = True
        for i in range(n):
            if i not in deps_coup:
                continue
            b = deps_coup[i].copy()
            if b != deps_1_all[i]:
                continue
            if b == set([]):
                b = set([i])
            shuffle += [b]
            for j in deps_1_all[i].copy():
                if j in deps_1_all:
                    del deps_1_all[j]
                    del deps_coup[j]
            for j in range(n):
                if j not in deps_coup:
                    continue
                deps_coup[j].difference_update(b)
                deps_1_all[j].difference_update(b)
            if i in deps_coup:
                del deps_coup[i]
            if i in deps_1_all:
                del deps_1_all[i]
            error = False
            break
    if error:
        raise FuchsiaError("Infinite loop")

    t = zero_matrix(SR,n)
    i = 0
    blocks = []
    for block in shuffle:
        blocks.append((i, len(block)))
        for j in list(block):
            t[j,i] = 1
            i += 1
    logger.info("found %d blocks" % len(blocks))
    mt = transform(m, None, t)
    if logger.is_verbose():
        logger.debug("matrix after transformation:\n%s\n" % matrix_mask_str(mt))
    logger.exit("block_triangular_form")
    return mt, t, blocks

def epsilon_form(m, x, eps, seed=0):
    logger.enter("epsilon_form")
    m, t1, b = block_triangular_form(m)
    m, t2 = reduce_diagonal_blocks(m, x, eps, b=b, seed=seed)
    m, t3 = fuchsify_off_diagonal_blocks(m, x, eps, b=b)
    m, t4 = factorize(m, x, eps, b=b, seed=seed)
    t = t1*t2*t3*t4
    logger.exit("epsilon_form")
    return m, t

def matrix_mask(m):
    n = m.nrows()
    return matrix(SR, n, n, [int(not bool(ex==0)) for ex in m.list()])

def matrix_mask_str(m):
    s = ''
    for row in matrix_mask(m).rows():
        s += ' '.join([str(ex) for ex in row]).replace('0', '.').replace('1','x') + "\n"
    return s

def matrix_str(m, n=2):
    buf = ""
    ind = " "*n
    for col in m.columns():
        for ex in col:
            buf += ind + str(ex) + "\n"
        buf += "\n"
    buf = buf[:-1]
    return buf

#==================================================================================================
# Step I: Fuchsify
#==================================================================================================

def is_fuchsian(m, x):
    for i, expr in enumerate(m.list()):
        if expr.is_zero():
            continue
        points = singularities_expr(expr, x)
        for x0, p in points.items():
            if p > 0:
                return False
    return True

def fuchsify(M, x, seed=0):
    """Given a system of differential equations of the form dF/dx=M*F,
    try to find a transformation T, which will reduce M to Fuchsian
    form. Return the transformed M and T. Raise FuchsiaError if
    M can not be transformed into Fuchsian form.

    Note that such transformations are not unique; you can obtain
    different ones by supplying different seeds.
    """
    logger.enter("fuchsify")
    assert M.is_square()
    rng = Random(seed)
    poincare_map = singularities(M, x)
    def iter_reductions(p1, U):
        for p2, prank2 in poincare_map.items():
            if bool(p2 == p1): continue
            while prank2 >= 0:
                B0 = matrix_c0(M, x, p2, prank2)
                if not B0.is_zero(): break
                poincare_map[p2] = prank2 = prank2 - 1
            if prank2 < 0: continue
            v = find_dual_basis_spanning_left_invariant_subspace(B0, U, rng)
            if v is None: continue
            P = fuchsia_simplify(U*v, x)
            M2 = fuchsia_simplify(balance_transform(M, P, p1, p2, x), x)
            yield p2, P, M2
    combinedT = identity_matrix(M.base_ring(), M.nrows())
    reduction_points = [pt for pt,p in poincare_map.items() if p >= 1]
    reduction_points.sort()
    if reduction_points == []:
        logger.info("already fuchsian")
    else:
        for pt in reduction_points:
            logger.info("rank = %d, x = %s" % (poincare_map[pt], pt))
    while reduction_points:
        pointidx = rng.randint(0, len(reduction_points) - 1)
        point = reduction_points[pointidx]
        prank = poincare_map[point]
        if prank < 1:
            del reduction_points[pointidx]
            continue
        while True:
            A0 = matrix_c0(M, x, point, prank)
            if A0.is_zero(): break
            A1 = matrix_c1(M, x, point, prank)
            try:
                U, V = alg1x(A0, A1, x)
            except FuchsiaError as e:
                logger.debug("Managed to fuchsify matrix to this state:\n"
                        "%s\nfurther reduction is pointless:\n%s" % (M, e))
                raise FuchsiaError("matrix cannot be reduced to Fuchsian form")
            try:
                point2, P, M = min(iter_reductions(point, U), \
                        key=lambda point2_P_M2: matrix_complexity(point2_P_M2[2]))
            except ValueError as e:
                point2 = any_integer(rng, M.base_ring(), poincare_map)
                P = fuchsia_simplify(U*V, x)
                M = balance_transform(M, P, point, point2, x)
                M = fuchsia_simplify(M, x)
                logger.info("Will introduce an apparent singularity at %s." % point2)
            logger.debug(
                "Applying balance between %s and %s with projector:\n%s" % (point, point2, P))
            combinedT = combinedT * balance(P, point, point2, x)
            if point2 not in poincare_map:
                poincare_map[point2] = 1
        poincare_map[point] = prank - 1
    combinedT = fuchsia_simplify(combinedT, x)
    logger.exit("fuchsify")
    return M, combinedT

def fuchsify_off_diagonal_blocks(m, x, eps, b=None):
    logger.enter("fuchsify_off_diagonal_blocks")
    n = m.nrows()
    if b is None:
        m, t, b = block_triangular_form(m)
    else:
        t = identity_matrix(SR, n)
    for i, (ki, ni) in enumerate(b):
        for j, (kj, nj) in enumerate(reversed(b[:i])):
            pts = singularities(m.submatrix(ki, kj, ni, nj), x)
            printed = False
            while any(pts.values()):
                bj = m.submatrix(ki, kj, ni, nj)
                if bj.is_zero():
                    break
                for x0, p in pts.items():
                    if p < 1:
                        continue
                    if not printed:
                        printed = True
                        logger.info("Fuchsifying block (%d, %d) (%d, %d)" % (ki, ni, kj, nj))
                        logger.debug("  singular points = %s" % (pts,))
                    a0 = matrix_residue(m.submatrix(ki, ki, ni, ni)/eps, x, x0)
                    b0 = matrix_c0(bj, x, x0, p)
                    c0 = matrix_residue(m.submatrix(kj, kj, nj, nj)/eps, x, x0)

                    d_vars = [gensym() for i in range(ni*nj)]
                    d = matrix(SR, ni, nj, d_vars)
                    eq = d + eps/p*(a0*d - d*c0) + b0/p
                    sol = solve(eq.list(), *d_vars, solution_dict=True)
                    d = d.subs(sol[0])

                    t0 = identity_matrix(SR, n)
                    t0[ki:ki+ni, kj:kj+nj] = \
                            d/(x-x0)**p if not (x0 == oo) else d*(x**p)
                    m = fuchsia_simplify(transform(m, x, t0), x)

                    t = fuchsia_simplify(t*t0, x)
                    pts[x0] -= 1
    logger.exit("fuchsify_off_diagonal_blocks")
    return m, t

def reduce_at_one_point(M, x, v, p, v2=oo):
    """Given a system of differential equations of the form dF/dx=M*F,
    with M having a singularity around x=v with Poincare rank p,
    try to find a transformation T, which will reduce p by 1 (but
    possibly introduce another singularity at x=v2). Return the
    transformed M and T.
    """
    assert M.is_square()
    assert p > 0
    n = M.nrows()
    combinedT = identity_matrix(M.base_ring(), n)
    while True:
        A0 = matrix_c0(M, x, v, p)
        if A0.is_zero(): break
        A1 = matrix_c1(M, x, v, p)
        U, V = alg1x(A0, A1, x)
        P = U*V
        M = balance_transform(M, P, v, v2, x)
        M = fuchsia_simplify(M, x)
        combinedT = combinedT * balance(P, v, v2, x)
    combinedT = fuchsia_simplify(combinedT, x)
    return M, combinedT

def find_dual_basis_spanning_left_invariant_subspace(A, U, rng):
    """Find matrix V, such that it's rows span a left invariant
    subspace of A and form a dual basis with columns of U (that
    is, V*U=I).
    """
    evlist = []
    for eigenval, eigenvects, evmult in A.eigenvectors_left():
        evlist.extend(eigenvects)
    rng.shuffle(evlist)
    W = matrix(evlist)
    try:
        M = solve_left_fixed(W*U, identity_matrix(U.ncols()))
        return M*W
    except ValueError:
        return None

def alg1x(A0, A1, x):
    #if not matrix_is_nilpotent(A0):
    #    raise FuchsiaError("matrix is irreducible (non-nilpotent residue)")
    # We will rely on undocumented behavior of jordan_form() call:
    # we will take it upon faith that it sorts Jordan cells by
    # their size, and puts the largest ones closer to (0, 0).
    A0J, U = A0.jordan_form(transformation=True)
    invU = U.inverse()
    A0J_cs = jordan_cell_sizes(A0J)
    assert all(A0J_cs[i] >= A0J_cs[i+1] for i in range(len(A0J_cs) - 1))
    ncells = len(A0J_cs)
    A0J_css = [sum(A0J_cs[:i]) for i in range(ncells + 1)]
    nsimplecells = sum(1 if s == 1 else 0 for s in A0J_cs)
    u0 = [U[:,A0J_css[i]] for i in range(ncells)]
    v0t = [invU[A0J_css[i+1]-1,:] for i in range(ncells)]
    L0 = matrix([
        [(v0t[k]*(A1)*u0[l])[0,0] for l in range(ncells)]
        for k in range(ncells)
    ])
    L0 = fuchsia_simplify(L0, x)
    L1 = matrix([
        [(v0t[k]*u0[l])[0,0] for l in range(ncells)]
        for k in range(ncells)
    ])
    assert (L1 - diagonal_matrix(
            [0]*(ncells-nsimplecells) + [1]*nsimplecells)).is_zero()
    #zero_rows = [i for i in range(A0.nrows()) if A0J[i,:].is_zero()]
    #zero_cols = [j for j in range(A0.nrows()) if A0J[:,j].is_zero()]
    #assert len(zero_rows) == len(zero_cols) == ncells
    #_L0 = (invU*A1*U)[zero_rows, zero_cols]
    #assert (L0 - _L0).is_zero()
    lam = SR.symbol()
    if not (L0 - lam*L1).determinant().is_zero():
        raise FuchsiaError("matrix is Moser-irreducible")
    S, D = alg1(L0, A0J_cs)
    I_E = identity_matrix(D.base_ring(), A0.nrows())
    for i in range(ncells):
        for j in range(ncells):
            if not D[i,j].is_zero():
                ni = A0J_css[i]
                nj = A0J_css[j]
                for k in range(min(A0J_cs[i], A0J_cs[j])):
                    I_E[ni+k,nj+k] += D[i,j]
    U_t = U*I_E
    invU_t = U_t.inverse()
    return \
        U_t[:, [A0J_css[i] for i in S]], \
        invU_t[[A0J_css[i] for i in S], :]

def alg1(L0, jordan_cellsizes):
    assert L0.is_square()
    ring = L0.base_ring()
    N = L0.nrows()
    I_N = identity_matrix(N)
    S = set()
    D = zero_matrix(ring, N)
    K = sum(1 if s > 1 else 0 for s in jordan_cellsizes)
    while True:
        Lx = copy(L0)
        for i in S:
            Lx[:, i] = 0
            Lx[i, :] = 0
        fi = None
        c = None
        for i in range(0, N):
            if i not in S:
                try:
                    c = solve_right_fixed(Lx[:,0:i], Lx[:,i])
                except ValueError:
                    # No solution found; vectors are independent.
                    continue
                fi = i
                break
        assert fi is not None
        assert c is not None
        D0 = matrix(ring, N)
        invD0 = matrix(ring, N)
        for j in range(fi):
            D0[j, fi] = -c[j, 0]
            invD0[j, fi] = -c[j, 0] \
                    if jordan_cellsizes[j] == jordan_cellsizes[fi] else 0
        L0 = (I_N - invD0)*L0*(I_N + D0)
        D = D + D0 + D*D0
        if fi < K:
            break
        S.add(fi)
    # Check Alg.1 promise
    for j in range(N):
        for k in range(N):
            if (j not in S) and ((k in S) or k == fi):
                assert L0[j,k] == 0
    S.add(fi)
    return S, D

#==================================================================================================
# Step II: Normalize
#==================================================================================================

def is_normalized(M, x, eps):
    """Return True if (a Fuchsian) matrix M is normalized, that
    is all the eigenvalues of it's residues in x lie in [-1/2, 1/2)
    range (in limit eps->0). Return False otherwise.

    Examples:
    >>> x, e = var("x epsilon")
    >>> is_normalized(matrix([[(1+e)/3/x, 0], [0, e/x]]), x, e)
    True
    """
    points = singularities(M, x)
    for x0, p in points.items():
        M0 = matrix_residue(M, x, x0)
        for ev in M0.eigenvalues():
            ev = limit_fixed(ev, eps, 0)
            if not (Rational((-1, 2)) <= ev and ev < Rational((1, 2))):
                return False
    return True

def are_diagonal_blocks_reduced(m, b, x, eps):
    """Return True if diagonal blocks of `m` are normalized, that is all the eigenvalues of theirs
    residues in `x` lie in the range [-1/2, 1/2) in eps->0 limit; return False otherwise. Diagonal
    blocks are defined by the list `b` which corresponds to the equivalent value returned by the
    `block_triangular_form` routine.

    Examples:
    >>> x, e = var("x epsilon")
    >>> are_diagonal_blocks_reduced(matrix([[(1+e)/3/x, 0], [0, e/x]]), [(0,1),(1,1)], x, e)
    True
    """
    for ki, ni in b:
        mi = fuchsia_simplify(m.submatrix(ki, ki, ni, ni), x)
        if not is_normalized(mi, x, eps):
            return False
    return True

def reduce_diagonal_blocks(m, x, eps, b=None, seed=0):
    """Given a lower block-triangular system of differential equations of the form dF/dx=m*F,
    find a transformation that will shift all eigenvalues of all residues of all its diagonal
    blocks into the range [-1/2, 1/2) in eps->0 limit. Return the transformed matrix m and the
    transformation; raise FuchsiaError if such transformation is not found. Diagonal blocks
     are defined by the list `b` which corresponds to the equivalent value returned by the
    `block_triangular_form` routine.
    """
    logger.enter("reduce_diagonal_blocks")
    n = m.nrows()
    if b is None:
        m, t, b = block_triangular_form(m)
    else:
        t = identity_matrix(SR, n)
    for i, (ki, ni) in enumerate(b):
        mi = fuchsia_simplify(m.submatrix(ki, ki, ni, ni), x)
        ti = identity_matrix(SR, ni)
        logger.info("reducing block #%d (%d,%d)" % (i,ki,ni))
        if logger.is_verbose():
            logger.debug("\n%s" % matrix_str(mi, 2))

        mi_fuchs, ti_fuchs = fuchsify(mi, x, seed=seed)
        ti = ti*ti_fuchs

        mi_norm, ti_norm = normalize(mi_fuchs, x, eps, seed=seed)
        ti = ti*ti_norm

        mi_eps, ti_eps = factorize(mi_norm, x, eps, seed=seed)
        ti = ti*ti_eps

        t[ki:ki+ni, ki:ki+ni] = ti
    mt = fuchsia_simplify(transform(m, x, t), x)
    logger.exit("reduce_diagonal_blocks")
    return mt, t


def cmap_add_div(Cmap, C, pi, ki, x0):
    # M += C*(x-pi)^ki/(x-x0)
    if pi == x0:
        Cmap[pi, ki-1] += C
    elif ki >= 0:
        assert pi == 0
        for k in range(0, ki):
            Cmap[pi, k] += C * x0**(ki-1-k)
        Cmap[x0, -1] += C * x0**ki
    else:
        for k in range(ki, 0):
            Cmap[pi, k] -= C * (x0-pi)**(ki-1-k)
        Cmap[x0, -1] += C * (x0-pi)**ki

def cmap_add_mul(Cmap, C, pi, ki, x0):
    # M += C*(x-pi)^ki*(x-x0) = C*(x-pi)^ki*{(x-pi) + (pi-x0)}
    Cmap[pi, ki] += C*(pi-x0)
    if ki == -1:
        Cmap[SR(0), ki+1] += C
    else:
        Cmap[pi, ki+1] += C

def memoize(f):
    @wraps(f)
    def wrapper(self):
        key = f.__name__+"_cache"
        if not hasattr(self, key):
#            print("cache miss: %s" % key)
            setattr(self, key, f(self))
#        else:
#            print("cache hit: %s" % key)
        return getattr(self, key)
    return wrapper

class FuchsianSystem(object):
    @staticmethod
    def from_M(m, x, eps):
        # singular points
        xs = sorted(singularities(m, x).keys())
        # residues in singular points
        rs = {}
        for xi in xs:
            rs[xi] = SystemResidue(matrix_residue(m, x, xi), eps)
        # TODO: check if system's residues have half-integer eigenvalues.
        # If yes, then inform user about that and let us remove resonances only.
        # Though, such a form does not allow to construct an eps-form, it
        # however can be used to find solutions as a power series using
        # semi-analytical methods.
        return FuchsianSystem(rs, x, eps)

    def __init__(self, rs, x, eps):
        self.rs = rs
        self.xs = sorted(rs.keys())
        self.x = x
        self.eps = eps

    def __str__(self):
        #s = "zeroes:\n  %s\n" % (self.xs,)
        s = ""
        a = 0
        for xi, ri in self.rs.items():
            s += "a = %s (x = %s)\n%s\n" % (a, xi, ri)
            a += 1
        return s

    def write(self, filename):
        m = SR(0)
        for xi, ri in self.rs.items():
            if xi == oo:
                continue
            m += ri.m/(self.x-xi)
        export_matrix_to_file(filename, m)

    def _balance(self, i1, j1, i2, j2):
        x1, x2 = self.xs[i1], self.xs[i2]
        r1, r2 = self.rs[x1], self.rs[x2]
        v1, v2 = r1.eigenvectors_right()[j1][1], r2.eigenvectors_left()[j2][1]
        print("x1 = %s" % x1)
        print("x2 = %s" % x2)
        print("v1 = %s" % v1)
        print("v2 = %s" % v2)

        sp = simple(dot_product(v1, v2))
        print("sp = %s" % sp)
        P = simple(cross_product(v1, v2) / sp)
        print("P  =\n%s" % P)
        m = balance_transform(self.m, P, x1, x2, self.x)
        return FuchsianSystem.from_M(m, self.x, self.eps)

    #@logcall
    def balance(self, i1, j1, i2, j2):
        x1, x2 = self.xs[i1], self.xs[i2]
        r1, r2 = self.rs[x1], self.rs[x2]
        v1, v2 = r1.eigenvectors_right()[j1][1], r2.eigenvectors_left()[j2][1]

        sp = simple(dot_product(v1, v2))
        print("sp = %s" % sp)
        P = simple(cross_product(v1, v2) / sp)
        print("P  =\n%s" % P)

        coP = 1-P
        if x1 == oo:
            Cmap = defaultdict(SR.zero)
            for pi, Ci in self.rs.items():
                if pi == oo:
                    continue
                # coP Ci coP (x-pi)^-1 + P Ci P (x-pi)^-1
                Cmap[pi, -1] += coP*Ci.m*coP + P*Ci.m*P
                # coP Ci P -(x-x2) (x-pi)^-1
                cmap_add_mul(Cmap, -coP*Ci.m*P, pi, -1, x2)
                # P Ci coP -1/(x-x2) (x-pi)^-1
                cmap_add_div(Cmap, -P*Ci.m*coP, pi, -1, x2)
            # -P/(x-x2)
            Cmap[x2, -1] -= P
        elif x2 == oo:
            Cmap = defaultdict(SR.zero)
            for pi, Ci in self.rs.items():
                if pi == oo:
                    continue
                # coP Ci coP (x-pi)^-1 + P Ci P (x-pi)^-1
                Cmap[pi, -1] += coP*Ci.m*coP + P*Ci.m*P
                # P Ci coP -(x-x1) (x-pi)^-1
                cmap_add_mul(Cmap, -P*Ci.m*coP, pi, -1, x1)
                # coP Ci P -1/(x-x1) (x-pi)^-1
                cmap_add_div(Cmap, -coP*Ci.m*P, pi, -1, x1)
            # P/(x-x1)
            Cmap[x1, -1] += P
        else:
            Cmap = defaultdict(SR.zero, [((pi,-1),ri.m) for pi,ri in self.rs.items()])
            for pi, Ci in self.rs.items():
                if pi == oo:
                    continue
                # coP Ci P (x1-x2)/(x-x1) (x-pi)^-1
                cmap_add_div(Cmap, coP*Ci.m*P*(x1-x2), pi, -1, x1)
                # P Ci coP (x2-x1)/(x-x2) (x-pi)^-1
                cmap_add_div(Cmap, P*Ci.m*coP*(x2-x1), pi, -1, x2)
            # P/(x-x1) - P/(x-x2)
            Cmap[x1, -1] += P
            Cmap[x2, -1] -= P
        rs = {}
        Coo = SR.zero()
        for key, C in Cmap.items():
            pi, ki = key
            if pi == oo:
                continue
            Cmap[key] = C = simple(C)
            if C.is_zero():
                del Cmap[key]
            else:
                assert(ki == -1)
                rs[pi] = SystemResidue(C)
                Coo -= C
        Coo = simple(Coo)
        rs[oo] = SystemResidue(Coo)
        #T = fuchsia_simplify(self.T * balance(P, x1, x2, self.x))
        #return RationalSystem(dict(Cmap), self.x, T, self.trace + [("balance", P, x1, x2)])
        return FuchsianSystem(rs, self.x, self.eps) 

    #@logcall
    @memoize
    def balances(self):
        res = []
        for i1,x1 in enumerate(self.xs):
            for i2,x2 in enumerate(self.xs):
                if i1 == i2:
                    continue
                for j1, j2, aux in self.balances2(x1, x2):
                    res.append((i1,j1, i2,j2, aux))
        return res

    #@logcall
    def balances2(self, x1, x2):
        #assert(not (x1-x2).is_zero())
        res = []
        r1, r2 = self.rs[x1], self.rs[x2]
        for j1, (a1,v1) in enumerate(r1.eigenvectors_right()):
            for j2, (a2,v2) in enumerate(r2.eigenvectors_left()):
                sp = simple(dot_product(v1, v2))
                if sp.is_zero():
                    continue
                res.append((j1, j2, (a1,a2, sp)))
        return res

    #@logcall
    @memoize
    def eigenvalues(self):
        res = []
        for xi in self.xs:
            res.append((xi, self.rs[xi].eigenvalues()))
        return res

#    @memoize
#    @logcall
#    def eigenvectors_left(self):
#        return self._eigenvectors(left=True)
#
#    @memoize
#    @logcall
#    def eigenvectors_right(self):
#        return self._eigenvectors(left=False)
#
#    def _eigenvectors(self, left=True):
#        res = []
#        for xi, ri in self.rs.items():
#            eigenvectors = ri.eigenvectors_left
#            if left is not True:
#                eigenvectors = ri.eigenvectors_right
#            res.append((xi, eigenvectors()))
#        return res

    def normalize(self):
        pass

class SystemResidue(object):
    def __init__(self, m, eps=None):
        self.m = m
        self.eps = eps

    def __str__(self):
        return str(self.m)

    def __repr__(self):
        return self.m.__repr__()

    @logcall
    @memoize
    def eigenvalues(self):
        return simple(mathematica.Eigenvalues(self.m).sage())

    @logcall
    @memoize
    def eigenvectors_left(self):
        return self._eigenvectors(left=True)

    @logcall
    @memoize
    def eigenvectors_right(self):
        return self._eigenvectors(left=False)

    def _eigenvectors(self, left=True):
        m = self.m.transpose() if left is True else self.m
        s = simple(mathematica.Eigensystem(m).sage())
        res = zip(s[0], s[1])
        return res

#    def _eigenvectors(self, left=True):
#        res = []
#        eigenvectors = self.m.eigenvectors_left
#        if left is not True:
#            eigenvectors = self.m.eigenvectors_right
#        res = [vvm[0:2] for vvm in eigenvectors()]
#        return res

# helpers
@logcall
def simple(obj, x=None):
    def _simplify(ex):
            return mathematica.Simplify(ex).sage()
    if hasattr(obj, "apply_map"):
        return obj.apply_map(_simplify)
    else:
        return _simplify(obj)

def ex_in(ex, lst):
    #if any((ex-t).is_zero() for t in printed):
    for t in lst:
        if (ex-t).is_zero():
            return True
    return False
        

class NormalizeAssistant(object):
    def __init__(self, m):
        self.m = m           # current matrix to balance
        self.history = []    # matrices balanced so far

    # Entry point
    def start(self):
#        cmd = self.cmd_main
#        while cmd != None:
#            cmd = cmd()
#        return
        try:
            cmd = self.cmd_main
            while cmd != None:
                cmd = cmd()
        except Exception as error:
            print("UNEXPECTED ERROR:\n%s" % (error,))
            return self.start()

    # Helpers
    def choose_balance(self, i1, i2):
        x1, x2 = self.m.xs[i1], self.m.xs[i2]
        bs = self.m.balances2(x1, x2)
        for i, (j1, j2, aux) in enumerate(bs):
            print("[%s] %s" % (i, aux))
        n = raw_input("(choose balance) > ")
        try:
            n = int(n)
        except ValueError:
            return None
        if n < 0 or n > len(bs):
            return None
        b = bs[n]
        return (i1, b[0], i2, b[1])

    # Commands
    def cmd_balance(self, b):
        print("balancing with %s" % (b,))

        mm = apply(self.m.balance, b)
        self.history.append(mm)
        self.m = mm
        return self.cmd_main

    def cmd_balance_res(self, args):
        if len(args) != 2:
            return self.cmd_error("needs exactly 2 arguments")
        try:
            i1, i2 = [int(n) for n in args]
        except ValueError:
            return self.cmd_error("arguments out of range (type e for valid values)")
        b = self.choose_balance(i1, i2)
        if b is None:
            return self.cmd_error("unknown balance")
        return self.cmd_balance(b[0:4])

    def cmd_error(self, msg):
        print("ERROR: %s" % msg)
        return self.cmd_main

    def cmd_main(self):
        cmd = raw_input("> ").split(' ')
        cmd, args = cmd[0], cmd[1:]
        if cmd == 's':
            return self.cmd_save(args)
        elif cmd == 'b':
            return self.cmd_balance_res(args)
        elif cmd == 'e':
            return self.cmd_print_eigenvalues
        elif cmd == 'h':
            return self.cmd_print_help
        elif cmd == 'q':
            return None
        elif cmd == '':
            return self.cmd_print_help
        return self.cmd_unknown_command(cmd)

    def cmd_print_eigenvalues(self):
        for i, (xi, a) in enumerate(self.m.eigenvalues()):
            print("[%s] x = %s\n    %s" % (i, xi, a))
        return self.cmd_main

    def cmd_print_help(self):
        """
        Available commands:
            b <i> <j>    balance xi, xj points
            e            show eigenvectors
            h            print this help
            q            quit
            s <file>     save matrix to file
        """
        print(inspect.cleandoc(self.cmd_print_help.__doc__))
        return self.cmd_main

    def cmd_save(self, args):
        if len(args) != 1:
            return self.cmd_error("needs exactly 1 arguments")
        name = args[0]
        self.m.write(name)
        return self.cmd_main

    def cmd_unknown_balance(self, n):
        print("Unknown balance: '%s'" % n)
        return self.cmd_main

    def cmd_unknown_command(self, cmd):
        print("Unknown command: '%s'" % cmd)
        return self.cmd_print_help

def normalize(m, x, eps, seed=0):
    """Given a Fuchsian system of differential equations of the
    form dF/dx=m*F, find a transformation that will shift all
    the eigenvalues of m's residues into [-1/2, 1/2) range (in
    the limit eps->0). Return the transformed matrix m and the
    transformation. Raise FuchsiaError if such transformation
    is not found.
    """
    class State(object):
        def __init__(self, m, x, eps, seed):
            # random
            self.random = Random(seed)
            # points
            self.points = singularities(m, x).keys()
            # ev_cum
            self.ev_cum = {}
            for x0 in self.points:
                pos, neg = 0, 0
                a0 = matrix_residue(m, x, x0)
                for ev in a0.eigenvalues():
                    ev0 = limit_fixed(ev, eps, 0)
                    if ev0 > 0:
                        pos += ev0
                    if ev0 < 0:
                        neg += ev0
                    logger.debug("x = %s, ev = %s" % (x0, ev0))
                self.ev_cum[x0] = [pos,neg]
            # x0
            logger.debug("state.ev_cum = %s" % self.ev_cum)
            self.x0 = None

        def is_normalized(self):
            for ev in self.ev_cum.values():
                if ev != [0,0]:
                    return False
            return True

        def pairs(self):
            points = [x0 for x0 in self.points if self.ev_cum[x0] != [0,0]]
            if len(points) == 1:
                self.select_x0(points[0])
                points.append(self.x0)
            return permutations(points, 2)

        def select_x0(self, x1):
            if self.x0 is not None:
                return
            for x0 in self.points:
                if x0 != x1:
                    self.x0 = x0
                    return

        def update_ev_cum(self, x2, x1):
            ev_cum_x1, ev_cum_x2 = self.ev_cum[x1], self.ev_cum[x2]
            if ev_cum_x1[0] > 0:
                ev_cum_x1[0] -= 1
            else:
                ev_cum_x1[1] -= 1
            if ev_cum_x2[1] < 0:
                ev_cum_x2[1] += 1
            else:
                ev_cum_x2[0] += 1

    logger.enter("normalize")
    state = State(m, x, eps, seed)
    T = identity_matrix(m.base_ring(), m.nrows())
    if state.is_normalized():
        logger.info("already normalized")
    i = 0
    while not state.is_normalized():
        i += 1
        m = fuchsia_simplify(m, x)
        logger.info("step %s" % i)
        balances = find_balances(m, x, eps, state)
        b = select_balance(balances, eps, state)
        if b is None:
            raise FuchsiaError("can not balance matrix")
        logger.info("  balancing x = %s and x = %s" % (b[1],b[2]))
        if logger.is_verbose():
            logger.debug("\n      use the balance:\n        %s\n" % b)

        cond, x1, x2, a0_eval, b0_eval, a0_evec, b0_evec, scale = b
        if cond == 1:
            P = cross_product(a0_evec, b0_evec) / scale
            m = balance_transform(m, P, x1, x2, x)
            T0 = balance(P, x1, x2, x)
            state.update_ev_cum(x1, x2)
        else:
            P = cross_product(b0_evec, a0_evec) / scale
            m = balance_transform(m, P, x2, x1, x)
            T0 = balance(P, x2, x1, x)
            state.update_ev_cum(x2, x1)

        T = fuchsia_simplify(T*T0, x)
    logger.exit("normalize")
    return m, T

def find_balances(m, x, eps, state={}):
    residues = {}
    for x1, x2 in state.pairs():
        logger.debug("trying to balance x = %s and x = %s" % (x1,x2))
        for xi in [x1,x2]:
            if xi not in residues:
                residues[xi] = matrix_residue(m, x, xi)
        a0, b0 = residues[x1], residues[x2]

        a0_evr, b0_evl = eigenvectors_right(a0), eigenvectors_left(b0)

        if logger.is_verbose():
            msg = "\n  Eigenvalues:\n"
            msg += "    x = %s:\n" % x1
            a0_evals = [];
            for ev, evec, emult in a0_evr:
                a0_evals += [ev]*emult
            msg += "        %s\n" % str(a0_evals)
            msg += "    x = %s:\n" % x2
            b0_evals = [];
            for ev, evec, emult in b0_evl:
                b0_evals += [ev]*emult
            msg += "        %s\n" % str(b0_evals)
            logger.debug(msg)

        balances_1 = find_balances_by_cond(a0_evr, b0_evl, lambda a0_eval, b0_eval: limit_fixed(a0_eval, eps, 0) < -0.5)
        for balance in balances_1:
            balance = [1, x1, x2] + balance
            yield balance

        a0_evl, b0_evr = eigenvectors_left(a0), eigenvectors_right(b0)
        balances_2 = find_balances_by_cond(a0_evl, b0_evr, lambda a0_eval, b0_eval: limit_fixed(a0_eval, eps, 0) >= 0.5)
        for balance in balances_2:
            balance = [2, x1, x2] + balance
            yield balance

def find_balances_by_cond(a0_ev, b0_ev, cond):
    res = []
    for a0_eval, a0_evecs, a0_evmult in a0_ev:
        for b0_eval, b0_evecs, b0_evmult in b0_ev:
            if not cond(a0_eval, b0_eval):
                logger.debug("Balance rejected:\n    a0_eval = %s\n    b0_eval = %s" % (a0_eval, b0_eval))
                continue
            for a0_evec in a0_evecs:
                for b0_evec in b0_evecs:
                    scale = fuchsia_simplify(dot_product(a0_evec, b0_evec))
                    balance = [a0_eval, b0_eval, a0_evec, b0_evec, scale]
                    if scale == 0:
                        logger.debug("Balance rejected:\n    a0_eval = %s\n    b0_eval = %s\n    a0_evec = %s\n    b0_evec = %s\n    scale   = %s" % tuple(balance))
                        continue
                    logger.debug("Balance found:\n    a0_eval = %s\n    b0_eval = %s\n    a0_evec = %s\n    b0_evec = %s\n    scale   = %s" % tuple(balance))
                    res.append(balance)
    return res

def select_balance(balances, eps, state={}):
    min_degree, min_balance = None, None
    bs = []
    for b in balances:
        cond, x1, x2, a0_eval, b0_eval, a0_evec, b0_evec, scale = b
        if (cond == 1) and limit_fixed(a0_eval, eps, 0) < -0.5 and \
                limit_fixed(b0_eval, eps, 0) >= 0.5:
            degree = max(scale.numerator().degree(eps), scale.denominator().degree(eps))
            if degree < 4:
                return b
            if (min_degree is None) or (min_degree > degree):
                 min_degree = degree
                 min_balance = b
        elif (cond == 2) and limit_fixed(a0_eval, eps, 0) >= 0.5 and \
                limit_fixed(b0_eval, eps, 0) < -0.5:
            degree = max(scale.numerator().degree(eps), scale.denominator().degree(eps))
            if degree < 4:
                return b
            if (min_degree is None) or (min_degree > degree):
                 min_degree = degree
                 min_balance = b
        bs.append(b)
    if min_balance is not None:
        return min_balance

    x0 = state.x0
    if x0 is None:
        for b in bs:
            cond, x1, x2, ev1, ev2 = b[:5]
            if cond == 1:
                x0 = x2
                break
            if cond == 2:
                x0 = x1
                break
        logger.info("  select x0 = %s" % x0)
        state.x0 = x0

    balances_x0 = [b for b in bs if (b[0] == 1 and b[2] == x0) or (b[0] == 2 and b[1] == x0)]
    b = state.random.choice(balances_x0) if balances_x0 else None
    return b

def eigenvectors_left(m):
    if m._cache is None:
        m._cache = {}
    if "eigenvectors_left" not in m._cache:
        res = simplify(m.eigenvectors_left())
        m._cache["eigenvectors_left"] = res
    return m._cache["eigenvectors_left"]

def eigenvectors_right(m):
    if m._cache is None:
        m._cache = {}
    if "eigenvectors_right" not in m._cache:
        res = m.eigenvectors_right()
        m._cache["eigenvectors_right"] = res
    return m._cache["eigenvectors_right"]

#==================================================================================================
# Step III: Factorize
#==================================================================================================

def gensym():
    sym = SR.symbol()
    SR.symbols[str(sym)] = sym
    return sym

def factorize(M, x, epsilon, b=None, seed=0):
    """Given a normalized Fuchsian system of differential equations:
        dF/dx = M(x,epsilon)*F,
    try to find a transformation that will factor out an epsilon
    from M. Return a transformed M (proportional to epsilon)
    and T. Raise FuchsiaError if epsilon can not be factored.
    """
    logger.info("-> factorize")
    n = M.nrows()
    M = fuchsia_simplify(M, x)
    if epsilon not in expand(M/epsilon).variables():
        logger.info("   already in epsilon form")
        logger.info("<- factorize")
        return M, identity_matrix(SR, n)
    rng = Random(seed)
    mu = gensym()
    if b is None:
        T_symbols = [gensym() for i in range(n*n)]
        T = matrix(SR, n, n, T_symbols)
    else:
        T, T_symbols = matrix(SR, n), []
        for ki,ni in b:
            for i in range(ki,ki+ni):
                for j in range(ki+ni):
                    sym = gensym()
                    T[i,j] = sym
                    T_symbols.append(sym)
    eqs = []
    for point, prank in singularities(M, x).items():
        assert prank == 0
        logger.debug("    processing point x = %s" % point)
        R = matrix_c0(M, x, point, 0)
        R = fuchsia_simplify(R)
        eq = (R/epsilon)*T-T*(R.subs({epsilon: mu})/mu)
        eq = fuchsia_simplify(eq)
        eqs.extend(eq.list())
    logger.info("    found %d equations with %d unknowns" % (len(eqs), len(T_symbols)))
    solutions = fuchsia_solve(eqs, T_symbols)
    for solution in solutions:
        S = T.subs(solution)
        # Right now S likely has a number of free variables in
        # it; we can set them to arbibtrary values, as long as
        # it'll make S invertible.
        rndrange = 0
        while True:
            try:
                sT = S.subs([
                    e==rng.randint(-rndrange, rndrange)
                    for e in S.variables() if e != epsilon
                ])
                sT = fuchsia_simplify(sT,x)
                M = fuchsia_simplify(transform(M, x, sT), x)
                # We're leaking a bunch of temprary variables here,
                # which accumulate in SR.variables, but who cares?
                logger.info("<- factorize")
                return M, sT
            except (ZeroDivisionError, ValueError):
                rndrange += 1 + rndrange//16
                # We've tried a bunch of substituions, and they didn't
                # work. Is the matrix at all invertible? Let's check.
                if rndrange == 16 and not S.is_invertible():
                    break
    raise FuchsiaError("can not factor epsilon")

#==================================================================================================
# Helpers
#==================================================================================================

def matrix_complexity(M):
    return len(str(M.list()))

def simplify_by_jordanification(M, x):
    """Try to simplify matrix M by constant transformations that
    transform M's residues into their Jordan forms. Return the
    simplified matrix and the transformation. If none of the
    attempted transformations reduce M's complexity (as measured
    by 'matrix_complexity()'), return the original matrix and
    the identity transformation.
    """
    minM = M
    minC = matrix_complexity(M)
    minT = identity_matrix(M.base_ring(), M.nrows())
    for point, prank in singularities(M, x).items():
        R = matrix_c0(M, x, point, prank)
        J, T = R.jordan_form(transformation=True)
        MM = fuchsia_simplify(transform(M, x, T), x)
        C = matrix_complexity(MM)
        if C < minC:
            minM = MM
            minC = C
            minT = T
    return minM, minT

def common_factor(expressions, filter):
    """Factorize given expressions, select those factors for
    which 'filter(factor)' is True, and return the product of
    factors common to all the expressions.

    Examples:
    >>> x = var("x")
    >>> common_factor([x*x-1, x+1], lambda f: True)
    x + 1
    >>> common_factor([1/x**2, 2/x**3, 3/x**4], lambda f: True)
    x^(-2)

    Note that if there is a mix of positive and negative exponents
    of a given factor, this function will use (one of) the most
    frequently occurring exponent:
    >>> common_factor([x, 1/x, 2/x**2, 3/x], lambda f: True)
    1/x
    """
    factor2exp2count = defaultdict(lambda: defaultdict(lambda: 0))
    for i, expr in enumerate(expressions):
        factors = dict(expr.factor_list())
        for factor, n in factors.items():
            if not filter(factor): continue
            if factor in factor2exp2count:
                factor2exp2count[factor][n] += 1
            else:
                if i > 0: factor2exp2count[factor][0] = i
                factor2exp2count[factor][n] = 1
        for factor, exps in factor2exp2count.items():
            if factor not in factors:
                exps[0] += 1
    result = SR(1)
    for factor, exp2count in factor2exp2count.items():
        exps = exp2count.keys()
        minn = min(exps)
        maxn = max(exps)
        if minn > 0: result *= factor**minn
        if maxn < 0: result *= factor**maxn
        if minn <= 0 and maxn >= 0:
            bestn = max(exps, key=lambda exp: exp2count[exp])
            result *= factor**bestn
    return result

def simplify_by_factorization(M, x):
    """Try to simplify matrix M by a constant transformation
    that extracts common factors found in M (if any). Return
    the simplified matrix and the transformation.
    """
    n = M.nrows()
    T = identity_matrix(M.base_ring(), n)
    factors = []
    for i in range(n):
        factor = common_factor(
            [M[i,k] for k in range(n) if i != k and not M[i,k].is_zero()] +
            [1/M[k,i] for k in range(n) if k != i and not M[k,i].is_zero()],
            lambda e: x not in e.variables())
        if factor != 1:
            dT = identity_matrix(M.base_ring(), n)
            dT[i,i] = T[i,i] = factor
            M = transform(M, x, dT)
    logger.debug(
            "stripping common factors with this transform:\ndiagonal_matrix([\n  %s\n])" % ",\n  ".join(str(e) for e in T.diagonal()))
    return fuchsia_simplify(M, x), T

#==============================================================================
# Import/Export routines
#==============================================================================

_parser = Parser(make_int=SR, make_float=SR,
        make_var=lambda v: I if v in "Ii" else SR.var(v),
        make_function={"sqrt": sqrt})

def parse(s):
    """Convert a given string to a Sage expression.

    >>> parse("(1 + I*eps)*(1 - I*eps)").simplify_rational()
    eps^2 + 1
    """
    return _parser.parse(s)

def import_matrix_from_file(filename):
    """Read and return a matrix from a named file. Both Mathematica
    and MatrixMarket array formats are supported. The exact format
    will be autodetected.
    """
    with open(filename, 'r') as f:
        hd = f.read(2)
        f.seek(0)
        if hd == "%%":
            return import_matrix_matrixmarket(f)
        elif hd == "{{":
            return import_matrix_mathematica(f)
        else:
            raise ValueError("File '%s' is not in a known matrix format" % filename)

def import_matrix_mathematica(f):
    """Read and return a matrix, stored in the Mathematica format,
    from a file-like object.
    """
    data =  f.read()
    data = data.replace(' ', '').replace('\n', '')
    nr = data.count('},{')+1
    data = data.replace('{', '').replace('}', '')
    data = data.replace(',', '\n')
    data = data.replace('Sqrt[', 'sqrt(').replace(']', ')')
    nc = int((data.count('\n')+1)/nr)

    import io
    sio = io.StringIO("%d %d\n" % (nc, nr) + data)
    return import_matrix_matrixmarket(sio).T

def import_matrix_matrixmarket(f):
    """Read and return a matrix, stored in the MatrixMarket
    array format, from a file-like object.
    """
    lineno = 0
    while True:
        s = f.readline()
        lineno += 1
        if not s.startswith('%'):
            break
    try:
        nrows, ncols = map(int, s.split())
    except ValueError:
        msg = ("Can not determine size of the matrix\n"
               "Make sure the file is in the MatrixMarket array format")
        raise FuchsiaError(msg)
    data = []
    for s in f.readlines():
        lineno += 1
        if s.startswith('%'):
            continue
        try:
            expr = parse(s)
        except SyntaxError:
            msg = ("Malformad expression at line %d:\n%s\n"
                   "Make sure the file is in the MatrixMarket array format")
            raise FuchsiaError(msg % (lineno, s))
        data.append(expr)
    m = matrix(ncols, nrows, data).T
    return m

def export_matrix_to_file(filename, m, fmt="m"):
    """Save matrix into a file in a given matrix format. Set
    'fmt' to "mtx" for MatrixMarket array format, or to "m"
    for Mathematica format.
    """
    with open(filename, 'w') as f:
        if fmt == "mtx":
            export_matrix_matrixmarket(f, m)
        elif fmt == "m":
            export_matrix_mathematica(f, m)
        else:
            raise FuchsiaError("Unknown matrix format '%s'" % fmt)

def export_matrix_mathematica(f, m):
    """Write matrix to a file-like object using Mathematica format."""
    i = 0
    f.write("{")
    nc, nr = m.ncols(), m.nrows()
    for row in m.rows():
        i += 1
        f.write("{")
        j = 0
        for mij in row:
            j += 1
            f.write(str(mij).replace(' ', ''))
            if j < nc:
                f.write(',')
        f.write("}")
        if i < nr:
            f.write(",")
    f.write('}')

def export_matrix_matrixmarket(f, m):
    """Write matrix to a file-like object using MatrixMarket
    array format.
    """
    f.write("%%MatrixMarket matrix array Fuchsia[symbolic] general\n")
    f.write("%d %d\n" % (m.nrows(), m.ncols()))
    i = 0
    for col in m.columns():
        i += 1
        f.write('%% #%d\n' % i)
        for mij in col:
            f.write(str(mij).replace(' ', ''))
            f.write('\n')

#==================================================================================================

def main():
    import getopt
    from   contextlib import contextmanager

    @contextmanager
    def profile(file_stats):
        if file_stats is None:
            yield
            return

        import cProfile, pstats
        profile = cProfile.Profile()
        profile.enable()
        try:
            yield
        finally:
            profile.disable()
            with open(file_stats, 'w') as f:
                stats = pstats.Stats(profile, stream=f).sort_stats("cumulative")
                stats.print_stats(50)

    def usage():
        print(__doc__)
        exit(0)

    try:
        print('\033[35;1mFuchsia v%s\033[0m' % (__version__))
        print("Authors: %s\n" % __author__)
        mpath = tpath = profpath = fmt = None
        M = T = None
        x, y, epsilon = SR.var("x y eps")
        fmt = "m"
        seed = 0
        logger.setLevel(logging.INFO)
        kwargs, args = getopt.gnu_getopt(sys.argv[1:],
                "hvl:f:P:x:y:e:m:t:s:", ["help", "use-maple"])
        for key, value in kwargs:
            if key in ["-h", "--help"]: usage()
            if key == "-f":
                if value not in ["m", "mtx"]:
                    raise getopt.GetoptError(
                            "'%s' is not a known matrix format; "
                            "'m' and 'mtx' are supported" % value)
                fmt = value
            if key == "-l":
                fh = logging.FileHandler(value, "w")
                logger_format = '%(levelname)-8s [%(asctime)s] %(message)s'
                fh.setFormatter(ElapsedTimeFormatter(logger_format))
                logger.addHandler(fh)
            if key == "-v": logger.setLevel(logging.DEBUG)
            if key == "-P": profpath = value
            if key == "-x": x = parse(value)
            if key == "-y": y = parse(value)
            if key == "-e": epsilon = SR.var(value)
            if key == "-m": mpath = value
            if key == "-t": tpath = value
            if key == "-s": seed = int(value)
            if key == "--use-maple":
                global USE_MAPLE
                USE_MAPLE = True
        with profile(profpath):
            if len(args) == 2 and args[0] == 'fuchsify':
                M = import_matrix_from_file(args[1])
                M, t1 = simplify_by_factorization(M, x)
                M, t2 = fuchsify(M, x, seed=seed)
                T = t1*t2
            elif len(args) == 2 and args[0] == 'normalize':
                M = import_matrix_from_file(args[1])
                M, t1 = simplify_by_factorization(M, x)
                M, t2 = normalize(M, x, epsilon, seed=seed)
                T = t1*t2
            elif len(args) == 2 and args[0] == 'factorize':
                M = import_matrix_from_file(args[1])
                M, T = factorize(M, x, epsilon, seed=seed)
            elif len(args) == 2 and args[0] == 'sort':
                M = import_matrix_from_file(args[1])
                M, T, B = block_triangular_form(M)
            elif len(args) == 2 and args[0] == 'reduce':
                m = import_matrix_from_file(args[1])
                M, T = epsilon_form(m, x, epsilon, seed=seed)
            elif len(args) == 2 and args[0] == 'cat':
                M = import_matrix_from_file(args[1])
                T = None
            elif len(args) == 3 and args[0] == 'transform':
                M = import_matrix_from_file(args[1])
                t = import_matrix_from_file(args[2])
                M = transform(M, x, t)
            elif len(args) == 3 and args[0] == 'changevar':
                M = import_matrix_from_file(args[1])
                fy = parse(args[2])
                extravars = set(fy.variables()) - set(M.variables())
                if y not in extravars:
                    raise getopt.GetoptError(
                            "'%s' is missing in '%s'" % (y, fy))
                M = change_variable(M, x, y, fy)
                x = y
            elif len(args) == 2 and args[0] == 'info':
                M = import_matrix_from_file(args[1])
                print("Matrix size: %s" % M.nrows())
                print("Matrix complexity: %s" % matrix_complexity(M))
                print("Matrix expansion:")
                for point, prank in singularities(M, x).items():
                    c0 = matrix_c0(M, x, point, prank)
                    print("  C0[%s=%s; %s] = {rank=%s, complexity=%s}" % (
                        x, point, prank, c0.rank(), matrix_complexity(c0)
                    ))
                    if logger.isEnabledFor(logging.DEBUG):
                        for eigenval in sorted(set(c0.eigenvalues()), key=str):
                            print("    eigenvalue: %s" % (eigenval))
                M = None
            elif len(args) == 0:
                usage()
            else:
                raise getopt.GetoptError("unknown command: %s" % (args))
        if M is not None:
            M = partial_fraction(M, x)
            if mpath is not None:
                export_matrix_to_file(mpath, M, fmt=fmt)
            else:
                print(M)
        if tpath is not None and T is not None:
            T = partial_fraction(T, x)
            export_matrix_to_file(tpath, T, fmt=fmt)
    except getopt.GetoptError as error:
        logger.error("%s", error)
        exit(1)
    except FuchsiaError as e:
        msg, tab = str(e), '  '
        msg = tab+msg.replace('\n', '\n'+tab)
        logger.error(msg)
        exit(1)

if __name__ == '__main__':
    main()

#!/usr/bin/env sage \ Usage:     fuchsia [options] <command> <args>...  Commands:     reduce [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>         find an epsilon form of the given matrix      fuchsify [-x <name>] [-m <path>] [-t <path>] <matrix>         find a transformation that will transform a given matrix         into Fuchsian form      normalize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>         find a transformation that will transform a given Fuchsian         matrix into normalized form      factorize [-x <name>] [-e <name>] [-m <path>] [-t <path>] <matrix>         find a transformation that will make a given normalized         matrix proportional to the infinitesimal parameter      info [-x <name>] <matrix>         show a description of a given matrix, listing singular         points and residue eigenvalues (in verbose more)      sort [-m <path>] [-t <path>] <matrix>         find a block-triangular form of the given matrix      transform [-x <name>] [-m <path>] <matrix> <transform>         transform a given matrix using a given transformation      changevar [-x <name>] [-y <name>] [-m <path>] <matrix> <expr>         transform a given matrix by susbtituting free variable         by a given expression  Options:     -h          show this help message     -f <fmt>    matrix file format: mtx or m (default: m)     -l <path>   write log to this file     -v          produce a more verbose log     -P <path>   save profile report into this file     -x <name>   use this name for the free variable (default: x)     -y <name>   use this name for the new free variable (default: y)     -e <name>   use this name for the infinitesimal parameter (default: eps)     -m <path>   save the resulting matrix into this file     -t <path>   save the resulting transformation into this file     -s <number> use this random seed when applicable (default: 0)     --use-maple speed up calculations by using Maple when possible  Arguments:     <matrix>    read the input matrix from this file     <transform> read the transformation matrix from this file     <expr>      arbitrary expression #@logcall Given a system of differential equations dF/dx = M*F,     and a transformation of base functions F = T*F', compute     and return M', such that dF'/dx = M'*F'.      Note: M' = inverse(T)*(M*T - dT/dx) Same thing as transform(M, x, balance(P, x1, x2, x)), but faster. #assert (P*P - P).is_zero() #@logcall #logger.debug("[%s -> %s] %s" % (x, x0, expr)) #logger.debug("res = %s" % res) Return a limit of expr when x->lim.      The standard 'limit()' function of SageMath does not allow     you to specify the variable, only it's name as a keyword     argument. If you have a variable, and not it's name, use     this function instead. Find values of x around which rational matrix M has     a singularity; return a dictionary with {val: p} entries,     where p is the Poincare rank of M at x=val.      Example:     >>> x = var("x")     >>> s = singularities(matrix([[1/x, 0], [1/(x+1)**3, 1/(x+1)]]), x)     >>> sorted(s.items())     [(-1, 2), (0, 0), (+Infinity, 0)]      >>> s = singularities(matrix([[x, 1, 1/x]]), x)     >>> sorted(s.items())     [(0, 0), (+Infinity, 2)] Return the 0-th coefficient of Taylor expansion of     a matrix M around a finite point x=point, assuming that     M(x->point)~1/(x-point)**exp.      Example:     >>> x = var('x')     >>> matrix_taylor0(matrix([[x/(x-1), 1/x, x, 1]]), x, 1, 1)     [1 0 0 0] Return the 1-th coefficient of Taylor expansion of     a matrix M around a finite point x=point, assuming that     M(x->point)~1/(x-point)**exp.      Example:     >>> x = var('x')     >>> matrix_taylor1(matrix([[x/(x-1), 1/x, x, 1]]), x, 0, 1)     [0 0 0 1] Return the 0-th coefficient of M's expansion at x=point,     assuming Poincare rank of M at that point is p. If point is     +Infinity, return minus the coefficient at the highest power of x.      Examples:     >>> x = var("x")     >>> m = matrix([[1/x, 1/x**2], [1, 1/(x-1)]])     >>> matrix_c0(m, x, 0, 1)     [0 1]     [0 0]     >>> matrix_c0(m, x, 1, 0)     [0 0]     [0 1]     >>> matrix_c0(m, x, oo, 1)     [ 0  0]     [-1  0]     >>> matrix_c0(m*x, x, oo, 2)     [ 0  0]     [-1  0] Return the 1-st coefficient of M's expansion at x=point,     assuming Poincare rank of M at that point is p. If point is     +Infinity, return minus the coefficient at the second-to-highest     power of x.      Examples:     >>> x = var("x")     >>> m = matrix([[1/x, 1/x**2], [1, 1/(x-1)]])     >>> matrix_c1(m, x, 0, 1)     [1 0]     [0 0]     >>> matrix_c1(m, x, oo, 1)     [-1  0]     [ 0 -1] Return matrix residue of M at x=x0, assuming that M's     Poincare rank at x=x0 is zero.      Example:     >>> x = var("x")     >>> m = matrix([[1/x, 2/x], [3/x, 4/x]])     >>> matrix_residue(m, x, 0)     [1 2]     [3 4]     >>> matrix_residue(m, x, oo)     [-1 -2]     [-3 -4] Return True if M is always nilpotent, False otherwise.      Examples:     >>> a, b, c = var('a b c')     >>> matrix_is_nilpotent(matrix([[0,a,b],[0,0,c],[0,0,0]]))     True     >>> matrix_is_nilpotent(matrix([[a,0,0],[0,b,0],[0,0,c]]))     False Return a tuple of Jordan cell sizes from a matrix J in Jordan     normal form.      Examples:     >>> jordan_cell_sizes(matrix([[1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,0,1]]))     (2, 2)     >>> jordan_cell_sizes(matrix([[1,1,0,0],[0,1,1,0],[0,0,1,0],[0,0,0,1]]))     (3, 1)     >>> jordan_cell_sizes(zero_matrix(5,5))     (1, 1, 1, 1, 1) As of SageMath 6.10, 'Matrix.solve_right' method uses a broken     check for solution correctness; this function corrects that check. # Sometimes 'solve_right' returns a giant expression, # which can be simplified into a tiny one; without this # simplification 'is_zero' call below may take forever # to finish. # # Note that the condition above is there to make sure # that an empty matrix is not passed into 'simplify_rational', # otherwise you'll get this error: #   TypeError: unable to make sense of Maxima expression #   'matrix()' in Sage As of SageMath 6.10, 'Matrix.solve_left' method uses a broken     check for solution correctness; this function corrects that check. #================================================================================================== # Transformation routines #================================================================================================== Find a lower block-triangular form of a given matrix.      Return a tuple `(M, T, B)` where:       * `M` is a new matrix;       * `T` is a transformation matrix;       * `B` is a list of tuples (ki, ni), such that M's i-th         diagonal block is given by `M.submatrix(ki, ki, ni, ni)`. #================================================================================================== # Step I: Fuchsify #================================================================================================== Given a system of differential equations of the form dF/dx=M*F,     try to find a transformation T, which will reduce M to Fuchsian     form. Return the transformed M and T. Raise FuchsiaError if     M can not be transformed into Fuchsian form.      Note that such transformations are not unique; you can obtain     different ones by supplying different seeds. Given a system of differential equations of the form dF/dx=M*F,     with M having a singularity around x=v with Poincare rank p,     try to find a transformation T, which will reduce p by 1 (but     possibly introduce another singularity at x=v2). Return the     transformed M and T. Find matrix V, such that it's rows span a left invariant     subspace of A and form a dual basis with columns of U (that     is, V*U=I). #if not matrix_is_nilpotent(A0): #    raise FuchsiaError("matrix is irreducible (non-nilpotent residue)") # We will rely on undocumented behavior of jordan_form() call: # we will take it upon faith that it sorts Jordan cells by # their size, and puts the largest ones closer to (0, 0). #zero_rows = [i for i in range(A0.nrows()) if A0J[i,:].is_zero()] #zero_cols = [j for j in range(A0.nrows()) if A0J[:,j].is_zero()] #assert len(zero_rows) == len(zero_cols) == ncells #_L0 = (invU*A1*U)[zero_rows, zero_cols] #assert (L0 - _L0).is_zero() # No solution found; vectors are independent. # Check Alg.1 promise #================================================================================================== # Step II: Normalize #================================================================================================== Return True if (a Fuchsian) matrix M is normalized, that     is all the eigenvalues of it's residues in x lie in [-1/2, 1/2)     range (in limit eps->0). Return False otherwise.      Examples:     >>> x, e = var("x epsilon")     >>> is_normalized(matrix([[(1+e)/3/x, 0], [0, e/x]]), x, e)     True Return True if diagonal blocks of `m` are normalized, that is all the eigenvalues of theirs     residues in `x` lie in the range [-1/2, 1/2) in eps->0 limit; return False otherwise. Diagonal     blocks are defined by the list `b` which corresponds to the equivalent value returned by the     `block_triangular_form` routine.      Examples:     >>> x, e = var("x epsilon")     >>> are_diagonal_blocks_reduced(matrix([[(1+e)/3/x, 0], [0, e/x]]), [(0,1),(1,1)], x, e)     True Given a lower block-triangular system of differential equations of the form dF/dx=m*F,     find a transformation that will shift all eigenvalues of all residues of all its diagonal     blocks into the range [-1/2, 1/2) in eps->0 limit. Return the transformed matrix m and the     transformation; raise FuchsiaError if such transformation is not found. Diagonal blocks      are defined by the list `b` which corresponds to the equivalent value returned by the     `block_triangular_form` routine. #%d (%d,%d)" % (i,ki,ni)) # M += C*(x-pi)^ki/(x-x0) # M += C*(x-pi)^ki*(x-x0) = C*(x-pi)^ki*{(x-pi) + (pi-x0)} #            print("cache miss: %s" % key) #        else: #            print("cache hit: %s" % key) # singular points # residues in singular points # TODO: check if system's residues have half-integer eigenvalues. # If yes, then inform user about that and let us remove resonances only. # Though, such a form does not allow to construct an eps-form, it # however can be used to find solutions as a power series using # semi-analytical methods. #s = "zeroes:\n  %s\n" % (self.xs,) #@logcall # coP Ci coP (x-pi)^-1 + P Ci P (x-pi)^-1 # coP Ci P -(x-x2) (x-pi)^-1 # P Ci coP -1/(x-x2) (x-pi)^-1 # -P/(x-x2) # coP Ci coP (x-pi)^-1 + P Ci P (x-pi)^-1 # P Ci coP -(x-x1) (x-pi)^-1 # coP Ci P -1/(x-x1) (x-pi)^-1 # P/(x-x1) # coP Ci P (x1-x2)/(x-x1) (x-pi)^-1 # P Ci coP (x2-x1)/(x-x2) (x-pi)^-1 # P/(x-x1) - P/(x-x2) #T = fuchsia_simplify(self.T * balance(P, x1, x2, self.x)) #return RationalSystem(dict(Cmap), self.x, T, self.trace + [("balance", P, x1, x2)]) #@logcall #@logcall #assert(not (x1-x2).is_zero()) #@logcall #    @memoize #    @logcall #    def eigenvectors_left(self): #        return self._eigenvectors(left=True) # #    @memoize #    @logcall #    def eigenvectors_right(self): #        return self._eigenvectors(left=False) # #    def _eigenvectors(self, left=True): #        res = [] #        for xi, ri in self.rs.items(): #            eigenvectors = ri.eigenvectors_left #            if left is not True: #                eigenvectors = ri.eigenvectors_right #            res.append((xi, eigenvectors())) #        return res #    def _eigenvectors(self, left=True): #        res = [] #        eigenvectors = self.m.eigenvectors_left #        if left is not True: #            eigenvectors = self.m.eigenvectors_right #        res = [vvm[0:2] for vvm in eigenvectors()] #        return res # helpers #if any((ex-t).is_zero() for t in printed): # current matrix to balance # matrices balanced so far # Entry point #        cmd = self.cmd_main #        while cmd != None: #            cmd = cmd() #        return # Helpers # Commands Available commands:             b <i> <j>    balance xi, xj points             e            show eigenvectors             h            print this help             q            quit             s <file>     save matrix to file Given a Fuchsian system of differential equations of the     form dF/dx=m*F, find a transformation that will shift all     the eigenvalues of m's residues into [-1/2, 1/2) range (in     the limit eps->0). Return the transformed matrix m and the     transformation. Raise FuchsiaError if such transformation     is not found. # random # points # ev_cum # x0 #================================================================================================== # Step III: Factorize #================================================================================================== Given a normalized Fuchsian system of differential equations:         dF/dx = M(x,epsilon)*F,     try to find a transformation that will factor out an epsilon     from M. Return a transformed M (proportional to epsilon)     and T. Raise FuchsiaError if epsilon can not be factored. # Right now S likely has a number of free variables in # it; we can set them to arbibtrary values, as long as # it'll make S invertible. # We're leaking a bunch of temprary variables here, # which accumulate in SR.variables, but who cares? # We've tried a bunch of substituions, and they didn't # work. Is the matrix at all invertible? Let's check. #================================================================================================== # Helpers #================================================================================================== Try to simplify matrix M by constant transformations that     transform M's residues into their Jordan forms. Return the     simplified matrix and the transformation. If none of the     attempted transformations reduce M's complexity (as measured     by 'matrix_complexity()'), return the original matrix and     the identity transformation. Factorize given expressions, select those factors for     which 'filter(factor)' is True, and return the product of     factors common to all the expressions.      Examples:     >>> x = var("x")     >>> common_factor([x*x-1, x+1], lambda f: True)     x + 1     >>> common_factor([1/x**2, 2/x**3, 3/x**4], lambda f: True)     x^(-2)      Note that if there is a mix of positive and negative exponents     of a given factor, this function will use (one of) the most     frequently occurring exponent:     >>> common_factor([x, 1/x, 2/x**2, 3/x], lambda f: True)     1/x Try to simplify matrix M by a constant transformation     that extracts common factors found in M (if any). Return     the simplified matrix and the transformation. #============================================================================== # Import/Export routines #============================================================================== Convert a given string to a Sage expression.      >>> parse("(1 + I*eps)*(1 - I*eps)").simplify_rational()     eps^2 + 1 Read and return a matrix from a named file. Both Mathematica     and MatrixMarket array formats are supported. The exact format     will be autodetected. Read and return a matrix, stored in the Mathematica format,     from a file-like object. Read and return a matrix, stored in the MatrixMarket     array format, from a file-like object. Save matrix into a file in a given matrix format. Set     'fmt' to "mtx" for MatrixMarket array format, or to "m"     for Mathematica format. Write matrix to a file-like object using Mathematica format. Write matrix to a file-like object using MatrixMarket     array format. #%d\n' % i) #==================================================================================================