def test_E():
    assert E(z, 0) == z
    assert E(0, m) == 0
    assert E(i*pi/2, m) == i*E(m)
    assert E(z, oo) == zoo
    assert E(z, -oo) == zoo
    assert E(0) == pi/2
    assert E(1) == 1
    assert E(oo) == I*oo
    assert E(-oo) == oo
    assert E(zoo) == zoo

    assert E(-z, m) == -E(z, m)

    assert E(z, m).diff(z) == sqrt(1 - m*sin(z)**2)
    assert E(z, m).diff(m) == (E(z, m) - F(z, m))/(2*m)
    assert E(z).diff(z) == (E(z) - K(z))/(2*z)
    r = randcplx()
    assert td(E(r, m), m)
    assert td(E(z, r), z)
    assert td(E(z), z)

    mi = Symbol('m', real=False)
    assert E(z, mi).conjugate() == E(z.conjugate(), mi.conjugate())
    mr = Symbol('m', real=True, negative=True)
    assert E(z, mr).conjugate() == E(z.conjugate(), mr)

    assert E(z).rewrite(hyper) == (pi/2)*hyper((-S.Half, S.Half), (S.One,), z)
    assert tn(E(z), (pi/2)*hyper((-S.Half, S.Half), (S.One,), z))
    assert E(z).rewrite(meijerg) == \
        -meijerg(((S.Half, S(3)/2), []), ((S.Zero,), (S.Zero,)), -z)/4
    assert tn(E(z), -meijerg(((S.Half, S(3)/2), []), ((S.Zero,), (S.Zero,)), -z)/4)

def test_polynomial():
    from sympy import oo
    assert hyperexpand(hyper([], [-1], z)) == oo
    assert hyperexpand(hyper([-2], [-1], z)) == oo
    assert hyperexpand(hyper([0, 0], [-1], z)) == 1
    assert can_do([-5, -2, randcplx(), randcplx()], [-10, randcplx()])
    assert hyperexpand(hyper((-1, 1), (-2,), z)) == 1 + z/2

def test_plan():
    assert devise_plan(Hyper_Function([0], ()),
            Hyper_Function([0], ()), z) == []
    with raises(ValueError):
        devise_plan(Hyper_Function([1], ()), Hyper_Function((), ()), z)
    with raises(ValueError):
        devise_plan(Hyper_Function([2], [1]), Hyper_Function([2], [2]), z)
    with raises(ValueError):
        devise_plan(Hyper_Function([2], []), Hyper_Function([S("1/2")], []), z)

    # We cannot use pi/(10000 + n) because polys is insanely slow.
    a1, a2, b1 = map(lambda n: randcplx(n), range(3))
    b1 += 2*I
    h = hyper([a1, a2], [b1], z)

    h2 = hyper((a1 + 1, a2), [b1], z)
    assert tn(apply_operators(h,
        devise_plan(Hyper_Function((a1 + 1, a2), [b1]),
            Hyper_Function((a1, a2), [b1]), z), op),
        h2, z)

    h2 = hyper((a1 + 1, a2 - 1), [b1], z)
    assert tn(apply_operators(h,
        devise_plan(Hyper_Function((a1 + 1, a2 - 1), [b1]),
            Hyper_Function((a1, a2), [b1]), z), op),
        h2, z)

