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

def test_clambdify():
    x = Symbol('x')
    y = Symbol('y')
    z = Symbol('z')
    f1 = sqrt(x*y)
    pf1 = lambdify((x, y), f1, 'math')
    cf1 = clambdify((x, y), f1)
    for i in xrange(10):
        if cf1(i, 10 - i) != pf1(i, 10 - i):
            raise ValueError("Values should be equal")
    f2 = (x - y) / z * pi
    pf2 = lambdify((x, y, z), f2, 'math')
    cf2 = clambdify((x, y, z), f2)
    if round(pf2(1, 2, 3), 14) != round(cf2(1, 2, 3), 14):
        raise ValueError("Values should be equal")
    # FIXME: slight difference in precision

def __init__(self, n, index):
        self.dim = n

        if index == '1-n':
            self.degree = 3
            data = [
                (fr(1, 2*n), fsd(n, (sqrt(fr(n, 3)), 1))),
                ]
        else:
            assert index == '2-n'
            self.degree = 5
            r = sqrt(fr(3, 5))
            data = [
                (fr(25*n**2 - 115*n + 162, 162), z(n)),
                (fr(70 - 25*n, 162), fsd(n, (r, 1))),
                (fr(25, 324), fsd(n, (r, 2))),
                ]

        self.points, self.weights = untangle(data)
        reference_volume = 2**n
        self.weights *= reference_volume
        return

def __init__(self, n):
        self.degree = 5

        r = sqrt(fr(7, 15))
        s, t = [sqrt((7 + i*sqrt(24)) / 15) for i in [+1, -1]]
        data = [
            (fr(5*n**2 - 15*n+14, 14), z(n)),
            (fr(25, 168), _s2(n, +r)),
            (fr(25, 168), _s2(n, -r)),
            (fr(-25*(n-2), 168), fsd(n, (r, 1))),
            (fr(5, 48), _s11(n, +s, -t)),
            (fr(5, 48), _s11(n, -s, +t)),
            (fr(-5*(n-2), 48), fsd(n, (s, 1))),
            (fr(-5*(n-2), 48), fsd(n, (t, 1))),
            ]

        self.points, self.weights = untangle(data)
        reference_volume = 2**n
        self.weights *= reference_volume
        return

def __init__(self):
        self.name = 'Maxwell'
        self.degree = 7

        r = sqrt(fr(12, 35))
        s, t = [sqrt((93 + i*3*sqrt(186)) / 155) for i in [+1, -1]]

        data = [
            (fr(1, 81), _z()),
            (fr(49, 324), _symm_r_0(r)),
            # ERR typo in Stroud: 648 vs 649
            (fr(31, 648), _symm_s_t(s, t))
            ]

        self.points, self.weights = untangle(data)
        self.weights *= 4
        return

def __init__(self):
        self.name = 'Phillips'

        c = 3*sqrt(385)
        r, s = [sqrt((105 + i*c) / 140) for i in [+1, -1]]
        t = sqrt(fr(3, 5))

        B1, B2 = [(77 - i*c) / 891 for i in [+1, -1]]
        B3 = fr(25, 324)

        self.degree = 7
        data = [
            (B1, _symm_r_0(r)),
            (B2, _symm_r_0(s)),
            (B3, _pm2(t, t))
            ]

        self.points, self.weights = untangle(data)
        self.weights *= 4
        return

def x():
    sqrt130 = sqrt(130)

    nu = sqrt((720 - 24*sqrt130) / 11)
    xi = sqrt(288 + 24*sqrt130)
    eta = sqrt((-216 + 24*sqrt130) / 7)
    A = (5175 - 13*sqrt130) / 8820
    B = (3870 + 283*sqrt130) / 493920
    C = (3204 - 281*sqrt130) / 197568
    # ERR in Stroud's book: 917568 vs. 197568
    D = (4239 + 373*sqrt130) / 197568

    data = [
        (A, numpy.array([[0, 0, 0]])),
        (B, fsd(3, (nu, 1))),
        (C, fsd(3, (xi, 2))),
        (D, pm(3, eta)),
        ]

    # ERR
    # TODO find out what's wrong
    warnings.warn('Stroud-Secrest\'s scheme X for E_3^r has degree 3, not 7.')
    return 3, data

