def test_mev():
    output("""\
    reim:{$[0>type x;1 0*x;2=count x;x;'`]};
    mc:{((x[0]*y 0)-x[1]*y 1;(x[0]*y 1)+x[1]*y 0)};
    mmc:{((.qml.mm[x 0]y 0)-.qml.mm[x 1]y 1;(.qml.mm[x 0]y 1)+.qml.mm[x 1]y 0)};
    mev_:{[b;x]
        if[2<>count wv:.qml.mev x;'`length];
        if[not all over prec>=abs
            mmc[flip vc;flip(flip')(reim'')flip x]-
            flip(w:reim'[wv 0])mc'vc:(flip')(reim'')(v:wv 1);'`check];
        / Normalize sign; LAPACK already normalized to real
        v*:1-2*0>{x a?max a:abs x}each vc[;0];
        (?'[prec>=abs w[;1];w[;0];w];?'[b;v;0n])};""")

    for A in eigenvalue_subjects:
        if A.rows <= 3:
            V = []
            for w, n, r in A.eigenvects():
                w = sp.simplify(sp.expand_complex(w))
                if len(r) == 1:
                    r = r[0]
                    r = sp.simplify(sp.expand_complex(r))
                    r = r.normalized() / sp.sign(max(r, key=abs))
                    r = sp.simplify(sp.expand_complex(r))
                else:
                    r = None
                V.extend([(w, r)]*n)
            V.sort(key=lambda (x, _): (-abs(x), -sp.im(x)))
        else:
            Am = mp.matrix(A)
            # extra precision for complex pairs to be equal in sort
            with mp.extradps(mp.mp.dps):
                W, R = mp.eig(Am)
            V = []
            for w, r in zip(W, (R.column(i) for i in range(R.cols))):
                w = mp.chop(w)
                with mp.extradps(mp.mp.dps):
                    _, S, _ = mp.svd(Am - w*mp.eye(A.rows))
                if sum(x == 0 for x in mp.chop(S)) == 1:
                    # nullity 1, so normalized eigenvector is unique
                    r /= mp.norm(r) * mp.sign(max(r, key=abs))
                    r = mp.chop(r)
                else:
                    r = None
                V.append((w, r))
            V.sort(key=lambda (x, _): (-abs(x), -x.imag))
        W, R = zip(*V)
        test("mev_[%sb" % "".join("0" if r is None else "1" for r in R), A,
             (W, [r if r is None else list(r) for r in R]), complex_pair=True)

def test_evalc():
    x = Symbol("x", real=True)
    y = Symbol("y", real=True)
    z = Symbol("z")
    assert ((x+I*y)**2).expand(complex=True) == x**2+2*I*x*y - y**2
    assert expand_complex(z**(2*I)) == I*im(z**(2*I)) + re(z**(2*I))

    assert exp(I*x) != cos(x)+I*sin(x)
    assert exp(I*x).expand(complex=True) == cos(x)+I*sin(x)
    assert exp(I*x+y).expand(complex=True) == exp(y)*cos(x)+I*sin(x)*exp(y)

    assert sin(I*x).expand(complex=True) == I * sinh(x)
    assert sin(x+I*y).expand(complex=True) == sin(x)*cosh(y) + \
            I * sinh(y) * cos(x)

    assert cos(I*x).expand(complex=True) == cosh(x)
    assert cos(x+I*y).expand(complex=True) == cos(x)*cosh(y) - \
            I * sinh(y) * sin(x)

    assert tan(I*x).expand(complex=True) == tanh(x) * I
    assert tan(x+I*y).expand(complex=True) == \
            ((sin(x)*cos(x) + I*cosh(y)*sinh(y)) / (cos(x)**2 + sinh(y)**2)).expand()

    assert sinh(I*x).expand(complex=True) == I * sin(x)
    assert sinh(x+I*y).expand(complex=True) == sinh(x)*cos(y) + \
            I * sin(y) * cosh(x)

    assert cosh(I*x).expand(complex=True) == cos(x)
    assert cosh(x+I*y).expand(complex=True) == cosh(x)*cos(y) + \
            I * sin(y) * sinh(x)

    assert tanh(I*x).expand(complex=True) == tan(x) * I
    assert tanh(x+I*y).expand(complex=True) == \
            ((sinh(x)*cosh(x) + I*cos(y)*sin(y)) / (sinh(x)**2 + cos(y)**2)).expand()

def expand_function(func):
    func = expand_complex(func).as_real_imag()
    result = []
    for part in func:
        part = str(part)
        part = part.replace('re(c)', 'c_real_t')
        part = part.replace('im(c)', 'c_imag_t')
        part = part.replace('re(z)', 'z_real_t')
        part = part.replace('im(z)', 'z_imag_t')
        part = part.replace('abs(', 'T.abs_(')
        part = part.replace('exp(', 'T.exp(')
        part = part.replace('inv(', 'T.inv(')
        part = part.replace('sqr(', 'T.sqr(')
        part = part.replace('sqrt(', 'T.sqrt(')
        part = part.replace('cos(', 'T.cos(')
        part = part.replace('sin(', 'T.sin(')
        part = part.replace('cosh(', 'T.cosh(')
        part = part.replace('sinh(', 'T.sinh(')
        part = part.replace('log(', 'T.log(')
        part = part.replace('arg(', 'T.angle(')
        part = part.replace('I', '1j')
        part = part.replace('im(', 'T.imag(')
        part = part.replace('re(', 'T.real(')
        result.append(part)
    return result

def test_poly():
    output("""\
    poly_:{
        r:{$[prec>=abs(x:reim x)1;x 0;x]} each .qml.poly x;
        r iasc {x:reim x;(abs x 1;neg x 1;x 0)} each r};""")

    for A in poly_subjects:
        A = sp.Poly(A, sp.Symbol("x"))
        if A.degree() <= 3 and all(a.is_real for a in A.all_coeffs()):
            R = []
            for r in sp.roots(A, multiple=True):
                r = sp.simplify(sp.expand_complex(r))
                R.append(r)
            R.sort(key=lambda x: (abs(sp.im(x)), -sp.im(x), sp.re(x)))
        else:
            R = mp.polyroots([mp.mpc(*(a.evalf(mp.mp.dps)).as_real_imag())
                              for a in A.all_coeffs()])
            R.sort(key=lambda x: (abs(x.imag), -x.imag, x.real))
        assert len(R) == A.degree()
        test("poly_", A.all_coeffs(), R, complex_pair=True)

 def simp(x):
     return simplify(expand_trig(expand_complex(expand(x))))

 def simp_exp(arg):
     return expand_complex(exp(arg))

 def simp(x): return simplify(expand_trig(expand_complex(expand(x))))
 def sinc(x): return sin(pi*x)/(pi*x)

def test_issue_4124():
    from sympy import oo
    assert expand_complex(I*oo) == oo*I