def test_airybiprime():
    z = Symbol('z', real=False)
    t = Symbol('t', negative=True)
    p = Symbol('p', positive=True)

    assert isinstance(airybiprime(z), airybiprime)

    assert airybiprime(0) == 3**(S(1)/6)/gamma(S(1)/3)
    assert airybiprime(oo) == oo
    assert airybiprime(-oo) == 0

    assert diff(airybiprime(z), z) == z*airybi(z)

    assert series(airybiprime(z), z, 0, 3) == (
        3**(S(1)/6)/gamma(S(1)/3) + 3**(S(5)/6)*z**2/(6*gamma(S(2)/3)) + O(z**3))

    assert airybiprime(z).rewrite(hyper) == (
        3**(S(5)/6)*z**2*hyper((), (S(5)/3,), z**S(3)/9)/(6*gamma(S(2)/3)) +
        3**(S(1)/6)*hyper((), (S(1)/3,), z**S(3)/9)/gamma(S(1)/3))

    assert isinstance(airybiprime(z).rewrite(besselj), airybiprime)
    assert airyai(t).rewrite(besselj) == (
        sqrt(-t)*(besselj(-S(1)/3, 2*(-t)**(S(3)/2)/3) +
                  besselj(S(1)/3, 2*(-t)**(S(3)/2)/3))/3)
    assert airybiprime(z).rewrite(besseli) == (
        sqrt(3)*(z**2*besseli(S(2)/3, 2*z**(S(3)/2)/3)/(z**(S(3)/2))**(S(2)/3) +
                 (z**(S(3)/2))**(S(2)/3)*besseli(-S(2)/3, 2*z**(S(3)/2)/3))/3)
    assert airybiprime(p).rewrite(besseli) == (
        sqrt(3)*p*(besseli(-S(2)/3, 2*p**(S(3)/2)/3) + besseli(S(2)/3, 2*p**(S(3)/2)/3))/3)

    assert expand_func(airybiprime(2*(3*z**5)**(S(1)/3))) == (
        sqrt(3)*(z**(S(5)/3)/(z**5)**(S(1)/3) - 1)*airyaiprime(2*3**(S(1)/3)*z**(S(5)/3))/2 +
        (z**(S(5)/3)/(z**5)**(S(1)/3) + 1)*airybiprime(2*3**(S(1)/3)*z**(S(5)/3))/2)

def test_to_expr():
    x = symbols('x')
    R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx')
    p = HolonomicFunction(Dx - 1, x, 0, [1]).to_expr()
    q = exp(x)
    assert p == q
    p = HolonomicFunction(Dx**2 + 1, x, 0, [1, 0]).to_expr()
    q = cos(x)
    assert p == q
    p = HolonomicFunction(Dx**2 - 1, x, 0, [1, 0]).to_expr()
    q = cosh(x)
    assert p == q
    p = HolonomicFunction(2 + (4*x - 1)*Dx + \
        (x**2 - x)*Dx**2, x, 0, [1, 2]).to_expr().expand()
    q = 1/(x**2 - 2*x + 1)
    assert p == q
    p = expr_to_holonomic(sin(x)**2/x).integrate((x, 0, x)).to_expr()
    q = (sin(x)**2/x).integrate((x, 0, x))
    assert p == q
    C_0, C_1, C_2, C_3 = symbols('C_0, C_1, C_2, C_3')
    p = expr_to_holonomic(log(1+x**2)).to_expr()
    q = C_2*log(x**2 + 1)
    assert p == q
    p = expr_to_holonomic(log(1+x**2)).diff().to_expr()
    q = C_0*x/(x**2 + 1)
    assert p == q
    p = expr_to_holonomic(erf(x) + x).to_expr()
    q = 3*C_3*x - 3*sqrt(pi)*C_3*erf(x)/2 + x + 2*x/sqrt(pi)
    assert p == q
    p = expr_to_holonomic(sqrt(x), x0=1).to_expr()
    assert p == sqrt(x)
    assert expr_to_holonomic(sqrt(x)).to_expr() == sqrt(x)
    p = expr_to_holonomic(sqrt(1 + x**2)).to_expr()
    assert p == sqrt(1+x**2)
    p = expr_to_holonomic((2*x**2 + 1)**(S(2)/3)).to_expr()
    assert p == (2*x**2 + 1)**(S(2)/3)
    p = expr_to_holonomic(sqrt(-x**2+2*x)).to_expr()
    assert p == sqrt(x)*sqrt(-x + 2)
    p = expr_to_holonomic((-2*x**3+7*x)**(S(2)/3)).to_expr()
    q = x**(S(2)/3)*(-2*x**2 + 7)**(S(2)/3)
    assert p == q
    p = from_hyper(hyper((-2, -3), (S(1)/2, ), x))
    s = hyperexpand(hyper((-2, -3), (S(1)/2, ), x))
    D_0 = Symbol('D_0')
    C_0 = Symbol('C_0')
    assert (p.to_expr().subs({C_0:1, D_0:0}) - s).simplify() == 0
    p.y0 = {0: [1], S(1)/2: [0]}
    assert p.to_expr() == s
    assert expr_to_holonomic(x**5).to_expr() == x**5
    assert expr_to_holonomic(2*x**3-3*x**2).to_expr().expand() == \
        2*x**3-3*x**2
    a = symbols("a")
    p = (expr_to_holonomic(1.4*x)*expr_to_holonomic(a*x, x)).to_expr()
    q = 1.4*a*x**2
    assert p == q
    p = (expr_to_holonomic(1.4*x)+expr_to_holonomic(a*x, x)).to_expr()
    q = x*(a + 1.4)
    assert p == q
    p = (expr_to_holonomic(1.4*x)+expr_to_holonomic(x)).to_expr()
    assert p == 2.4*x