def xi_():
    sqrt5 = sqrt(5)
    sqrt39 = sqrt(39)
    sqrt195 = sqrt(195)

    nu, xi = [
        sqrt(-50 + p_m*10*sqrt5 + 10*sqrt39 - p_m*2*sqrt195)
        for p_m in [+1, -1]
        ]
    eta = sqrt(36 + 4*sqrt39)
    mu, lmbda = [
        sqrt(54 + p_m * 18*sqrt5 + 6*sqrt39 + p_m * 2*sqrt195)
        for p_m in [+1, -1]
        ]
    A = (1725 - 26*sqrt39) / 2940
    B = (1065 + 171*sqrt39) / 54880
    C = (297 - 47*sqrt39) / 32928

    data = [
        (A, numpy.array([[0, 0, 0]])),
        (B, pm_roll(3, [xi, nu])),
        (C, pm(3, eta)),
        (C, pm_roll(3, [lmbda, mu])),
        ]
    return 7, data

def vii():
    # article:
    # nu, xi = numpy.sqrt((15 + plus_minus * 3*numpy.sqrt(5)))
    # A = 3/5
    # B = 1/30

    # book:
    nu, xi = [sqrt((5 - p_m * sqrt(5)) / 4) for p_m in [+1, -1]]
    A = fr(2, 5)
    B = fr(1, 20)

    data = [
        (A, numpy.array([[0, 0, 0]])),
        (B, pm_roll(3, [nu, xi])),
        ]
    return 5, data

def x(plus_minus):
    degree = 7

    sqrt15 = sqrt(15)

    r = sqrt((15 + plus_minus * sqrt15) / 4)
    s = sqrt((6 - plus_minus * sqrt15) / 2)
    t = sqrt((9 + plus_minus * 2*sqrt15) / 2)
    A = (720 + plus_minus * 8*sqrt15) / 2205
    B = (270 - plus_minus * 46*sqrt15) / 15435
    C = (162 + plus_minus * 41*sqrt15) / 6174
    D = (783 - plus_minus * 202*sqrt15) / 24696

    data = [
        (A, numpy.array([[0, 0, 0]])),
        (B, fsd(3, (r, 1))),
        (C, fsd(3, (s, 2))),
        (D, pm(3, t)),
        ]
    return degree, data

def xi_(p_m):
    degree = 7

    sqrt2 = sqrt(2)
    sqrt5 = sqrt(5)
    sqrt10 = sqrt(10)

    r = sqrt((25 + p_m * 15*sqrt2 + 5*sqrt5 + p_m * 3*sqrt10)/4)
    s = sqrt((25 + p_m * 15*sqrt2 - 5*sqrt5 - p_m * 3*sqrt10)/4)
    t = sqrt((3 - p_m * sqrt2) / 2)
    u = sqrt((9 - p_m * 3*sqrt2 - 3*sqrt5 + p_m * sqrt10) / 4)
    v = sqrt((9 - p_m * 3*sqrt2 + 3*sqrt5 - p_m * sqrt10) / 4)

    A = (80 + p_m * 8*sqrt2) / 245
    B = (395 - p_m * 279*sqrt2) / 13720
    C = (45 + p_m * 29*sqrt2) / 2744

    data = [
        (A, numpy.array([[0, 0, 0]])),
        (B, pm_roll(3, [r, s])),
        (C, pm_roll(3, [u, v])),
        (C, pm(3, t)),
        ]
    return degree, data

