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_factorial2_rewrite():
    n = Symbol('n', integer=True)
    assert factorial2(n).rewrite(gamma) == \
        2**(n/2)*Piecewise((1, Eq(Mod(n, 2), 0)), (sqrt(2)/sqrt(pi), Eq(Mod(n, 2), 1)))*gamma(n/2 + 1)
    assert factorial2(2*n).rewrite(gamma) == 2**n*gamma(n + 1)
    assert factorial2(2*n + 1).rewrite(gamma) == \
        sqrt(2)*2**(n + 1/2)*gamma(n + 3/2)/sqrt(pi)

def test_f_distribution():
    d1 = Symbol("d1", positive=True)
    d2 = Symbol("d2", positive=True)

    X = FDistribution("x", d1, d2)
    assert density(X)(x) == (d2**(d2/2)*sqrt((x*d1)**d1 *
        (x*d1 + d2)**(-d1 - d2))*gamma(d1/2 + d2/2)/(x*gamma(d1/2)*gamma(d2/2)))

def test_hankel_transform():
    from sympy import sinh, cosh, gamma, sqrt, exp

    r = Symbol("r")
    k = Symbol("k")
    nu = Symbol("nu")
    m = Symbol("m")
    a = symbols("a")

    assert hankel_transform(1/r, r, k, 0) == 1/k
    assert inverse_hankel_transform(1/k, k, r, 0) == 1/r

    assert hankel_transform(
        1/r**m, r, k, 0) == 2**(-m + 1)*k**(m - 2)*gamma(-m/2 + 1)/gamma(m/2)
    assert inverse_hankel_transform(
        2**(-m + 1)*k**(m - 2)*gamma(-m/2 + 1)/gamma(m/2), k, r, 0) == r**(-m)

    assert hankel_transform(1/r**m, r, k, nu) == (
        2*2**(-m)*k**(m - 2)*gamma(-m/2 + nu/2 + 1)/gamma(m/2 + nu/2))
    assert inverse_hankel_transform(2**(-m + 1)*k**(
        m - 2)*gamma(-m/2 + nu/2 + 1)/gamma(m/2 + nu/2), k, r, nu) == r**(-m)

    assert hankel_transform(r**nu*exp(-a*r), r, k, nu) == \
        2**(nu + 1)*a*k**(-nu - 3)*(a**2/k**2 + 1)**(-nu - S(
                                                     3)/2)*gamma(nu + S(3)/2)/sqrt(pi)
    assert inverse_hankel_transform(
        2**(nu + 1)*a*k**(-nu - 3)*(a**2/k**2 + 1)**(-nu - S(3)/2)*gamma(
        nu + S(3)/2)/sqrt(pi), k, r, nu) == r**nu*exp(-a*r)

def test_uppergamma():
    from sympy import meijerg, exp_polar, I, expint
    assert uppergamma(4, 0) == 6
    assert uppergamma(x, y).diff(y) == -y**(x-1)*exp(-y)
    assert td(uppergamma(randcplx(), y), y)
    assert uppergamma(x, y).diff(x) == \
           uppergamma(x, y)*log(y) + meijerg([], [1, 1], [0, 0, x], [], y)
    assert td(uppergamma(x, randcplx()), x)

    assert uppergamma(S.Half, x) == sqrt(pi)*(1 - erf(sqrt(x)))
    assert not uppergamma(S.Half - 3, x).has(uppergamma)
    assert not uppergamma(S.Half + 3, x).has(uppergamma)
    assert uppergamma(S.Half, x, evaluate=False).has(uppergamma)
    assert tn(uppergamma(S.Half + 3, x, evaluate=False),
              uppergamma(S.Half + 3, x), x)
    assert tn(uppergamma(S.Half - 3, x, evaluate=False),
              uppergamma(S.Half - 3, x), x)

    assert uppergamma(x, y).rewrite(lowergamma) == gamma(x) - lowergamma(x, y)

    assert tn_branch(-3, uppergamma)
    assert tn_branch(-4, uppergamma)
    assert tn_branch(S(1)/3, uppergamma)
    assert tn_branch(pi, uppergamma)
    assert uppergamma(3, exp_polar(4*pi*I)*x) == uppergamma(3, x)
    assert uppergamma(y, exp_polar(5*pi*I)*x) == \
           exp(4*I*pi*y)*uppergamma(y, x*exp_polar(pi*I)) + gamma(y)*(1-exp(4*pi*I*y))
    assert uppergamma(-2, exp_polar(5*pi*I)*x) == \
           uppergamma(-2, x*exp_polar(I*pi)) - 2*pi*I

    assert uppergamma(-2, x) == expint(3, x)/x**2
    assert uppergamma(x, y).rewrite(expint) == y**x*expint(-x + 1, y)