def test_K():
    assert K(0) == pi / 2
    assert K(S(1) / 2) == 8 * pi ** (S(3) / 2) / gamma(-S(1) / 4) ** 2
    assert K(1) == zoo
    assert K(-1) == gamma(S(1) / 4) ** 2 / (4 * sqrt(2 * pi))
    assert K(oo) == 0
    assert K(-oo) == 0
    assert K(I * oo) == 0
    assert K(-I * oo) == 0
    assert K(zoo) == 0

    assert K(z).diff(z) == (E(z) - (1 - z) * K(z)) / (2 * z * (1 - z))
    assert td(K(z), z)

    zi = Symbol("z", real=False)
    assert K(zi).conjugate() == K(zi.conjugate())
    zr = Symbol("z", real=True, negative=True)
    assert K(zr).conjugate() == K(zr)

    assert K(z).rewrite(hyper) == (pi / 2) * hyper((S.Half, S.Half), (S.One,), z)
    assert tn(K(z), (pi / 2) * hyper((S.Half, S.Half), (S.One,), z))
    assert K(z).rewrite(meijerg) == meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z) / 2
    assert tn(K(z), meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z) / 2)

    assert K(z).series(
        z
    ) == pi / 2 + pi * z / 8 + 9 * pi * z ** 2 / 128 + 25 * pi * z ** 3 / 512 + 1225 * pi * z ** 4 / 32768 + 3969 * pi * z ** 5 / 131072 + O(
        z ** 6
    )

def test_to_hyper():
    x = symbols('x')
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
    p = HolonomicFunction(Dx - 2, x, 0, [3]).to_hyper()
    q = 3 * hyper([], [], 2*x)
    assert p == q
    p = hyperexpand(HolonomicFunction((1 + x) * Dx - 3, x, 0, [2]).to_hyper()).expand()
    q = 2*x**3 + 6*x**2 + 6*x + 2
    assert p == q
    p = HolonomicFunction((1 + x)*Dx**2 + Dx, x, 0, [0, 1]).to_hyper()
    q = -x**2*hyper((2, 2, 1), (2, 3), -x)/2 + x
    assert p == q
    p = HolonomicFunction(2*x*Dx + Dx**2, x, 0, [0, 2/sqrt(pi)]).to_hyper()
    q = 2*x*hyper((1/2,), (3/2,), -x**2)/sqrt(pi)
    assert p == q
    p = hyperexpand(HolonomicFunction(2*x*Dx + Dx**2, x, 0, [1, -2/sqrt(pi)]).to_hyper())
    q = erfc(x)
    assert p.rewrite(erfc) == q
    p =  hyperexpand(HolonomicFunction((x**2 - 1) + x*Dx + x**2*Dx**2,
        x, 0, [0, S(1)/2]).to_hyper())
    q = besselj(1, x)
    assert p == q
    p = hyperexpand(HolonomicFunction(x*Dx**2 + Dx + x, x, 0, [1, 0]).to_hyper())
    q = besselj(0, x)
    assert p == q