def __init__(self, n):
        self.degree = 2
        self.dim = n

        n2 = fr(n, 2) if n % 2 == 0 else fr(n-1, 2)
        pts = [[
            [sqrt(fr(2, n+2)) * cos(2*k*i*pi / (n+1)) for i in range(n+1)],
            [sqrt(fr(2, n+2)) * sin(2*k*i*pi / (n+1)) for i in range(n+1)],
            ] for k in range(1, n2+1)]
        if n % 2 == 1:
            sqrt3pm = numpy.full(n+1, 1/sqrt(n+2))
            sqrt3pm[1::2] *= -1
            pts.append(sqrt3pm)
        pts = numpy.vstack(pts).T

        data = [(fr(1, n+1), pts)]

        self.points, self.weights = untangle(data)
        self.weights *= volume_unit_ball(n)
        return

def __init__(self, degree):
        self.degree = degree
        if degree == 2:
            data = [
                (fr(1, 4), _s31((5 - sqrt(5))/20)),
                ]
        else:
            assert degree == 3
            data = [
                (-fr(4, 5), _s4()),
                (+fr(9, 20), _s31(fr(1, 6))),
                ]

        self.bary, self.weights = untangle(data)
        self.points = self.bary[:, 1:]
        return

def _gen5_3(n):
    '''Spherical product Lobatto formula.
    '''
    data = []
    s = sqrt(n+3)
    for k in range(1, n+1):
        rk = sqrt((k+2) * (n+3))
        Bk = fr(2**(k-n) * (n+1), (k+1) * (k+2) * (n+3))
        arr = [rk] + (n-k) * [s]
        data += [
            (Bk, pm_array0(n, arr, range(k-1, n)))
            ]
    B0 = 1 - sum([item[0]*len(item[1]) for item in data])
    data += [
        (B0, numpy.full((1, n), 0))
        ]
    return 5, data

def __init__(self, n):
        self.dim = n
        self.degree = 7

        r = 1
        s = sqrt(fr(1, n))
        t = sqrt(fr(1, 2))
        B = fr(8-n, n * (n+2) * (n+4))
        C = fr(n**3, 2**n * n * (n+2) * (n+4))
        D = fr(4, n * (n+2) * (n+4))

        data = [
            (B, fsd(n, (r, 1))),
            (C, pm(n, s)),
            # ERR Stroud's book wrongly states (t, t,..., t)_FS instead of
            # (t, t, 0, ..., 0)_FS.
            (D, fsd(n, (t, 2))),
            ]

        self.points, self.weights = untangle(data)
        self.weights *= integrate_monomial_over_unit_nsphere(n * [0])
        return

def vi():
    sqrt74255 = sqrt(74255)

    nu = sqrt(42)
    xi, eta = [sqrt((6615 - p_m * 21 * sqrt74255) / 454) for p_m in [+1, -1]]
    A = fr(5, 588)
    B, C = [
        (5272105 + p_m * 18733 * sqrt74255) / 43661940
        for p_m in [+1, -1]
        ]

    data = [
        (A, fsd(2, (nu, 1))),
        (B, pm(2, xi)),
        (C, pm(2, eta)),
        ]
    return 7, data

def standard_deviation(X, condition=None, **kwargs):
    """
    Standard Deviation of a random expression

    Square root of the Expectation of (X-E(X))**2

    Examples
    ========

    >>> from sympy.stats import Bernoulli, std
    >>> from sympy import Symbol, simplify

    >>> p = Symbol('p')
    >>> B = Bernoulli('B', p, 1, 0)

    >>> simplify(std(B))
    sqrt(p*(-p + 1))
    """
    return sqrt(variance(X, condition, **kwargs))

