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
python类sqrt()的实例源码
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
#--------------------------------------------------------------------------