def test_gamma_series():
    assert gamma(x + 1).series(x, 0, 3) == \
        1 - EulerGamma*x + x**2*(EulerGamma**2/2 + pi**2/12) + O(x**3)
    assert gamma(x).series(x, -1, 3) == \
        -1/(x + 1) + EulerGamma - 1 + (x + 1)*(-1 - pi**2/12 - EulerGamma**2/2 + \
       EulerGamma) + (x + 1)**2*(-1 - pi**2/12 - EulerGamma**2/2 + EulerGamma**3/6 - \
       polygamma(2, 1)/6 + EulerGamma*pi**2/12 + EulerGamma) + O((x + 1)**3, (x, -1))

def test_cosine_transform():
    from sympy import sinh, cosh, Si, Ci

    t = symbols("t")
    w = symbols("w")
    a = symbols("a")
    f = Function("f")

    # Test unevaluated form
    assert cosine_transform(f(t), t, w) == CosineTransform(f(t), t, w)
    assert inverse_cosine_transform(f(w), w, t) == InverseCosineTransform(f(w), w, t)

    assert cosine_transform(1/sqrt(t), t, w) == 1/sqrt(w)
    assert inverse_cosine_transform(1/sqrt(w), w, t) == 1/sqrt(t)

    assert cosine_transform(1/(a**2+t**2), t, w) == sqrt(2)*sqrt(pi)*(-sinh(a*w) + cosh(a*w))/(2*a)

    assert cosine_transform(t**(-a), t, w) == 2**(-a + S(1)/2)*w**(a - 1)*gamma((-a + 1)/2)/gamma(a/2)
    assert inverse_cosine_transform(2**(-a + S(1)/2)*w**(a - 1)*gamma(-a/2 + S(1)/2)/gamma(a/2), w, t) == t**(-a)

    assert cosine_transform(exp(-a*t), t, w) == sqrt(2)*a/(sqrt(pi)*(a**2 + w**2))
    assert inverse_cosine_transform(sqrt(2)*a/(sqrt(pi)*(a**2 + w**2)), w, t) == -sinh(a*t) + cosh(a*t)

    assert cosine_transform(exp(-a*sqrt(t))*cos(a*sqrt(t)), t, w) == a*(-sinh(a**2/(2*w)) + cosh(a**2/(2*w)))/(2*w**(S(3)/2))

    assert cosine_transform(1/(a+t), t, w) == -sqrt(2)*((2*Si(a*w) - pi)*sin(a*w) + 2*cos(a*w)*Ci(a*w))/(2*sqrt(pi))
    assert inverse_cosine_transform(sqrt(2)*meijerg(((S(1)/2, 0), ()), ((S(1)/2, 0, 0), (S(1)/2,)), a**2*w**2/4)/(2*pi), w, t) == 1/(a + t)

    assert cosine_transform(1/sqrt(a**2+t**2), t, w) == sqrt(2)*meijerg(((S(1)/2,), ()), ((0, 0), (S(1)/2,)), a**2*w**2/4)/(2*sqrt(pi))
    assert inverse_cosine_transform(sqrt(2)*meijerg(((S(1)/2,), ()), ((0, 0), (S(1)/2,)), a**2*w**2/4)/(2*sqrt(pi)), w, t) == 1/(t*sqrt(a**2/t**2 + 1))

def test_simplify_other():
    assert simplify(sin(x)**2 + cos(x)**2) == 1
    assert simplify(gamma(x + 1)/gamma(x)) == x
    assert simplify(sin(x)**2 + cos(x)**2 + factorial(x)/gamma(x)) == 1 + x
    assert simplify(Eq(sin(x)**2 + cos(x)**2, factorial(x)/gamma(x))) == Eq(1, x)
    nc = symbols('nc', commutative=False)
    assert simplify(x + x*nc) == x*(1 + nc)

