Python sympy 模块，exp() 实例源码

def transmutagen_cram_error(degree, t0, prec=200):
    from ..partialfrac import (thetas_alphas, thetas_alphas_to_expr_complex,
        customre, t)
    from ..cram import get_CRAM_from_cache

    from sympy import re, exp, nsolve, diff

    expr = get_CRAM_from_cache(degree, prec)
    thetas, alphas, alpha0 = thetas_alphas(expr, prec)
    part_frac = thetas_alphas_to_expr_complex(thetas, alphas, alpha0)
    part_frac = part_frac.replace(customre, re)


    E = part_frac - exp(-t)
    E = E.evalf(20)

    return E.subs(t, nsolve(diff(E, t), t0, prec=prec))

def can_do(ap, bq, numerical=True, div=1, lowerplane=False):
    from sympy import exp_polar, exp
    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)/div
    [a, b, c, d] = [2, -1, 3, 1]
    if lowerplane:
        [a, b, c, d] = [2, -2, 3, -1]
    return tn(
        hyper(ap, bq, z).subs(repl),
        r.replace(exp_polar, exp).subs(repl),
        z, a=a, b=b, c=c, d=d)

def test_hyperexpand_bases():
    assert hyperexpand(hyper([2], [a], z)) == \
        a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
        lowergamma(a - 1, z) - 1
    # TODO [a+1, a-S.Half], [2*a]
    assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
    assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
        -1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
    assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
        (-3*z + 3)/4/(z*sqrt(-z + 1)) \
        + (6*z - 3)*asin(sqrt(z))/(4*z**(S(3)/2))
    assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
        - asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
    assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
        sqrt(z)*(6*z/7 - S(6)/5)*atanh(sqrt(z)) \
        + (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
    assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
        -4*log(sqrt(-z + 1)/2 + S(1)/2)/z
    # TODO hyperexpand(hyper([a], [2*a + 1], z))
    # TODO [S.Half, a], [S(3)/2, a+1]
    assert hyperexpand(hyper([2], [b, 1], z)) == \
        z**(-b/2 + S(1)/2)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
        + z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
    # TODO [a], [a - S.Half, 2*a]

def test_quantum_fourier():
    assert QFT(0, 3).decompose() == \
        SwapGate(0, 2)*HadamardGate(0)*CGate((0,), PhaseGate(1)) * \
        HadamardGate(1)*CGate((0,), TGate(2))*CGate((1,), PhaseGate(2)) * \
        HadamardGate(2)

    assert IQFT(0, 3).decompose() == \
        HadamardGate(2)*CGate((1,), RkGate(2, -2))*CGate((0,), RkGate(2, -3)) * \
        HadamardGate(1)*CGate((0,), RkGate(1, -2))*HadamardGate(0)*SwapGate(0, 2)

    assert represent(QFT(0, 3), nqubits=3) == \
        Matrix([[exp(2*pi*I/8)**(i*j % 8)/sqrt(8) for i in range(8)] for j in range(8)])

    assert QFT(0, 4).decompose()  # non-trivial decomposition
    assert qapply(QFT(0, 3).decompose()*Qubit(0, 0, 0)).expand() == qapply(
        HadamardGate(0)*HadamardGate(1)*HadamardGate(2)*Qubit(0, 0, 0)
    ).expand()

def test_wavefunction():
    a = 1/Z
    R = {
        (1, 0): 2*sqrt(1/a**3) * exp(-r/a),
        (2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
        (2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
        (3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
        (1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
        (3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
        (1 - r/(6*a)) * r/a,
        (3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
        (4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
        (1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
        (4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
        (1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
        (4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
        (1 - r/(12*a)) * (r/a)**2,
        (4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
    }
    for n, l in R:
        assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0

def test_Znm():
    # http://en.wikipedia.org/wiki/Solid_harmonics#List_of_lowest_functions
    th, ph = Symbol("theta", real=True), Symbol("phi", real=True)
    from sympy.abc import n,m

    assert Znm(0, 0, th, ph) == Ynm(0, 0, th, ph)
    assert Znm(1, -1, th, ph) == (-sqrt(2)*I*(Ynm(1, 1, th, ph)
                                  - exp(-2*I*ph)*Ynm(1, 1, th, ph))/2)
    assert Znm(1, 0, th, ph) == Ynm(1, 0, th, ph)
    assert Znm(1, 1, th, ph) == (sqrt(2)*(Ynm(1, 1, th, ph)
                                 + exp(-2*I*ph)*Ynm(1, 1, th, ph))/2)
    assert Znm(0, 0, th, ph).expand(func=True) == 1/(2*sqrt(pi))
    assert Znm(1, -1, th, ph).expand(func=True) == (sqrt(3)*I*sqrt(-cos(th)**2 + 1)*exp(I*ph)/(4*sqrt(pi))
                                                    - sqrt(3)*I*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi)))
    assert Znm(1, 0, th, ph).expand(func=True) == sqrt(3)*cos(th)/(2*sqrt(pi))
    assert Znm(1, 1, th, ph).expand(func=True) == (-sqrt(3)*sqrt(-cos(th)**2 + 1)*exp(I*ph)/(4*sqrt(pi))
                                                   - sqrt(3)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi)))
    assert Znm(2, -1, th, ph).expand(func=True) == (sqrt(15)*I*sqrt(-cos(th)**2 + 1)*exp(I*ph)*cos(th)/(4*sqrt(pi))
                                                    - sqrt(15)*I*sqrt(-cos(th)**2 + 1)*exp(-I*ph)*cos(th)/(4*sqrt(pi)))
    assert Znm(2, 0, th, ph).expand(func=True) == 3*sqrt(5)*cos(th)**2/(4*sqrt(pi)) - sqrt(5)/(4*sqrt(pi))
    assert Znm(2, 1, th, ph).expand(func=True) == (-sqrt(15)*sqrt(-cos(th)**2 + 1)*exp(I*ph)*cos(th)/(4*sqrt(pi))
                                                   - sqrt(15)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)*cos(th)/(4*sqrt(pi)))

def eval(cls, nu, z):
        from sympy import (unpolarify, expand_mul, uppergamma, exp, gamma,
                           factorial)
        nu2 = unpolarify(nu)
        if nu != nu2:
            return expint(nu2, z)
        if nu.is_Integer and nu <= 0 or (not nu.is_Integer and (2*nu).is_Integer):
            return unpolarify(expand_mul(z**(nu - 1)*uppergamma(1 - nu, z)))

        # Extract branching information. This can be deduced from what is
        # explained in lowergamma.eval().
        z, n = z.extract_branch_factor()
        if n == 0:
            return
        if nu.is_integer:
            if (nu > 0) != True:
                return
            return expint(nu, z) \
                - 2*pi*I*n*(-1)**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)
        else:
            return (exp(2*I*pi*nu*n) - 1)*z**(nu - 1)*gamma(1 - nu) + expint(nu, z)

def _getunbranched(cls, ar):
        from sympy import exp_polar, log, polar_lift
        if ar.is_Mul:
            args = ar.args
        else:
            args = [ar]
        unbranched = 0
        for a in args:
            if not a.is_polar:
                unbranched += arg(a)
            elif a.func is exp_polar:
                unbranched += a.exp.as_real_imag()[1]
            elif a.is_Pow:
                re, im = a.exp.as_real_imag()
                unbranched += re*unbranched_argument(
                    a.base) + im*log(abs(a.base))
            elif a.func is polar_lift:
                unbranched += arg(a.args[0])
            else:
                return None
        return unbranched

def _mul_args(f):
    """
    Return a list ``L`` such that Mul(*L) == f.

    If f is not a Mul or Pow, L=[f].
    If f=g**n for an integer n, L=[g]*n.
    If f is a Mul, L comes from applying _mul_args to all factors of f.
    """
    args = Mul.make_args(f)
    gs = []
    for g in args:
        if g.is_Pow and g.exp.is_Integer:
            n = g.exp
            base = g.base
            if n < 0:
                n = -n
                base = 1/base
            gs += [base]*n
        else:
            gs.append(g)
    return gs

def test_as_integral():
    from sympy import Function, Integral
    f = Function('f')
    assert mellin_transform(f(x), x, s).rewrite('Integral') == \
        Integral(x**(s - 1)*f(x), (x, 0, oo))
    assert fourier_transform(f(x), x, s).rewrite('Integral') == \
        Integral(f(x)*exp(-2*I*pi*s*x), (x, -oo, oo))
    assert laplace_transform(f(x), x, s).rewrite('Integral') == \
        Integral(f(x)*exp(-s*x), (x, 0, oo))
    assert str(inverse_mellin_transform(f(s), s, x, (a, b)).rewrite('Integral')) \
        == "Integral(x**(-s)*f(s), (s, _c - oo*I, _c + oo*I))"
    assert str(inverse_laplace_transform(f(s), s, x).rewrite('Integral')) == \
        "Integral(f(s)*exp(s*x), (s, _c - oo*I, _c + oo*I))"
    assert inverse_fourier_transform(f(s), s, x).rewrite('Integral') == \
        Integral(f(s)*exp(2*I*pi*s*x), (s, -oo, oo))

# NOTE this is stuck in risch because meijerint cannot handle it

def test_recursive():
    from sympy import symbols, exp_polar, expand
    a, b, c = symbols('a b c', positive=True)
    r = exp(-(x - a)**2)*exp(-(x - b)**2)
    e = integrate(r, (x, 0, oo), meijerg=True)
    assert simplify(e.expand()) == (
        sqrt(2)*sqrt(pi)*(
        (erf(sqrt(2)*(a + b)/2) + 1)*exp(-a**2/2 + a*b - b**2/2))/4)
    e = integrate(exp(-(x - a)**2)*exp(-(x - b)**2)*exp(c*x), (x, 0, oo), meijerg=True)
    assert simplify(e) == (
        sqrt(2)*sqrt(pi)*(erf(sqrt(2)*(2*a + 2*b + c)/4) + 1)*exp(-a**2 - b**2
        + (2*a + 2*b + c)**2/8)/4)
    assert simplify(integrate(exp(-(x - a - b - c)**2), (x, 0, oo), meijerg=True)) == \
        sqrt(pi)/2*(1 + erf(a + b + c))
    assert simplify(integrate(exp(-(x + a + b + c)**2), (x, 0, oo), meijerg=True)) == \
        sqrt(pi)/2*(1 - erf(a + b + c))

def main():
    x = Symbol("x")
    a = Symbol("a")
    h = Symbol("h")

    show( limit(sqrt(x**2 - 5*x + 6) - x, x, oo), -Rational(5)/2 )

    show( limit(x*(sqrt(x**2 + 1) - x), x, oo), Rational(1)/2 )

    show( limit(x - sqrt3(x**3 - 1), x, oo), Rational(0) )

    show( limit(log(1 + exp(x))/x, x, -oo), Rational(0) )

    show( limit(log(1 + exp(x))/x, x, oo), Rational(1) )

    show( limit(sin(3*x)/x, x, 0), Rational(3) )

    show( limit(sin(5*x)/sin(2*x), x, 0), Rational(5)/2 )

    show( limit(((x - 1)/(x + 1))**x, x, oo), exp(-2))

def real(self, nested_scope=None):
        """Return the correspond floating point number."""
        op = self.children[0].name
        expr = self.children[1]
        dispatch = {
            'sin': sympy.sin,
            'cos': sympy.cos,
            'tan': sympy.tan,
            'asin': sympy.asin,
            'acos': sympy.acos,
            'atan': sympy.atan,
            'exp': sympy.exp,
            'ln': sympy.log,
            'sqrt': sympy.sqrt
        }
        if op in dispatch:
            arg = expr.real(nested_scope)
            return dispatch[op](arg)
        else:
            raise NodeException("internal error: undefined external")

def sym(self, nested_scope=None):
        """Return the corresponding symbolic expression."""
        op = self.children[0].name
        expr = self.children[1]
        dispatch = {
            'sin': sympy.sin,
            'cos': sympy.cos,
            'tan': sympy.tan,
            'asin': sympy.asin,
            'acos': sympy.acos,
            'atan': sympy.atan,
            'exp': sympy.exp,
            'ln': sympy.log,
            'sqrt': sympy.sqrt
        }
        if op in dispatch:
            arg = expr.sym(nested_scope)
            return dispatch[op](arg)
        else:
            raise NodeException("internal error: undefined external")

def _print_Pow(self, expr):
        PREC = sympy.printing.precedence.precedence(expr)
        if expr.exp == 0:
            return "1.0"
        elif expr.exp == 0.5:
            return 'sqrt(%s)' % self._print(expr.base)
        elif expr.exp == 1:
            return self.parenthesize(expr.base,PREC)
        elif expr.exp.is_constant and int(expr.exp) == expr.exp:
            base = self.parenthesize(expr.base,PREC)
            if expr.exp > 0:
                return "(" + "*".join(repeat(base,int(expr.exp))) + ")"
            else:
                return "(1.0/" + "/".join(repeat(base,int(expr.exp))) + ")"
        return 'pow(%s, %s)' % (self._print(expr.base),
                                self._print(expr.exp))

def __init__(self, space, nuclearHamiltonian, overRideDT=None):
        self.myHamiltonian = nuclearHamiltonian
        self.mySpace = space
        if overRideDT is None:
            self.dt = self.mySpace.dt
        else:
            self.dt = overRideDT

        #set up the kinetic propagator,
        kineticSurface = self.myHamiltonian.myKineticOperator.surface
        kineticSurface = np.exp(-1.0j *self.dt * kineticSurface / self.mySpace.hbar)
        self.myKineticOperator = momentumOperator(self.mySpace, kineticSurface)

        #set up the potential propagator
        potentialSurface = self.myHamiltonian.myPotentialOperator.surface
        potentialSurface = np.exp(-1.0j * self.dt * potentialSurface / self.mySpace.hbar)

        self.myPotentialOperator = positionOperator(self.mySpace, potentialSurface)

def __init__(self, space, nuclearHamiltonian, overRideDT=None):
        self.myHamiltonian = nuclearHamiltonian
        self.mySpace = space
        if overRideDT is None:
            self.dt = self.mySpace.dt
        else:
            self.dt = overRideDT

        #set up the kinetic propagator,
        kineticSurface = self.myHamiltonian.myKineticOperator.surface
        kineticSurface = np.exp(-1.0j * .5 *self.dt * kineticSurface / self.mySpace.hbar)
        self.myKineticOperator = momentumOperator(self.mySpace, kineticSurface)

        #set up the potential propagator
        potentialSurface = self.myHamiltonian.myPotentialOperator.surface
        potentialSurface = np.exp(-1.0j * self.dt * potentialSurface / self.mySpace.hbar)

        self.myPotentialOperator = positionOperator(self.mySpace, potentialSurface)

def __init__(self, space, nuclearHamiltonian, overRideDT=None):
        #raise Exception("DO NOT USE HIGH-ACCURACY PROPAGATOR; FUNCTIONALITY NOT CODED YET")
        self.myHamiltonian = nuclearHamiltonian
        self.mySpace = space
        if overRideDT is None:
            self.dt = self.mySpace.dt
        else:
            self.dt = overRideDT

        gamma = 1.0 / (2.0 - 2.0**(1.0/3.0))
        lam = -1.0j * self.dt / self.mySpace.hbar
        #set up the kinetic propagator,
        kineticSurfaceLambda = self.myHamiltonian.myKineticOperator.surface * lam * .5
        kineticSurface1 = np.exp( gamma * kineticSurfaceLambda )
        kineticSurface2 = np.exp( (1.0 - gamma) * kineticSurfaceLambda )
        self.myKineticOperator1 = momentumOperator(self.mySpace, kineticSurface1)
        self.myKineticOperator2 = momentumOperator(self.mySpace, kineticSurface2)

        #set up the potential propagator
        potentialSurfaceLambda = self.myHamiltonian.myPotentialOperator.surface * lam
        potentialSurface1 = np.exp( gamma * potentialSurfaceLambda )
        potentialSurface2 = np.exp( (1.0 - 2.0 *gamma) * potentialSurfaceLambda )
        self.myPotentialOperator1 = positionOperator(self.mySpace, potentialSurface1)
        self.myPotentialOperator2 = positionOperator(self.mySpace, potentialSurface2)

def potentialFunction(self, x):
        naturalPotential = self.De * (1 - np.exp(-self.a * (x - self.center)))**2 + self.energyOffset
        imaginaryPotential = 0
        try:
            imaginaryPotentialZeros = np.zeros(x.shape)
        except:
            if x > self.startOfAbsorbingPotential:
                imaginaryPotential = -1.0j *self.strength * np.cosh( (np.real(x) - self.mySpace.xMax)**2 / (self.mySpace.Dx*30)**2)**(-2)
            else:
                imaginaryPotential = 0
            return naturalPotential + imaginaryPotential

        if self.absorbingPotential:
            imaginaryPotential = -1.0j *self.strength * np.cosh( (np.real(x) - self.mySpace.xMax)**2 / (self.mySpace.Dx*30)**2)**(-2)
            imaginaryPotential = np.where(x > self.startOfAbsorbingPotential, imaginaryPotential, imaginaryPotentialZeros)
        return naturalPotential + imaginaryPotential

def can_do(ap, bq, numerical=True, div=1, lowerplane=False):
    from sympy import exp_polar, exp
    r = hyperexpand(hyper(ap, bq, z))
    if r.has(hyper):
        return False
    if not numerical:
        return True
    repl = {}
    randsyms = r.free_symbols - {z}
    while randsyms:
        # Only randomly generated parameters are checked.
        for n, a in enumerate(randsyms):
            repl[a] = randcplx(n)/div
        if not any([b.is_Integer and b <= 0 for b in Tuple(*bq).subs(repl)]):
            break
    [a, b, c, d] = [2, -1, 3, 1]
    if lowerplane:
        [a, b, c, d] = [2, -2, 3, -1]
    return tn(
        hyper(ap, bq, z).subs(repl),
        r.replace(exp_polar, exp).subs(repl),
        z, a=a, b=b, c=c, d=d)

def test_hyperexpand_bases():
    assert hyperexpand(hyper([2], [a], z)) == \
        a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
        lowergamma(a - 1, z) - 1
    # TODO [a+1, a-S.Half], [2*a]
    assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
    assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
        -1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
    assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
        (-3*z + 3)/4/(z*sqrt(-z + 1)) \
        + (6*z - 3)*asin(sqrt(z))/(4*z**(S(3)/2))
    assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
        - asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
    assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
        sqrt(z)*(6*z/7 - S(6)/5)*atanh(sqrt(z)) \
        + (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
    assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
        -4*log(sqrt(-z + 1)/2 + S(1)/2)/z
    # TODO hyperexpand(hyper([a], [2*a + 1], z))
    # TODO [S.Half, a], [S(3)/2, a+1]
    assert hyperexpand(hyper([2], [b, 1], z)) == \
        z**(-b/2 + S(1)/2)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
        + z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
    # TODO [a], [a - S.Half, 2*a]

def test_meijerg_lookup():
    from sympy import uppergamma, Si, Ci
    assert hyperexpand(meijerg([a], [], [b, a], [], z)) == \
        z**b*exp(z)*gamma(-a + b + 1)*uppergamma(a - b, z)
    assert hyperexpand(meijerg([0], [], [0, 0], [], z)) == \
        exp(z)*uppergamma(0, z)
    assert can_do_meijer([a], [], [b, a + 1], [])
    assert can_do_meijer([a], [], [b + 2, a], [])
    assert can_do_meijer([a], [], [b - 2, a], [])

    assert hyperexpand(meijerg([a], [], [a, a, a - S(1)/2], [], z)) == \
        -sqrt(pi)*z**(a - S(1)/2)*(2*cos(2*sqrt(z))*(Si(2*sqrt(z)) - pi/2)
                                   - 2*sin(2*sqrt(z))*Ci(2*sqrt(z))) == \
        hyperexpand(meijerg([a], [], [a, a - S(1)/2, a], [], z)) == \
        hyperexpand(meijerg([a], [], [a - S(1)/2, a, a], [], z))
    assert can_do_meijer([a - 1], [], [a + 2, a - S(3)/2, a + 1], [])

def _eval_power(self, expt):
        """
        b is I = sqrt(-1)
        e is symbolic object but not equal to 0, 1

        I**r -> (-1)**(r/2) -> exp(r/2*Pi*I) -> sin(Pi*r/2) + cos(Pi*r/2)*I, r is decimal
        I**0 mod 4 -> 1
        I**1 mod 4 -> I
        I**2 mod 4 -> -1
        I**3 mod 4 -> -I
        """

        if isinstance(expt, Number):
            if isinstance(expt, Integer):
                expt = expt.p % 4
                if expt == 0:
                    return S.One
                if expt == 1:
                    return S.ImaginaryUnit
                if expt == 2:
                    return -S.One
                return -S.ImaginaryUnit
            return (S.NegativeOne)**(expt*S.Half)
        return

def test_quantum_fourier():
    assert QFT(0, 3).decompose() == \
        SwapGate(0, 2)*HadamardGate(0)*CGate((0,), PhaseGate(1)) * \
        HadamardGate(1)*CGate((0,), TGate(2))*CGate((1,), PhaseGate(2)) * \
        HadamardGate(2)

    assert IQFT(0, 3).decompose() == \
        HadamardGate(2)*CGate((1,), RkGate(2, -2))*CGate((0,), RkGate(2, -3)) * \
        HadamardGate(1)*CGate((0,), RkGate(1, -2))*HadamardGate(0)*SwapGate(0, 2)

    assert represent(QFT(0, 3), nqubits=3) == \
        Matrix([[exp(2*pi*I/8)**(i*j % 8)/sqrt(8) for i in range(8)] for j in range(8)])

    assert QFT(0, 4).decompose()  # non-trivial decomposition
    assert qapply(QFT(0, 3).decompose()*Qubit(0, 0, 0)).expand() == qapply(
        HadamardGate(0)*HadamardGate(1)*HadamardGate(2)*Qubit(0, 0, 0)
    ).expand()

def test_wavefunction():
    a = 1/Z
    R = {
        (1, 0): 2*sqrt(1/a**3) * exp(-r/a),
        (2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
        (2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
        (3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
        (1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
        (3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
        (1 - r/(6*a)) * r/a,
        (3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
        (4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
        (1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
        (4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
        (1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
        (4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
        (1 - r/(12*a)) * (r/a)**2,
        (4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
    }
    for n, l in R:
        assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0

def chapman(z, Nm, Hm, H_O, exp=NP.exp):
    """
    Return Chapman function electron density at height *z*, maximum
    electron density at the F-peak *Nm*, height at the maximum *Hm*,
    and the scale height of atomic oxygen *H_O*. The exponential
    function can be overridden with *exp*.
    """
    return Nm * exp((1 - ((z - Hm) / H_O) - exp(-(z - Hm) / H_O)) / 2)

def chapman_sym(z, Nm, Hm, H_O):
    """
    Return symbolic (i.e., :module:`sympy`) Chapman electron density
    profile function.
    """
    return chapman(z, Nm, Hm, H_O, exp=SYM.exp)

def paper_cram_error(degree, t0, prec=200):
    from ..partialfrac import t, customre

    from sympy import re, exp, nsolve, diff

    paper_part_frac = get_paper_part_frac(degree).replace(customre, re)

    E = paper_part_frac - exp(-t)

    return E.subs(t, nsolve(diff(E, t), t0, prec=prec))

def _plot(args):
    n = args.degree
    from sympy import exp

    from .. import plot_in_terminal

    rat_func = create_expression(n)

    print(rat_func)

    plot_in_terminal(rat_func - exp(-t), (0, 100), prec=20)

def main():
    setup_matplotlib_rc()

    expr = get_CRAM_from_cache(degree, prec)

    c = 0.6*degree

    # Get the translated approximation on [-1, 1]. This is similar logic from CRAM_exp().
    n, d = map(Poly, fraction(expr))
    inv = -c*(t + 1)/(t - 1)
    p, q = map(lambda i: Poly(i, t), fraction(inv))
    n, d = n.transform(p, q), d.transform(p, q)
    rat_func = n/d.TC()/(d/d.TC())
    rat_func = rat_func.evalf(prec)

    plt.clf()
    fig, ax = plt.subplots()

    fig.set_size_inches(4, 4)

    plot_in_terminal(expr - exp(-t), (0, 100), prec=prec, points=points,
        axes=ax)

    # Hide grid lines
    ax.grid(False)

    # Hide axes ticks
    ax.set_xticks([])
    ax.set_yticks([])

    plt.axis('off')
    plt.tight_layout()
    plt.savefig('logo.png', transparent=True, bbox_inches='tight', pad_inches=0)

def ContinuousRV(symbol, density, set=Interval(-oo, oo)):
    """
    Create a Continuous Random Variable given the following:

    -- a symbol
    -- a probability density function
    -- set on which the pdf is valid (defaults to entire real line)

    Returns a RandomSymbol.

    Many common continuous random variable types are already implemented.
    This function should be necessary only very rarely.

    Examples
    ========

    >>> from sympy import Symbol, sqrt, exp, pi
    >>> from sympy.stats import ContinuousRV, P, E

    >>> x = Symbol("x")

    >>> pdf = sqrt(2)*exp(-x**2/2)/(2*sqrt(pi)) # Normal distribution
    >>> X = ContinuousRV(x, pdf)

    >>> E(X)
    0
    >>> P(X>0)
    1/2
    """
    pdf = Lambda(symbol, density)
    dist = ContinuousDistributionHandmade(pdf, set)
    return SingleContinuousPSpace(symbol, dist).value

def pdf(self, x):
        alpha, beta, sigma = self.alpha, self.beta, self.sigma
        return (exp(-alpha*log(x/sigma) - beta*log(x/sigma)**2)
               *(alpha/x + 2*beta*log(x/sigma)/x))

def pdf(self, x):
        return 2**(1 - self.k/2)*x**(self.k - 1)*exp(-x**2/2)/gamma(self.k/2)

def pdf(self, x):
        k, l = self.k, self.l
        return exp(-(x**2+l**2)/2)*x**k*l / (l*x)**(k/2) * besseli(k/2-1, l*x)

def pdf(self, x):
        k = self.k
        return 1/(2**(k/2)*gamma(k/2))*x**(k/2 - 1)*exp(-x/2)

def pdf(self, x):
        return self.rate * exp(-self.rate*x)

def pdf(self, x):
        d1, d2 = self.d1, self.d2
        return (2*d1**(d1/2)*d2**(d2/2) / beta_fn(d1/2, d2/2) *
               exp(d1*x) / (d1*exp(2*x)+d2)**((d1+d2)/2))

def pdf(self, x):
        a, s, m = self.a, self.s, self.m
        return a/s * ((x-m)/s)**(-1-a) * exp(-((x-m)/s)**(-a))

def pdf(self, x):
        k, theta = self.k, self.theta
        return x**(k - 1) * exp(-x/theta) / (gamma(k)*theta**k)

def pdf(self, x):
        a, b = self.a, self.b
        return b**a/gamma(a) * x**(-a-1) * exp(-b/x)

def pdf(self, x):
        mu, b = self.mu, self.b
        return 1/(2*b)*exp(-Abs(x - mu)/b)

def pdf(self, x):
        mean, std = self.mean, self.std
        return exp(-(log(x) - mean)**2 / (2*std**2)) / (x*sqrt(2*pi)*std)

def pdf(self, x):
        a = self.a
        return sqrt(2/pi)*x**2*exp(-x**2/(2*a**2))/a**3

def pdf(self, x):
        mu, omega = self.mu, self.omega
        return 2*mu**mu/(gamma(mu)*omega**mu)*x**(2*mu - 1)*exp(-mu/omega*x**2)

def pdf(self, x):
        sigma = self.sigma
        return x/sigma**2*exp(-x**2/(2*sigma**2))