Python sympy.hyperexpand函数代码示例

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_mellin_transform2():
    MT = mellin_transform
    # TODO we cannot currently do these (needs summation of 3F2(-1))
    #      this also implies that they cannot be written as a single g-function
    #      (although this is possible)
    mt = MT(log(x)/(x + 1), x, s)
    assert mt[1:] == ((0, 1), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)
    mt = MT(log(x)**2/(x + 1), x, s)
    assert mt[1:] == ((0, 1), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)
    mt = MT(log(x)/(x + 1)**2, x, s)
    assert mt[1:] == ((0, 2), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)

 def _eval_expand_func(self, **hints):
     from sympy import gamma, hyperexpand
     if len(self.ap) == 2 and len(self.bq) == 1 and self.argument == 1:
         a, b = self.ap
         c = self.bq[0]
         return gamma(c)*gamma(c - a - b)/gamma(c - a)/gamma(c - b)
     return hyperexpand(self)

def test_lookup_table():
    from random import uniform, randrange
    from sympy import Add
    from sympy.integrals.meijerint import z as z_dummy

    table = {}
    _create_lookup_table(table)
    for _, l in sorted(table.items()):
        for formula, terms, cond, hint in sorted(l, key=default_sort_key):
            subs = {}
            for a in list(formula.free_symbols) + [z_dummy]:
                if hasattr(a, "properties") and a.properties:
                    # these Wilds match positive integers
                    subs[a] = randrange(1, 10)
                else:
                    subs[a] = uniform(1.5, 2.0)
            if not isinstance(terms, list):
                terms = terms(subs)

            # First test that hyperexpand can do this.
            expanded = [hyperexpand(g) for (_, g) in terms]
            assert all(x.is_Piecewise or not x.has(meijerg) for x in expanded)

            # Now test that the meijer g-function is indeed as advertised.
            expanded = Add(*[f * x for (f, x) in terms])
            a, b = formula.n(subs=subs), expanded.n(subs=subs)
            r = min(abs(a), abs(b))
            if r < 1:
                assert abs(a - b).n() <= 1e-10
            else:
                assert (abs(a - b) / r).n() <= 1e-10

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 _inverse_mellin_transform(F, s, x_, strip, as_meijerg=False):
    """ A helper for the real inverse_mellin_transform function, this one here
        assumes x to be real and positive. """
    from sympy import (expand, expand_mul, hyperexpand, meijerg, And, Or,
                       arg, pi, re, factor, Heaviside, gamma, Add)
    x = _dummy('t', 'inverse-mellin-transform', F, positive=True)
    # Actually, we won't try integration at all. Instead we use the definition
    # of the Meijer G function as a fairly general inverse mellin transform.
    F = F.rewrite(gamma)
    for g in [factor(F), expand_mul(F), expand(F)]:
        if g.is_Add:
            # do all terms separately
            ress = [_inverse_mellin_transform(G, s, x, strip, as_meijerg,
                                              noconds=False) \
                    for G in g.args]
            conds = [p[1] for p in ress]
            ress = [p[0] for p in ress]
            res = Add(*ress)
            if not as_meijerg:
                res = factor(res, gens=res.atoms(Heaviside))
            return res.subs(x, x_), And(*conds)

        try:
            a, b, C, e, fac = _rewrite_gamma(g, s, strip[0], strip[1])
        except IntegralTransformError:
            continue
        G = meijerg(a, b, C/x**e)
        if as_meijerg:
            h = G
        else:
            h = hyperexpand(G)
            if h.is_Piecewise and len(h.args) == 3:
                # XXX we break modularity here!
                h = Heaviside(x - abs(C))*h.args[0].args[0] \
                  + Heaviside(abs(C) - x)*h.args[1].args[0]
        # We must ensure that the intgral along the line we want converges,
        # and return that value.
        # See [L], 5.2
        cond = [abs(arg(G.argument)) < G.delta*pi]
        # Note: we allow ">=" here, this corresponds to convergence if we let
        # limits go to oo symetrically. ">" corresponds to absolute convergence.
        cond += [And(Or(len(G.ap) != len(G.bq), 0 >= re(G.nu) + 1),
                     abs(arg(G.argument)) == G.delta*pi)]
        cond = Or(*cond)
        if cond is False:
            raise IntegralTransformError('Inverse Mellin', F, 'does not converge')
        return (h*fac).subs(x, x_), cond

    raise IntegralTransformError('Inverse Mellin', F, '')