def test_airyai():
    z = Symbol('z', real=False)
    t = Symbol('t', negative=True)
    p = Symbol('p', positive=True)

    assert isinstance(airyai(z), airyai)

    assert airyai(0) == 3**(S(1)/3)/(3*gamma(S(2)/3))
    assert airyai(oo) == 0
    assert airyai(-oo) == 0

    assert diff(airyai(z), z) == airyaiprime(z)

    assert series(airyai(z), z, 0, 3) == (
        3**(S(5)/6)*gamma(S(1)/3)/(6*pi) - 3**(S(1)/6)*z*gamma(S(2)/3)/(2*pi) + O(z**3))

    assert airyai(z).rewrite(hyper) == (
        -3**(S(2)/3)*z*hyper((), (S(4)/3,), z**S(3)/9)/(3*gamma(S(1)/3)) +
         3**(S(1)/3)*hyper((), (S(2)/3,), z**S(3)/9)/(3*gamma(S(2)/3)))

    assert isinstance(airyai(z).rewrite(besselj), airyai)
    assert airyai(t).rewrite(besselj) == (
        sqrt(-t)*(besselj(-S(1)/3, 2*(-t)**(S(3)/2)/3) +
                  besselj(S(1)/3, 2*(-t)**(S(3)/2)/3))/3)
    assert airyai(z).rewrite(besseli) == (
        -z*besseli(S(1)/3, 2*z**(S(3)/2)/3)/(3*(z**(S(3)/2))**(S(1)/3)) +
         (z**(S(3)/2))**(S(1)/3)*besseli(-S(1)/3, 2*z**(S(3)/2)/3)/3)
    assert airyai(p).rewrite(besseli) == (
        sqrt(p)*(besseli(-S(1)/3, 2*p**(S(3)/2)/3) -
                 besseli(S(1)/3, 2*p**(S(3)/2)/3))/3)

    assert expand_func(airyai(2*(3*z**5)**(S(1)/3))) == (
        -sqrt(3)*(-1 + (z**5)**(S(1)/3)/z**(S(5)/3))*airybi(2*3**(S(1)/3)*z**(S(5)/3))/6 +
         (1 + (z**5)**(S(1)/3)/z**(S(5)/3))*airyai(2*3**(S(1)/3)*z**(S(5)/3))/2)

def test_K():
    assert K(0) == pi/2
    assert K(S(1)/2) == 8*pi**(S(3)/2)/gamma(-S(1)/4)**2
    assert K(1) == zoo
    assert K(-1) == gamma(S(1)/4)**2/(4*sqrt(2*pi))
    assert K(oo) == 0
    assert K(-oo) == 0
    assert K(I*oo) == 0
    assert K(-I*oo) == 0
    assert K(zoo) == 0

    assert K(z).diff(z) == (E(z) - (1 - z)*K(z))/(2*z*(1 - z))
    assert td(K(z), z)

    zi = Symbol('z', real=False)
    assert K(zi).conjugate() == K(zi.conjugate())
    zr = Symbol('z', real=True, negative=True)
    assert K(zr).conjugate() == K(zr)

    assert K(z).rewrite(hyper) == \
        (pi/2)*hyper((S.Half, S.Half), (S.One,), z)
    assert tn(K(z), (pi/2)*hyper((S.Half, S.Half), (S.One,), z))
    assert K(z).rewrite(meijerg) == \
        meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z)/2
    assert tn(K(z), meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z)/2)