def test_sqrtdenest3():
    z = sqrt(13 - 2*r10 + 2*r2*sqrt(-2*r10 + 11))
    assert sqrtdenest(z) == -1 + r2 + r10
    assert sqrtdenest(z, max_iter=1) == -1 + sqrt(2) + sqrt(10)
    n = sqrt(2*r6/7 + 2*r7/7 + 2*sqrt(42)/7 + 2)
    d = sqrt(16 - 2*r29 + 2*sqrt(55 - 10*r29))
    assert sqrtdenest(n/d).equals(
        r7*(1 + r6 + r7)/(7*(sqrt(-2*r29 + 11) + r5)))
    z = sqrt(sqrt(r2 + 2) + 2)
    assert sqrtdenest(z) == z
    assert sqrtdenest(sqrt(-2*r10 + 4*r2*sqrt(-2*r10 + 11) + 20)) == \
        sqrt(-2*r10 - 4*r2 + 8*r5 + 20)
    assert sqrtdenest(sqrt((112 + 70*r2) + (46 + 34*r2)*r5)) == \
        r10 + 5 + 4*r2 + 3*r5
    z = sqrt(5 + sqrt(2*r6 + 5)*sqrt(-2*r29 + 2*sqrt(-10*r29 + 55) + 16))
    r = sqrt(-2*r29 + 11)
    assert sqrtdenest(z) == sqrt(r2*r + r3*r + r10 + r15 + 5)

def test_sqrt_symbolic_denest():
    x = Symbol('x')
    z = sqrt(((1 + sqrt(sqrt(2 + x) + 3))**2).expand())
    assert sqrtdenest(z) == sqrt((1 + sqrt(sqrt(2 + x) + 3))**2)
    z = sqrt(((1 + sqrt(sqrt(2 + cos(1)) + 3))**2).expand())
    assert sqrtdenest(z) == 1 + sqrt(sqrt(2 + cos(1)) + 3)
    z = ((1 + cos(2))**4 + 1).expand()
    assert sqrtdenest(z) == z
    z = sqrt(((1 + sqrt(sqrt(2 + cos(3*x)) + 3))**2 + 1).expand())
    assert sqrtdenest(z) == z
    c = cos(3)
    c2 = c**2
    assert sqrtdenest(sqrt(2*sqrt(1 + r3)*c + c2 + 1 + r3*c2)) == \
        -1 - sqrt(1 + r3)*c
    ra = sqrt(1 + r3)
    z = sqrt(20*ra*sqrt(3 + 3*r3) + 12*r3*ra*sqrt(3 + 3*r3) + 64*r3 + 112)
    assert sqrtdenest(z) == z

def test_special_assumptions():
    x = Symbol('x')
    z2 = z = Symbol('z', zero=True)
    assert z2 == z == S.Zero
    assert (2*z).is_positive is False
    assert (2*z).is_negative is False
    assert (2*z).is_zero is True
    assert (z2*z).is_positive is False
    assert (z2*z).is_negative is False
    assert (z2*z).is_zero is True

    e = -3 - sqrt(5) + (-sqrt(10)/2 - sqrt(2)/2)**2
    assert (e < 0) is S.false
    assert (e > 0) is S.false
    assert (e == 0) is False  # it's not a literal 0
    assert e.equals(0) is True

def test_pow_eval():
    # XXX Pow does not fully support conversion of negative numbers
    #     to their complex equivalent

    assert sqrt(-1) == I

    assert sqrt(-4) == 2*I
    assert sqrt( 4) == 2
    assert (8)**Rational(1, 3) == 2
    assert (-8)**Rational(1, 3) == 2*((-1)**Rational(1, 3))

    assert sqrt(-2) == I*sqrt(2)
    assert (-1)**Rational(1, 3) != I
    assert (-10)**Rational(1, 3) != I*((10)**Rational(1, 3))
    assert (-2)**Rational(1, 4) != (2)**Rational(1, 4)

    assert 64**Rational(1, 3) == 4
    assert 64**Rational(2, 3) == 16
    assert 24/sqrt(64) == 3
    assert (-27)**Rational(1, 3) == 3*(-1)**Rational(1, 3)

    assert (cos(2) / tan(2))**2 == (cos(2) / tan(2))**2

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_twave():
    A1, phi1, A2, phi2, f = symbols('A1, phi1, A2, phi2, f')
    c = Symbol('c')  # Speed of light in vacuum
    n = Symbol('n')  # Refractive index
    t = Symbol('t')  # Time
    w1 = TWave(A1, f, phi1)
    w2 = TWave(A2, f, phi2)
    assert w1.amplitude == A1
    assert w1.frequency == f
    assert w1.phase == phi1
    assert w1.wavelength == c/(f*n)
    assert w1.time_period == 1/f
    w3 = w1 + w2
    assert w3.amplitude == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)
    assert w3.frequency == f
    assert w3.wavelength == c/(f*n)
    assert w3.time_period == 1/f
    assert w3.angular_velocity == 2*pi*f
    assert w3.equation() == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)*cos(2*pi*f*t + phi1 + phi2)