def test_laplace_transform():
    LT = laplace_transform
    a, b, c, = symbols('a b c', positive=True)
    t = symbols('t')

    # basic tests from wikipedia

    assert LT((t-a)**b*exp(-c*(t-a))*Heaviside(t-a), t, s) \
           == ((s + c)**(-b - 1)*exp(-a*s)*gamma(b + 1), -c, True)
    assert LT(t**a, t, s) == (s**(-a - 1)*gamma(a + 1), 0, True)
    assert LT(Heaviside(t), t, s) == (1/s, 0, True)
    assert LT(Heaviside(t - a), t, s) == (exp(-a*s)/s, 0, True)
    assert LT(1 - exp(-a*t), t, s) == (a/(s*(a + s)), 0, True)

    assert LT((exp(2*t)-1)*exp(-b - t)*Heaviside(t)/2, t, s, noconds=True) \
           == exp(-b)/(s**2 - 1)

    assert LT(exp(t), t, s)[0:2] == (1/(s-1), 1)
    assert LT(exp(2*t), t, s)[0:2] == (1/(s-2), 2)
    assert LT(exp(a*t), t, s)[0:2] == (1/(s-a), a)

    lt = LT(log(t/a), t, s)
    assert lt[1:] == (0, True)
    # TODO hyperexpand is not clever enough to recognise this on its own
    assert hyperexpand(lt[0], allow_hyper=True) == (-log(a*s) - EulerGamma)/s

    assert LT(erf(t), t, s) == ((-erf(s/2) + 1)*exp(s**2/4)/s, 0, True)

    assert LT(sin(a*t), t, s) == (a/(a**2 + s**2), 0, True)
    assert LT(cos(a*t), t, s) == (s/(a**2 + s**2), 0, True)
    # TODO would be nice to have these come out better
    assert LT(exp(-a*t)*sin(b*t), t, s) == (1/b/(1 + (a + s)**2/b**2), -a, True)
    assert LT(exp(-a*t)*cos(b*t), t, s) == \
           (sqrt((a + s)**2)/b**2/(1 + (a + s)**2/b**2), -a, True)
    # TODO sinh, cosh have delicate cancellation

    # TODO conditions are a mess
    assert LT(besselj(0, t), t, s, noconds=True) == 1/sqrt(1 + s**2)
    assert LT(besselj(1, t), t, s, noconds=True) == 1 - 1/sqrt(1 + 1/s**2)

 def _eval_simplify(self, ratio, measure, rational, inverse):
     from sympy.simplify.hyperexpand import hyperexpand
     return hyperexpand(self)