def test_hyper_unpolarify():
    from sympy import exp_polar
    a = exp_polar(2*pi*I)*x
    b = x
    assert hyper([], [], a).argument == b
    assert hyper([0], [], a).argument == a
    assert hyper([0], [0], a).argument == b
    assert hyper([0, 1], [0], a).argument == a

def test_hyperexpand_parametric():
    assert (
        hyperexpand(hyper([a, S(1) / 2 + a], [S(1) / 2], z))
        == (1 + sqrt(z)) ** (-2 * a) / 2 + (1 - sqrt(z)) ** (-2 * a) / 2
    )
    assert hyperexpand(hyper([a, -S(1) / 2 + a], [2 * a], z)) == 2 ** (2 * a - 1) * ((-z + 1) ** (S(1) / 2) + 1) ** (
        -2 * a + 1
    )

def test_hyperexpand():
    # Luke, Y. L. (1969), The Special Functions and Their Approximations,
    # Volume 1, section 6.2

    assert hyperexpand(hyper([], [], z)) == exp(z)
    assert hyperexpand(hyper([1, 1], [2], -z) * z) == log(1 + z)
    assert hyperexpand(hyper([], [S.Half], -z ** 2 / 4)) == cos(z)
    assert hyperexpand(z * hyper([], [S("3/2")], -z ** 2 / 4)) == sin(z)
    assert hyperexpand(hyper([S("1/2"), S("1/2")], [S("3/2")], z ** 2) * z) == asin(z)

def test_limits():
    k, x = symbols('k, x')
    assert hyper((1,), (S(4)/3, S(5)/3), k**2).series(k) == \
           hyper((1,), (S(4)/3, S(5)/3), 0) + \
           9*k**2*hyper((2,), (S(7)/3, S(8)/3), 0)/20 + \
           81*k**4*hyper((3,), (S(10)/3, S(11)/3), 0)/1120 + \
           O(k**6) # issue 6350
    assert limit(meijerg((), (), (1,), (0,), -x), x, 0) == \
            meijerg(((), ()), ((1,), (0,)), 0) # issue 6052

def test_hyper_rewrite_sum():
    from sympy import RisingFactorial, factorial, Dummy, Sum
    _k = Dummy("k")
    assert replace_dummy(hyper((1, 2), (1, 3), x).rewrite(Sum), _k) == \
        Sum(x**_k / factorial(_k) * RisingFactorial(2, _k) /
            RisingFactorial(3, _k), (_k, 0, oo))

    assert hyper((1, 2, 3), (-1, 3), z).rewrite(Sum) == \
        hyper((1, 2, 3), (-1, 3), z)

def test_hyperexpand_special():
    assert hyperexpand(hyper([a, b], [c], 1)) == \
           gamma(c)*gamma(c - a - b)/gamma(c - a)/gamma(c - b)
    assert hyperexpand(hyper([a, b], [1 + a - b], -1)) == \
           gamma(1 + a/2)*gamma(1 + a - b)/gamma(1 + a)/gamma(1 + a/2 - b)
    assert hyperexpand(hyper([a, b], [1 + b - a], -1)) == \
           gamma(1 + b/2)*gamma(1 + b - a)/gamma(1 + b)/gamma(1 + b/2 - a)
    assert hyperexpand(meijerg([1 - z - a/2], [1 - z + a/2], [b/2], [-b/2], 1)) == \
           gamma(1 - 2*z)*gamma(z + a/2 + b/2)/gamma(1 - z + a/2 - b/2) \
           /gamma(1 - z - a/2 + b/2)/gamma(1 - z + a/2 + b/2)

def test_hyperexpand():
    # Luke, Y. L. (1969), The Special Functions and Their Approximations,
    # Volume 1, section 6.2

    assert hyperexpand(hyper([], [], z)) == exp(z)
    assert hyperexpand(hyper([1, 1], [2], -z)*z) == log(1 + z)
    assert hyperexpand(hyper([], [S.Half], -z**2/4)) == cos(z)
    assert hyperexpand(z*hyper([], [S('3/2')], -z**2/4)) == sin(z)
    assert hyperexpand(hyper([S('1/2'), S('1/2')], [S('3/2')], z**2)*z) \
        == asin(z)
    assert isinstance(Sum(binomial(2, z)*z**2, (z, 0, a)).doit(), Expr)