def _apply_operator_Qubit(self, qubits, **options):
        """
        qubits: a set of qubits (Qubit)
        Returns: quantum object (quantum expression - QExpr)
        """
        if qubits.nqubits != self.nqubits:
            raise QuantumError(
                'WGate operates on %r qubits, got: %r'
                % (self.nqubits, qubits.nqubits)
            )

        # See 'Quantum Computer Science' by David Mermin p.92 -> W|a> result
        # Return (2/(sqrt(2^n)))|phi> - |a> where |a> is the current basis
        # state and phi is the superposition of basis states (see function
        # create_computational_basis above)
        basis_states = superposition_basis(self.nqubits)
        change_to_basis = (2/sqrt(2**self.nqubits))*basis_states
        return change_to_basis - qubits

def test_p():
    assert Px.hilbert_space == L2(Interval(S.NegativeInfinity, S.Infinity))
    assert qapply(Px*PxKet(px)) == px*PxKet(px)
    assert PxKet(px).dual_class() == PxBra
    assert PxBra(x).dual_class() == PxKet
    assert (Dagger(PxKet(py))*PxKet(px)).doit() == DiracDelta(px - py)
    assert (XBra(x)*PxKet(px)).doit() == \
        exp(I*x*px/hbar)/sqrt(2*pi*hbar)
    assert represent(PxKet(px)) == DiracDelta(px - px_1)

    rep_x = represent(PxOp(), basis=XOp)
    assert rep_x == -hbar*I*DiracDelta(x_1 - x_2)*DifferentialOperator(x_1)
    assert rep_x == represent(PxOp(), basis=XOp())
    assert rep_x == represent(PxOp(), basis=XKet)
    assert rep_x == represent(PxOp(), basis=XKet())

    assert represent(PxOp()*XKet(), basis=XKet) == \
        -hbar*I*DiracDelta(x - x_2)*DifferentialOperator(x)
    assert represent(XBra("y")*PxOp()*XKet(), basis=XKet) == \
        -hbar*I*DiracDelta(x - y)*DifferentialOperator(x)

def test_tensorproduct():
    a = BosonOp("a")
    b = BosonOp("b")
    ket1 = TensorProduct(BosonFockKet(1), BosonFockKet(2))
    ket2 = TensorProduct(BosonFockKet(0), BosonFockKet(0))
    ket3 = TensorProduct(BosonFockKet(0), BosonFockKet(2))
    bra1 = TensorProduct(BosonFockBra(0), BosonFockBra(0))
    bra2 = TensorProduct(BosonFockBra(1), BosonFockBra(2))
    assert qapply(TensorProduct(a, b ** 2) * ket1) == sqrt(2) * ket2
    assert qapply(TensorProduct(a, Dagger(b) * b) * ket1) == 2 * ket3
    assert qapply(bra1 * TensorProduct(a, b * b),
                  dagger=True) == sqrt(2) * bra2
    assert qapply(bra2 * ket1).doit() == TensorProduct(1, 1)
    assert qapply(TensorProduct(a, b * b) * ket1) == sqrt(2) * ket2
    assert qapply(Dagger(TensorProduct(a, b * b) * ket1),
                  dagger=True) == sqrt(2) * Dagger(ket2)

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_boson_states():
    a = BosonOp("a")

    # Fock states
    n = 3
    assert (BosonFockBra(0) * BosonFockKet(1)).doit() == 0
    assert (BosonFockBra(1) * BosonFockKet(1)).doit() == 1
    assert qapply(BosonFockBra(n) * Dagger(a)**n * BosonFockKet(0)) \
        == sqrt(prod(range(1, n+1)))

    # Coherent states
    alpha1, alpha2 = 1.2, 4.3
    assert (BosonCoherentBra(alpha1) * BosonCoherentKet(alpha1)).doit() == 1
    assert (BosonCoherentBra(alpha2) * BosonCoherentKet(alpha2)).doit() == 1
    assert abs((BosonCoherentBra(alpha1) * BosonCoherentKet(alpha2)).doit() -
               exp(-S(1) / 2 * (alpha1 - alpha2) ** 2)) < 1e-12
    assert qapply(a * BosonCoherentKet(alpha1)) == \
        alpha1 * BosonCoherentKet(alpha1)