def test_mellin_transform():
    from sympy import Max, Min, Ne
    MT = mellin_transform

    bpos = symbols('b', positive=True)

    # 8.4.2
    assert MT(x**nu*Heaviside(x - 1), x, s) == \
        (-1/(nu + s), (-oo, -re(nu)), True)
    assert MT(x**nu*Heaviside(1 - x), x, s) == \
        (1/(nu + s), (-re(nu), oo), True)

    assert MT((1 - x)**(beta - 1)*Heaviside(1 - x), x, s) == \
        (gamma(beta)*gamma(s)/gamma(beta + s), (0, oo), re(-beta) < 0)
    assert MT((x - 1)**(beta - 1)*Heaviside(x - 1), x, s) == \
        (gamma(beta)*gamma(1 - beta - s)/gamma(1 - s),
            (-oo, -re(beta) + 1), re(-beta) < 0)

    assert MT((1 + x)**(-rho), x, s) == \
        (gamma(s)*gamma(rho - s)/gamma(rho), (0, re(rho)), True)

    # TODO also the conditions should be simplified
    assert MT(abs(1 - x)**(-rho), x, s) == (
        cos(pi*(rho/2 - s))*gamma(s)*gamma(rho - s)/(cos(pi*rho/2)*gamma(rho)),
        (0, re(rho)), And(re(rho) - 1 < 0, re(rho) < 1))
    mt = MT((1 - x)**(beta - 1)*Heaviside(1 - x)
            + a*(x - 1)**(beta - 1)*Heaviside(x - 1), x, s)
    assert mt[1], mt[2] == ((0, -re(beta) + 1), True)

    assert MT((x**a - b**a)/(x - b), x, s)[0] == \
        pi*b**(a + s - 1)*sin(pi*a)/(sin(pi*s)*sin(pi*(a + s)))
    assert MT((x**a - bpos**a)/(x - bpos), x, s) == \
        (pi*bpos**(a + s - 1)*sin(pi*a)/(sin(pi*s)*sin(pi*(a + s))),
            (Max(-re(a), 0), Min(1 - re(a), 1)), True)

    expr = (sqrt(x + b**2) + b)**a
    assert MT(expr.subs(b, bpos), x, s) == \
        (-a*(2*bpos)**(a + 2*s)*gamma(s)*gamma(-a - 2*s)/gamma(-a - s + 1),
         (0, -re(a)/2), True)
    expr = (sqrt(x + b**2) + b)**a/sqrt(x + b**2)
    assert MT(expr.subs(b, bpos), x, s) == \
        (2**(a + 2*s)*bpos**(a + 2*s - 1)*gamma(s)
                                         *gamma(1 - a - 2*s)/gamma(1 - a - s),
            (0, -re(a)/2 + S(1)/2), True)

    # 8.4.2
    assert MT(exp(-x), x, s) == (gamma(s), (0, oo), True)
    assert MT(exp(-1/x), x, s) == (gamma(-s), (-oo, 0), True)

    # 8.4.5
    assert MT(log(x)**4*Heaviside(1 - x), x, s) == (24/s**5, (0, oo), True)
    assert MT(log(x)**3*Heaviside(x - 1), x, s) == (6/s**4, (-oo, 0), True)
    assert MT(log(x + 1), x, s) == (pi/(s*sin(pi*s)), (-1, 0), True)
    assert MT(log(1/x + 1), x, s) == (pi/(s*sin(pi*s)), (0, 1), True)
    assert MT(log(abs(1 - x)), x, s) == (pi/(s*tan(pi*s)), (-1, 0), True)
    assert MT(log(abs(1 - 1/x)), x, s) == (pi/(s*tan(pi*s)), (0, 1), True)

    # TODO we cannot currently do these (needs summation of 3F2(-1))
    #      this also implies that they cannot be written as a single g-function
    #      (although this is possible)
    mt = MT(log(x)/(x + 1), x, s)
    assert mt[1:] == ((0, 1), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)
    mt = MT(log(x)**2/(x + 1), x, s)
    assert mt[1:] == ((0, 1), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)
    mt = MT(log(x)/(x + 1)**2, x, s)
    assert mt[1:] == ((0, 2), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)

    # 8.4.14
    assert MT(erf(sqrt(x)), x, s) == \
        (-gamma(s + S(1)/2)/(sqrt(pi)*s), (-S(1)/2, 0), True)