def test_studentt():
    nu = Symbol("nu", positive=True)
    x = Symbol("x")

    X = StudentT(nu, symbol=x)
    assert Density(X) == (Lambda(_x, (_x**2/nu + 1)**(-nu/2 - S.Half)
                          *gamma(nu/2 + S.Half)/(sqrt(pi)*sqrt(nu)*gamma(nu/2))))

def test_binomial_rewrite():
    n = Symbol("n", integer=True)
    k = Symbol("k", integer=True)

    assert binomial(n, k).rewrite(factorial) == factorial(n) / (factorial(k) * factorial(n - k))
    assert binomial(n, k).rewrite(gamma) == gamma(n + 1) / (gamma(k + 1) * gamma(n - k + 1))
    assert binomial(n, k).rewrite(ff) == ff(n, k) / factorial(k)

def test_simplify_other():
    assert simplify(sin(x) ** 2 + cos(x) ** 2) == 1
    assert simplify(gamma(x + 1) / gamma(x)) == x
    assert simplify(sin(x) ** 2 + cos(x) ** 2 + factorial(x) / gamma(x)) == 1 + x
    assert simplify(Eq(sin(x) ** 2 + cos(x) ** 2, factorial(x) / gamma(x))) == Eq(x, 1)
    nc = symbols("nc", commutative=False)
    assert simplify(x + x * nc) == x * (1 + nc)
    # issue 6123
    # f = exp(-I*(k*sqrt(t) + x/(2*sqrt(t)))**2)
    # ans = integrate(f, (k, -oo, oo), conds='none')
    ans = I * (
        -pi
        * x
        * exp(-3 * I * pi / 4 + I * x ** 2 / (4 * t))
        * erf(x * exp(-3 * I * pi / 4) / (2 * sqrt(t)))
        / (2 * sqrt(t))
        + pi * x * exp(-3 * I * pi / 4 + I * x ** 2 / (4 * t)) / (2 * sqrt(t))
    ) * exp(-I * x ** 2 / (4 * t)) / (sqrt(pi) * x) - I * sqrt(pi) * (
        -erf(x * exp(I * pi / 4) / (2 * sqrt(t))) + 1
    ) * exp(
        I * pi / 4
    ) / (
        2 * sqrt(t)
    )
    assert simplify(ans) == -(-1) ** (S(3) / 4) * sqrt(pi) / sqrt(t)
    # issue 6370
    assert simplify(2 ** (2 + x) / 4) == 2 ** x

def test_rewrite():
    x, y, z = symbols('x y z')
    f1 = sin(x) + cos(x)
    assert f1.rewrite(cos,exp) == exp(I*x)/2 + sin(x) + exp(-I*x)/2
    assert f1.rewrite([cos],sin) == sin(x) + sin(x + pi/2, evaluate=False)
    f2 = sin(x) + cos(y)/gamma(z)
    assert f2.rewrite(sin,exp) == -I*(exp(I*x) - exp(-I*x))/2 + cos(y)/gamma(z)

def test_jacobi():
    n = Symbol("n")
    a = Symbol("a")
    b = Symbol("b")

    assert jacobi(0, a, b, x) == 1
    assert jacobi(1, a, b, x) == a/2 - b/2 + x*(a/2 + b/2 + 1)

    assert jacobi(n, a, a, x) == RisingFactorial(a + 1, n)*gegenbauer(n, a + S(1)/2, x)/RisingFactorial(2*a + 1, n)
    assert jacobi(n, a, -a, x) == ((-1)**a*(-x + 1)**(-a/2)*(x + 1)**(a/2)*assoc_legendre(n, a, x)*
                                   factorial(-a + n)*gamma(a + n + 1)/(factorial(a + n)*gamma(n + 1)))
    assert jacobi(n, -b, b, x) == ((-x + 1)**(b/2)*(x + 1)**(-b/2)*assoc_legendre(n, b, x)*
                                   gamma(-b + n + 1)/gamma(n + 1))
    assert jacobi(n, 0, 0, x) == legendre(n, x)
    assert jacobi(n, S.Half, S.Half, x) == RisingFactorial(S(3)/2, n)*chebyshevu(n, x)/factorial(n + 1)
    assert jacobi(n, -S.Half, -S.Half, x) == RisingFactorial(S(1)/2, n)*chebyshevt(n, x)/factorial(n)

    X = jacobi(n, a, b, x)
    assert isinstance(X, jacobi)

    assert jacobi(n, a, b, -x) == (-1)**n*jacobi(n, b, a, x)
    assert jacobi(n, a, b, 0) == 2**(-n)*gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(factorial(n)*gamma(a + 1))
    assert jacobi(n, a, b, 1) == RisingFactorial(a + 1, n)/factorial(n)

    m = Symbol("m", positive=True)
    assert jacobi(m, a, b, oo) == oo*RisingFactorial(a + b + m + 1, m)

    assert conjugate(jacobi(m, a, b, x)) == jacobi(m, conjugate(a), conjugate(b), conjugate(x))

    assert diff(jacobi(n,a,b,x), n) == Derivative(jacobi(n, a, b, x), n)
    assert diff(jacobi(n,a,b,x), x) == (a/2 + b/2 + n/2 + S(1)/2)*jacobi(n - 1, a + 1, b + 1, x)