def test_grover_iteration_2():
    numqubits = 4
    basis_states = superposition_basis(numqubits)
    v = OracleGate(numqubits, return_one_on_two)
    # After (pi/4)sqrt(pow(2, n)), IntQubit(2) should have highest prob
    # In this case, after around pi times (3 or 4)
    # print ''
    # print basis_states
    iterated = grover_iteration(basis_states, v)
    iterated = qapply(iterated)
    # print iterated
    iterated = grover_iteration(iterated, v)
    iterated = qapply(iterated)
    # print iterated
    iterated = grover_iteration(iterated, v)
    iterated = qapply(iterated)
    # print iterated
    # In this case, probability was highest after 3 iterations
    # Probability of Qubit('0010') was 251/256 (3) vs 781/1024 (4)
    # Ask about measurement
    expected = (-13*basis_states)/64 + 264*IntQubit(2, numqubits)/256
    assert qapply(expected) == iterated

def test_measure_partial():
    #Basic test of collapse of entangled two qubits (Bell States)
    state = Qubit('01') + Qubit('10')
    assert measure_partial(state, (0,)) == \
        [(Qubit('10'), Rational(1, 2)), (Qubit('01'), Rational(1, 2))]
    assert measure_partial(state, (0,)) == \
        measure_partial(state, (1,))[::-1]

    #Test of more complex collapse and probability calculation
    state1 = sqrt(2)/sqrt(3)*Qubit('00001') + 1/sqrt(3)*Qubit('11111')
    assert measure_partial(state1, (0,)) == \
        [(sqrt(2)/sqrt(3)*Qubit('00001') + 1/sqrt(3)*Qubit('11111'), 1)]
    assert measure_partial(state1, (1, 2)) == measure_partial(state1, (3, 4))
    assert measure_partial(state1, (1, 2, 3)) == \
        [(Qubit('00001'), Rational(2, 3)), (Qubit('11111'), Rational(1, 3))]

    #test of measuring multiple bits at once
    state2 = Qubit('1111') + Qubit('1101') + Qubit('1011') + Qubit('1000')
    assert measure_partial(state2, (0, 1, 3)) == \
        [(Qubit('1000'), Rational(1, 4)), (Qubit('1101'), Rational(1, 4)),
         (Qubit('1011')/sqrt(2) + Qubit('1111')/sqrt(2), Rational(1, 2))]
    assert measure_partial(state2, (0,)) == \
        [(Qubit('1000'), Rational(1, 4)),
         (Qubit('1111')/sqrt(3) + Qubit('1101')/sqrt(3) +
          Qubit('1011')/sqrt(3), Rational(3, 4))]

def _represent_NumberOp(self, basis, **options):
        ndim_info = options.get('ndim', 4)
        format = options.get('format','sympy')
        spmatrix = options.get('spmatrix', 'csr')
        matrix = matrix_zeros(ndim_info, ndim_info, **options)
        for i in range(ndim_info - 1):
            value = sqrt(i + 1)
            if format == 'scipy.sparse':
                value = float(value)
            matrix[i + 1, i] = value
        if format == 'scipy.sparse':
            matrix = matrix.tocsr()
        return matrix

    #--------------------------------------------------------------------------
    # Printing Methods
    #--------------------------------------------------------------------------

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_units():
    assert (5*m/s * day) / km == 432
    assert foot / meter == Rational('0.3048')
    # amu is a pure mass so mass/mass gives a number, not an amount (mol)
    assert str(grams/(amu).n(2)) == '6.0e+23'

    # Light from the sun needs about 8.3 minutes to reach earth
    t = (1*au / speed_of_light).evalf() / minute
    assert abs(t - 8.31) < 0.1

    assert sqrt(m**2) == m
    assert (sqrt(m))**2 == m

    t = Symbol('t')
    assert integrate(t*m/s, (t, 1*s, 5*s)) == 12*m*s
    assert (t * m/s).integrate((t, 1*s, 5*s)) == 12*m*s