def test_to_meijerg():
    x = symbols('x')
    assert hyperexpand(expr_to_holonomic(sin(x)).to_meijerg()) == sin(x)
    assert hyperexpand(expr_to_holonomic(cos(x)).to_meijerg()) == cos(x)
    assert hyperexpand(expr_to_holonomic(exp(x)).to_meijerg()) == exp(x)
    assert hyperexpand(expr_to_holonomic(log(x)).to_meijerg()).simplify() == log(x)
    assert expr_to_holonomic(4*x**2/3 + 7).to_meijerg() == 4*x**2/3 + 7
    assert hyperexpand(expr_to_holonomic(besselj(2, x), lenics=3).to_meijerg()) == besselj(2, x)
    p = hyper((-S(1)/2, -3), (), x)
    assert from_hyper(p).to_meijerg() == hyperexpand(p)
    p = hyper((S(1), S(3)), (S(2), ), x)
    assert (hyperexpand(from_hyper(p).to_meijerg()) - hyperexpand(p)).expand() == 0
    p = from_hyper(hyper((-2, -3), (S(1)/2, ), x))
    s = hyperexpand(hyper((-2, -3), (S(1)/2, ), x))
    C_0 = Symbol('C_0')
    C_1 = Symbol('C_1')
    D_0 = Symbol('D_0')
    assert (hyperexpand(p.to_meijerg()).subs({C_0:1, D_0:0}) - s).simplify() == 0
    p.y0 = {0: [1], S(1)/2: [0]}
    assert (hyperexpand(p.to_meijerg()) - s).simplify() == 0
    p = expr_to_holonomic(besselj(S(1)/2, x), initcond=False)
    assert (p.to_expr() - (D_0*sin(x) + C_0*cos(x) + C_1*sin(x))/sqrt(x)).simplify() == 0
    p = expr_to_holonomic(besselj(S(1)/2, x), y0={S(-1)/2: [sqrt(2)/sqrt(pi), sqrt(2)/sqrt(pi)]})
    assert (p.to_expr() - besselj(S(1)/2, x) - besselj(S(-1)/2, x)).simplify() == 0

def test_J14():
    p = hyper([S(1)/2, S(1)/2], [S(3)/2], z**2)
    assert hyperexpand(p) == asin(z)/z