def test_logcombine_1():
    x, y = symbols("x,y")
    a = Symbol("a")
    z, w = symbols("z,w", positive=True)
    b = Symbol("b", real=True)
    assert logcombine(log(x)+2*log(y)) == log(x) + 2*log(y)
    assert logcombine(log(x)+2*log(y), force=True) == log(x*y**2)
    assert logcombine(a*log(w)+log(z)) == a*log(w) + log(z)
    assert logcombine(b*log(z)+b*log(x)) == log(z**b) + b*log(x)
    assert logcombine(b*log(z)-log(w)) == log(z**b/w)
    assert logcombine(log(x)*log(z)) == log(x)*log(z)
    assert logcombine(log(w)*log(x)) == log(w)*log(x)
    assert logcombine(cos(-2*log(z)+b*log(w))) in [cos(log(w**b/z**2)),
                                                   cos(log(z**2/w**b))]
    assert logcombine(log(log(x)-log(y))-log(z), force=True) == \
        log(log((x/y)**(1/z)))
    assert logcombine((2+I)*log(x), force=True) == I*log(x)+log(x**2)
    assert logcombine((x**2+log(x)-log(y))/(x*y), force=True) == \
        log(x**(1/(x*y))*y**(-1/(x*y)))+x/y
    assert logcombine(log(x)*2*log(y)+log(z), force=True) == \
        log(z*y**log(x**2))
    assert logcombine((x*y+sqrt(x**4+y**4)+log(x)-log(y))/(pi*x**Rational(2, 3)*\
        sqrt(y)**3), force=True) == \
        log(x**(1/(pi*x**Rational(2, 3)*sqrt(y)**3))*y**(-1/(pi*\
        x**Rational(2, 3)*sqrt(y)**3))) + sqrt(x**4 + y**4)/(pi*\
        x**Rational(2, 3)*sqrt(y)**3) + x**Rational(1, 3)/(pi*sqrt(y))
    assert logcombine(Eq(log(x), -2*log(y)), force=True) == \
        Eq(log(x*y**2), Integer(0))
    assert logcombine(Eq(y, x*acos(-log(x/y))), force=True) == \
        Eq(y, x*acos(log(y/x)))
    assert logcombine(gamma(-log(x/y))*acos(-log(x/y)), force=True) == \
        acos(log(y/x))*gamma(log(y/x))
    assert logcombine((2+3*I)*log(x), force=True) == \
        log(x**2)+3*I*log(x)
    assert logcombine(Eq(y, -log(x)), force=True) == Eq(y, log(1/x))