def test_dot_different_frames():
    assert dot(N.x, A.x) == cos(q1)
    assert dot(N.x, A.y) == -sin(q1)
    assert dot(N.x, A.z) == 0
    assert dot(N.y, A.x) == sin(q1)
    assert dot(N.y, A.y) == cos(q1)
    assert dot(N.y, A.z) == 0
    assert dot(N.z, A.x) == 0
    assert dot(N.z, A.y) == 0
    assert dot(N.z, A.z) == 1

    assert dot(N.x, A.x + A.y) == sqrt(2)*cos(q1 + pi/4) == dot(A.x + A.y, N.x)

    assert dot(A.x, C.x) == cos(q3)
    assert dot(A.x, C.y) == 0
    assert dot(A.x, C.z) == sin(q3)
    assert dot(A.y, C.x) == sin(q2)*sin(q3)
    assert dot(A.y, C.y) == cos(q2)
    assert dot(A.y, C.z) == -sin(q2)*cos(q3)
    assert dot(A.z, C.x) == -cos(q2)*sin(q3)
    assert dot(A.z, C.y) == sin(q2)
    assert dot(A.z, C.z) == cos(q2)*cos(q3)

def real_imag(ba, bd, gen):
    """
    Helper function, to get the real and imaginary part of a rational function
    evaluated at sqrt(-1) without actually evaluating it at sqrt(-1)

    Separates the even and odd power terms by checking the degree of terms wrt
    mod 4. Returns a tuple (ba[0], ba[1], bd) where ba[0] is real part
    of the numerator ba[1] is the imaginary part and bd is the denominator
    of the rational function.
    """
    bd = bd.as_poly(gen).as_dict()
    ba = ba.as_poly(gen).as_dict()
    denom_real = [value if key[0] % 4 == 0 else -value if key[0] % 4 == 2 else 0 for key, value in bd.items()]
    denom_imag = [value if key[0] % 4 == 1 else -value if key[0] % 4 == 3 else 0 for key, value in bd.items()]
    bd_real = sum(r for r in denom_real)
    bd_imag = sum(r for r in denom_imag)
    num_real = [value if key[0] % 4 == 0 else -value if key[0] % 4 == 2 else 0 for key, value in ba.items()]
    num_imag = [value if key[0] % 4 == 1 else -value if key[0] % 4 == 3 else 0 for key, value in ba.items()]
    ba_real = sum(r for r in num_real)
    ba_imag = sum(r for r in num_imag)
    ba = ((ba_real*bd_real + ba_imag*bd_imag).as_poly(gen), (ba_imag*bd_real - ba_real*bd_imag).as_poly(gen))
    bd = (bd_real*bd_real + bd_imag*bd_imag).as_poly(gen)
    return (ba[0], ba[1], bd)

def test_Domain_unify_algebraic():
    sqrt5 = QQ.algebraic_field(sqrt(5))
    sqrt7 = QQ.algebraic_field(sqrt(7))
    sqrt57 = QQ.algebraic_field(sqrt(5), sqrt(7))

    assert sqrt5.unify(sqrt7) == sqrt57

    assert sqrt5.unify(sqrt5[x, y]) == sqrt5[x, y]
    assert sqrt5[x, y].unify(sqrt5) == sqrt5[x, y]

    assert sqrt5.unify(sqrt5.frac_field(x, y)) == sqrt5.frac_field(x, y)
    assert sqrt5.frac_field(x, y).unify(sqrt5) == sqrt5.frac_field(x, y)

    assert sqrt5.unify(sqrt7[x, y]) == sqrt57[x, y]
    assert sqrt5[x, y].unify(sqrt7) == sqrt57[x, y]

    assert sqrt5.unify(sqrt7.frac_field(x, y)) == sqrt57.frac_field(x, y)
    assert sqrt5.frac_field(x, y).unify(sqrt7) == sqrt57.frac_field(x, y)