def test_mellin_transform():
    from sympy import Max, Min
    MT = mellin_transform

    bpos = symbols('b', positive=True)

    # 8.4.2
    assert MT(x**nu*Heaviside(x - 1), x, s) \
           == (1/(-nu - s), (-oo, -re(nu)), True)
    assert MT(x**nu*Heaviside(1 - x), x, s) \
           == (1/(nu + s), (-re(nu), oo), True)

    assert MT((1-x)**(beta - 1)*Heaviside(1-x), x, s) \
           == (gamma(beta)*gamma(s)/gamma(beta + s),
               (0, oo), re(-beta) < 0)
    assert MT((x-1)**(beta - 1)*Heaviside(x-1), x, s) \
           == (gamma(beta)*gamma(1 - beta - s)/gamma(1 - s),
               (-oo, -re(beta) + 1), re(-beta) < 0)

    assert MT((1+x)**(-rho), x, s) == (gamma(s)*gamma(rho-s)/gamma(rho),
                                       (0, re(rho)), True)

    # TODO also the conditions should be simplified
    assert MT(abs(1-x)**(-rho), x, s) == \
        (cos(pi*rho/2 - pi*s)*gamma(s)*gamma(rho-s)/(cos(pi*rho/2)*gamma(rho)),\
         (0, re(rho)), And(re(rho) - 1 < 0, re(rho) < 1))

    mt = MT((1-x)**(beta-1)*Heaviside(1-x)
            + a*(x-1)**(beta-1)*Heaviside(x-1), x, s)
    assert mt[1], mt[2] == ((0, -re(beta) + 1), True)

    assert MT((x**a-b**a)/(x-b), x, s)[0] == \
           pi*b**(a+s-1)*sin(pi*a)/(sin(pi*s)*sin(pi*(a + s)))
    assert MT((x**a-bpos**a)/(x-bpos), x, s) == \
           (pi*bpos**(a+s-1)*sin(pi*a)/(sin(pi*s)*sin(pi*(a + s))),
            (Max(-re(a), 0), Min(1 - re(a), 1)), True)

    expr = (sqrt(x+b**2)+b)**a/sqrt(x+b**2)
    assert MT(expr.subs(b, bpos), x, s) == \
           (2**(a + 2*s)*bpos**(a + 2*s - 1)*gamma(s) \
                                         *gamma(1 - a - 2*s)/gamma(1 - a - s),
            (0, -re(a)/2 + S(1)/2), True)
    # TODO does not work with bneg, argument wrong. Needs changes to matching.
    #assert MT(expr.subs(b, -bpos), x, s) == \
    #       ((-1)**(a+1)*2**(a + 2*s)*bpos**(a + 2*s - 1)*gamma(a + s) \
    #                                      *gamma(1 - a - 2*s)/gamma(1 - s),
    #        (-re(a), -re(a)/2 + S(1)/2), True)
    expr = (sqrt(x+b**2)+b)**a
    assert MT(expr.subs(b, bpos), x, s) == \
           (-2**(a + 2*s)*a*bpos**(a + 2*s)*gamma(s)*gamma(-a - 2*s)/gamma(-a - s + 1),
            (0, -re(a)/2), True)
    #assert MT(expr.subs(b, -bpos), x, s) == \
    #       (2**(a + 2*s)*a*bpos**(a + 2*s)*gamma(-a - 2*s)*gamma(a + s)/gamma(-s + 1),
    #        (-re(a), -re(a)/2), True)
    # Test exponent 1:
    #assert MT(expr.subs({b: -bpos, a:1}), x, s) == \
    #       (-bpos**(2*s + 1)*gamma(s)*gamma(-s - S(1)/2)/(2*sqrt(pi)),
    #        (-1, -S(1)/2), True)

    # 8.4.2
    assert MT(exp(-x), x, s) == (gamma(s), (0, oo), True)
    assert MT(exp(-1/x), x, s) == (gamma(-s), (-oo, 0), True)

    # 8.4.5
    assert MT(log(x)**4*Heaviside(1-x), x, s) == (24/s**5, (0, oo), True)
    assert MT(log(x)**3*Heaviside(x-1), x, s) == (6/s**4, (-oo, 0), True)
    assert MT(log(x + 1), x, s) == (pi/(s*sin(pi*s)), (-1, 0), True)
    assert MT(log(1/x + 1), x, s) == (pi/(s*sin(pi*s)), (0, 1), True)
    assert MT(log(abs(1 - x)), x, s) == (pi/(s*tan(pi*s)), (-1, 0), True)
    assert MT(log(abs(1 - 1/x)), x, s) == (pi/(s*tan(pi*s)), (0, 1), True)

    # TODO we cannot currently do these (needs summation of 3F2(-1))
    #      this also implies that they cannot be written as a single g-function
    #      (although this is possible)
    mt = MT(log(x)/(x+1), x, s)
    assert mt[1:] == ((0, 1), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)
    mt = MT(log(x)**2/(x+1), x, s)
    assert mt[1:] == ((0, 1), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)
    mt = MT(log(x)/(x+1)**2, x, s)
    assert mt[1:] == ((0, 2), True)
    assert not hyperexpand(mt[0], allow_hyper=True).has(meijerg)

    # 8.4.14
    assert MT(erf(sqrt(x)), x, s) == \
           (-gamma(s + S(1)/2)/(sqrt(pi)*s), (-S(1)/2, 0), True)

    # 8.4.19
    assert MT(besselj(a, 2*sqrt(x)), x, s) == \
           (gamma(a/2 + s)/gamma(a/2 - s + 1), (-re(a)/2, S(3)/4), True)
    assert MT(sin(sqrt(x))*besselj(a, sqrt(x)), x, s) == \
           (2**a*gamma(S(1)/2 - 2*s)*gamma((a+1)/2 + s) \
                / (gamma(1 - s- a/2)*gamma(1 + a - 2*s)),
            (-(re(a) + 1)/2, S(1)/4), True)
    assert MT(cos(sqrt(x))*besselj(a, sqrt(x)), x, s) == \
           (2**a*gamma(a/2 + s)*gamma(S(1)/2 - 2*s)
                / (gamma(S(1)/2 - s - a/2)*gamma(a - 2*s + 1)),
            (-re(a)/2, S(1)/4), True)
    assert MT(besselj(a, sqrt(x))**2, x, s) == \
#.........这里部分代码省略.........