def test_sine_transform():
    from sympy import EulerGamma

    t = symbols("t")
    w = symbols("w")
    a = symbols("a")
    f = Function("f")

    # Test unevaluated form
    assert sine_transform(f(t), t, w) == SineTransform(f(t), t, w)
    assert inverse_sine_transform(f(w), w, t) == InverseSineTransform(f(w), w, t)

    assert sine_transform(1 / sqrt(t), t, w) == 1 / sqrt(w)
    assert inverse_sine_transform(1 / sqrt(w), w, t) == 1 / sqrt(t)

    assert sine_transform((1 / sqrt(t)) ** 3, t, w) == sqrt(w) * gamma(S(1) / 4) / (2 * gamma(S(5) / 4))

    assert sine_transform(t ** (-a), t, w) == 2 ** (-a + S(1) / 2) * w ** (a - 1) * gamma(-a / 2 + 1) / gamma(
        (a + 1) / 2
    )
    assert inverse_sine_transform(
        2 ** (-a + S(1) / 2) * w ** (a - 1) * gamma(-a / 2 + 1) / gamma(a / 2 + S(1) / 2), w, t
    ) == t ** (-a)

    assert sine_transform(exp(-a * t), t, w) == sqrt(2) * w / (sqrt(pi) * (a ** 2 + w ** 2))
    assert inverse_sine_transform(sqrt(2) * w / (sqrt(pi) * (a ** 2 + w ** 2)), w, t) == exp(-a * t)

    assert sine_transform(log(t) / t, t, w) == -sqrt(2) * sqrt(pi) * (log(w ** 2) + 2 * EulerGamma) / 4

    assert sine_transform(t * exp(-a * t ** 2), t, w) == sqrt(2) * w * exp(-w ** 2 / (4 * a)) / (4 * a ** (S(3) / 2))
    assert inverse_sine_transform(sqrt(2) * w * exp(-w ** 2 / (4 * a)) / (4 * a ** (S(3) / 2)), w, t) == t * exp(
        -a * t ** 2
    )

def test_gamma_series():
    assert gamma(x + 1).series(x, 0, 3) == \
        1 - EulerGamma*x + x**2*(EulerGamma**2/2 + pi**2/12) + O(x**3)
    assert gamma(x).series(x, -1, 3) == \
        -1/x + EulerGamma - 1 + x*(-1 - pi**2/12 - EulerGamma**2/2 + EulerGamma) \
        + x**2*(-1 - pi**2/12 - EulerGamma**2/2 + EulerGamma**3/6 -
        polygamma(2, 1)/6 + EulerGamma*pi**2/12 + EulerGamma) + O(x**3)

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_lowergamma():
    from sympy import meijerg, exp_polar, I, expint
    assert lowergamma(x, y).diff(y) == y**(x-1)*exp(-y)
    assert td(lowergamma(randcplx(), y), y)
    assert lowergamma(x, y).diff(x) == \
           gamma(x)*polygamma(0, x) - uppergamma(x, y)*log(y) \
           + meijerg([], [1, 1], [0, 0, x], [], y)

    assert lowergamma(S.Half, x) == sqrt(pi)*erf(sqrt(x))
    assert not lowergamma(S.Half - 3, x).has(lowergamma)
    assert not lowergamma(S.Half + 3, x).has(lowergamma)
    assert lowergamma(S.Half, x, evaluate=False).has(lowergamma)
    assert tn(lowergamma(S.Half + 3, x, evaluate=False),
              lowergamma(S.Half + 3, x), x)
    assert tn(lowergamma(S.Half - 3, x, evaluate=False),
              lowergamma(S.Half - 3, x), x)

    assert lowergamma(x, y).rewrite(uppergamma) == gamma(x) - uppergamma(x, y)

    assert tn_branch(-3, lowergamma)
    assert tn_branch(-4, lowergamma)
    assert tn_branch(S(1)/3, lowergamma)
    assert tn_branch(pi, lowergamma)
    assert lowergamma(3, exp_polar(4*pi*I)*x) == lowergamma(3, x)
    assert lowergamma(y, exp_polar(5*pi*I)*x) == \
           exp(4*I*pi*y)*lowergamma(y, x*exp_polar(pi*I))
    assert lowergamma(-2, exp_polar(5*pi*I)*x) == \
           lowergamma(-2, x*exp_polar(I*pi)) + 2*pi*I

    assert lowergamma(x, y).rewrite(expint) == -y**x*expint(-x + 1, y) + gamma(x)
    k = Symbol('k', integer=True)
    assert lowergamma(k, y).rewrite(expint) == -y**k*expint(-k + 1, y) + gamma(k)
    k = Symbol('k', integer=True, positive=False)
    assert lowergamma(k, y).rewrite(expint) == lowergamma(k, y)

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_betaprime():
    alpha = Symbol("alpha", positive=True)
    beta = Symbol("beta", positive=True)

    X = BetaPrime('x', alpha, beta)
    assert density(X)(x) == (x**(alpha - 1)*(x + 1)**(-alpha - beta)
                          *gamma(alpha + beta)/(gamma(alpha)*gamma(beta)))