'Return `(P, J)` where `J` is a Jordan block
matrix and `P` is a matrix such that
`self == P*J*P**-1`
Parameters
calc_transform : bool
If ``False``, then only `J` is returned.
chop : bool
All matrices are convered to exact types when computing
eigenvalues and eigenvectors.  As a result, there may be
approximation errors.  If ``chop==True``, these errors
will be truncated.
Examples
>>> from sympy import Matrix
>>> m = Matrix([[ 6,  5, -2, -3], [-3, -1,  3,  3], [ 2,  1, -2, -3], [-1,  1,  5,  5]])
>>> P, J = m.jordan_form()
>>> J
Matrix([
[2, 1, 0, 0],
[0, 2, 0, 0],
[0, 0, 2, 1],
[0, 0, 0, 2]])
See Also
jordan_block'

def jordan_form(self, calc_transform=True, **kwargs):

    if (not self.is_square):
        raise NonSquareMatrixError('Only    square    matrices    have    Jordan    forms')
    chop = kwargs.pop('chop', False)
    mat = self
    has_floats = any((v.has(Float) for v in self))
    def restore_floats(*args):
        'If    `has_floats`    is    `True`,    cast    all    `args`    as\n                                                matrices    of    floats.'
        if has_floats:
            args = [m.evalf(chop=chop) for m in args]
        if (len(args) == 1):
            return args[0]
        return args
    mat_cache = {}
    def eig_mat(val, pow):
        'Cache    computations    of    (self    -    val*I)**pow    for    quick\n                                                retrieval'
        if ((val, pow) in mat_cache):
            return mat_cache[(val, pow)]
        if ((val, (pow - 1)) in mat_cache):
            mat_cache[(val, pow)] = (mat_cache[(val, (pow - 1))] * mat_cache[(val, 1)])
        else:
            mat_cache[(val, pow)] = ((mat - (val * self.eye(self.rows))) ** pow)
        return mat_cache[(val, pow)]
    def nullity_chain(val):
        'Calculate    the    sequence        [0,    nullity(E),    nullity(E**2),    ...]\n                                                until    it    is    constant    where    `E    =    self    -    val*I`'
        cols = self.cols
        ret = [0]
        nullity = (cols - eig_mat(val, 1).rank())
        i = 2
        while (nullity != ret[(-1)]):
            ret.append(nullity)
            nullity = (cols - eig_mat(val, i).rank())
            i += 1
        return ret
    def blocks_from_nullity_chain(d):
        'Return    a    list    of    the    size    of    each    Jordan    block.\n                                                If    d_n    is    the    nullity    of    E**n,    then    the    number\n                                                of    Jordan    blocks    of    size    n    is\n\n                                                                2*d_n    -    d_(n-1)    -    d_(n+1)'
        mid = [(((2 * d[n]) - d[(n - 1)]) - d[(n + 1)]) for n in range(1, (len(d) - 1))]
        end = ([(d[(-1)] - d[(-2)])] if (len(d) > 1) else [d[0]])
        return (mid + end)
    def pick_vec(small_basis, big_basis):
        "Picks    a    vector    from    big_basis    that    isn't    in\n                                                the    subspace    spanned    by    small_basis"
        if (len(small_basis) == 0):
            return big_basis[0]
        for v in big_basis:
            (_, pivots) = self.hstack(*(small_basis + [v])).echelon_form(with_pivots=True)
            if (pivots[(-1)] == len(small_basis)):
                return v
    if has_floats:
        mat = mat.applyfunc((lambda x: nsimplify(x, rational=True)))
    eigs = mat.eigenvals()
    if (sum((m for m in eigs.values())) != mat.cols):
        raise MatrixError('Could    not    compute    eigenvalues    for    {}'.format(mat))
    if (len(eigs.keys()) == mat.cols):
        blocks = list(sorted(eigs.keys(), key=default_sort_key))
        jordan_mat = mat.diag(*blocks)
        if (not calc_transform):
            return restore_floats(jordan_mat)
        jordan_basis = [eig_mat(eig, 1).nullspace()[0] for eig in blocks]
        basis_mat = mat.hstack(*jordan_basis)
        return restore_floats(basis_mat, jordan_mat)
    block_structure = []
    for eig in sorted(eigs.keys(), key=default_sort_key):
        chain = nullity_chain(eig)
        block_sizes = blocks_from_nullity_chain(chain)
        size_nums = [((i + 1), num) for (i, num) in enumerate(block_sizes)]
        size_nums.reverse()
        block_structure.extend(((eig, size) for (size, num) in size_nums for _ in range(num)))
    blocks = (mat.jordan_block(size=size, eigenvalue=eig) for (eig, size) in block_structure)
    jordan_mat = mat.diag(*blocks)
    if (not calc_transform):
        return restore_floats(jordan_mat)
    jordan_basis = []
    for eig in sorted(eigs.keys(), key=default_sort_key):
        eig_basis = []
        for (block_eig, size) in block_structure:
            if (block_eig != eig):
                continue
            null_big = eig_mat(eig, size).nullspace()
            null_small = eig_mat(eig, (size - 1)).nullspace()
            vec = pick_vec((null_small + eig_basis), null_big)
            new_vecs = [(eig_mat(eig, i) * vec) for i in range(size)]
            eig_basis.extend(new_vecs)
            jordan_basis.extend(reversed(new_vecs))
    basis_mat = mat.hstack(*jordan_basis)
    return restore_floats(basis_mat, jordan_mat)