def test_to_ZZ_ANP_poly():
    A = AlgebraicField(QQ, sqrt(2))
    R, x = ring("x", A)
    f = x*(sqrt(2) + 1)

    T, x_, z_ = ring("x_, z_", ZZ)
    f_ = x_*z_ + x_

    assert _to_ZZ_poly(f, T) == f_
    assert _to_ANP_poly(f_, R) == f

    R, x, t, s = ring("x, t, s", A)
    f = x*t**2 + x*s + sqrt(2)

    D, t_, s_ = ring("t_, s_", ZZ)
    T, x_, z_ = ring("x_, z_", D)
    f_ = (t_**2 + s_)*x_ + z_

    assert _to_ZZ_poly(f, T) == f_
    assert _to_ANP_poly(f_, R) == f

def test_PolyElement_from_expr():
    x, y, z = symbols("x,y,z")
    R, X, Y, Z = ring((x, y, z), ZZ)

    f = R.from_expr(1)
    assert f == 1 and isinstance(f, R.dtype)

    f = R.from_expr(x)
    assert f == X and isinstance(f, R.dtype)

    f = R.from_expr(x*y*z)
    assert f == X*Y*Z and isinstance(f, R.dtype)

    f = R.from_expr(x*y*z + x*y + x)
    assert f == X*Y*Z + X*Y + X and isinstance(f, R.dtype)

    f = R.from_expr(x**3*y*z + x**2*y**7 + 1)
    assert f == X**3*Y*Z + X**2*Y**7 + 1 and isinstance(f, R.dtype)

    raises(ValueError, lambda: R.from_expr(1/x))
    raises(ValueError, lambda: R.from_expr(2**x))
    raises(ValueError, lambda: R.from_expr(7*x + sqrt(2)))

def test_imps():
    # Here we check if the default returned functions are anonymous - in
    # the sense that we can have more than one function with the same name
    f = implemented_function('f', lambda x: 2*x)
    g = implemented_function('f', lambda x: math.sqrt(x))
    l1 = lambdify(x, f(x))
    l2 = lambdify(x, g(x))
    assert str(f(x)) == str(g(x))
    assert l1(3) == 6
    assert l2(3) == math.sqrt(3)
    # check that we can pass in a Function as input
    func = sympy.Function('myfunc')
    assert not hasattr(func, '_imp_')
    my_f = implemented_function(func, lambda x: 2*x)
    assert hasattr(func, '_imp_')
    # Error for functions with same name and different implementation
    f2 = implemented_function("f", lambda x: x + 101)
    raises(ValueError, lambda: lambdify(x, f(f2(x))))

def test_special_printers():
    class IntervalPrinter(LambdaPrinter):
        """Use ``lambda`` printer but print numbers as ``mpi`` intervals. """

        def _print_Integer(self, expr):
            return "mpi('%s')" % super(IntervalPrinter, self)._print_Integer(expr)

        def _print_Rational(self, expr):
            return "mpi('%s')" % super(IntervalPrinter, self)._print_Rational(expr)

    def intervalrepr(expr):
        return IntervalPrinter().doprint(expr)

    expr = sympy.sqrt(sympy.sqrt(2) + sympy.sqrt(3)) + sympy.S(1)/2

    func0 = lambdify((), expr, modules="mpmath", printer=intervalrepr)
    func1 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter)
    func2 = lambdify((), expr, modules="mpmath", printer=IntervalPrinter())

    mpi = type(mpmath.mpi(1, 2))

    assert isinstance(func0(), mpi)
    assert isinstance(func1(), mpi)
    assert isinstance(func2(), mpi)

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))