def can_do(ap, bq, numerical=True):
    r = hyperexpand(hyper(ap, bq, z))
    if r.has(hyper):
        return False

    if not numerical:
        return True

    repl = {}
    for n, a in enumerate(r.free_symbols - set([z])):
        repl[a] = randcplx(n)
    return tn(hyper(ap, bq, z).subs(repl), r.subs(repl), z)

def test_plan_derivatives():
    a1, a2, a3 = 1, 2, S("1/2")
    b1, b2 = 3, S("5/2")
    h = hyper((a1, a2, a3), (b1, b2), z)
    h2 = hyper((a1 + 1, a2 + 1, a3 + 2), (b1 + 1, b2 + 1), z)
    ops = devise_plan(IndexPair((a1 + 1, a2 + 1, a3 + 2), (b1 + 1, b2 + 1)), IndexPair((a1, a2, a3), (b1, b2)), z)
    f = Formula((a1, a2, a3), (b1, b2), z, h, [])
    deriv = make_derivative_operator(f.M, z)
    assert tn((apply_operators(f.C, ops, deriv) * f.B)[0], h2, z)

    h2 = hyper((a1, a2 - 1, a3 - 2), (b1 - 1, b2 - 1), z)
    ops = devise_plan(IndexPair((a1, a2 - 1, a3 - 2), (b1 - 1, b2 - 1)), IndexPair((a1, a2, a3), (b1, b2)), z)
    assert tn((apply_operators(f.C, ops, deriv) * f.B)[0], h2, z)

def test_from_hyper():
    x = symbols('x')
    R, Dx = DifferentialOperators(QQ.old_poly_ring(x), 'Dx')
    p = hyper([1, 1], [S(3)/2], x**2/4)
    q = HolonomicFunction((4*x) + (5*x**2 - 8)*Dx + (x**3 - 4*x)*Dx**2, x, 1, [2*sqrt(3)*pi/9, -4*sqrt(3)*pi/27 + 4/3])
    r = from_hyper(p)
    assert r == q
    p = from_hyper(hyper([1], [S(3)/2], x**2/4))
    q = HolonomicFunction(-x + (-x**2/2 + 2)*Dx + x*Dx**2, x)
    x0 = 1
    y0 = '[sqrt(pi)*exp(1/4)*erf(1/2), -sqrt(pi)*exp(1/4)*erf(1/2)/2 + 1]'
    assert sstr(p.y0) == y0
    assert q.annihilator == p.annihilator

def test_hyper():
    raises(TypeError, 'hyper(1, 2, z)')

    assert hyper((1, 2),(1,), z) == hyper(Tuple(1, 2), Tuple(1), z)

    h = hyper((1, 2), (3, 4, 5), z)
    assert h.ap == Tuple(1, 2)
    assert h.bq == Tuple(3, 4, 5)
    assert h.argument == z
    assert h.is_commutative is True

    # just a few checks to make sure that all arguments go where they should
    assert tn(hyper(Tuple(), Tuple(), z), exp(z), z)
    assert tn(z*hyper((1, 1), Tuple(2), -z), log(1 + z), z)

    # differentiation
    h = hyper((randcplx(), randcplx(), randcplx()), (randcplx(), randcplx()), z)
    assert td(h, z)

    a1, a2, b1, b2, b3 = symbols('a1:3, b1:4')
    assert hyper((a1, a2), (b1, b2, b3), z).diff(z) == \
             a1*a2/(b1*b2*b3) * hyper((a1+1, a2+1), (b1+1, b2+1, b3+1), z)

    # differentiation wrt parameters is not supported
    raises(NotImplementedError, 'hyper((z,), (), z).diff(z)')