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
python类Rational()的实例源码
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 = '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 __init__(self):
self.name = 'Meister'
self.degree = 7
r = fr(2, 3)
s = fr(1, 3)
data = [
(fr(1024, 6720), _z()),
(fr(576, 6720), _symm_s(r)),
(fr(576, 6720), _symm_r_0(r)),
(-fr(9, 6720), _symm_s(s)),
(fr(117, 6720), _symm_s_t(1, s)),
(fr(47, 6720), _symm_s(1)),
]
self.points, self.weights = untangle(data)
self.weights *= 4
return
def __init__(self, index):
self.name = 'Irwin({})'.format(index)
if index == 1:
self.degree = 3
data = [
(fr(14, 48), _symm_s(1)),
(-fr(1, 48), _symm_s_t(3, 1))
]
else:
assert index == 2
self.degree = 5
data = [
(fr(889, 2880), _symm_s(1)),
(-fr(98, 2880), _symm_s_t(3, 1)),
(fr(5, 2880), _symm_s(3)),
(fr(11, 2880), _symm_s_t(5, 1)),
]
self.points, self.weights = untangle(data)
self.weights *= 4
return
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 __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 __init__(self, index):
if index == 4:
self.degree = 2
data = [
(fr(1, 4), _s31(0.1381966011250105))
]
else:
assert index == 5
self.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 __init__(self, index):
if index == 1:
self.degree = 2
data = [
(fr(1, 4), _r(1/sqrt(5))),
]
elif index == 2:
self.degree = 2
data = [
(fr(1, 4), _r(-1/sqrt(5))),
]
else:
assert index == 3
self.degree = 3
data = [
(-fr(4, 5), _s4()),
(fr(9, 20), _r(fr(1, 3))),
]
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 _xi11(d, a):
assert d > 1
b = fr(1 - (d-1) * a, 2)
if d == 2:
out = numpy.array([
[b, b, a],
[b, a, b],
[a, b, b],
])
else:
assert d == 3
out = numpy.array([
[b, b, a, a],
[b, a, b, a],
[b, a, a, b],
[a, b, a, b],
[a, a, b, b],
[a, b, b, a],
])
return out
def __init__(self, n):
self.dim = n
self.degree = 3
r = fr(1, n+1)
s = fr(1, n)
prod = (n+1) * (n+2) * (n+3)
A = fr((3-n) * (n+1)**3, prod)
B = fr(3, prod)
C = fr(n**3, prod)
data = [
(A, [(n+1) * [r]]),
(B, rd(n+1, [(1, 1)])),
(C, rd(n+1, [(s, n)])),
]
self.bary, self.weights = untangle(data)
self.points = self.bary[:, 1:]
return
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 integrate_monomial_over_unit_nsphere(alpha):
'''
Gerald B. Folland,
How to Integrate a Polynomial over a Sphere,
The American Mathematical Monthly,
Vol. 108, No. 5 (May, 2001), pp. 446-448,
<https://doi.org/10.2307/2695802>.
'''
if any(a % 2 == 1 for a in alpha):
return 0
# Use lgamma since other with ordinary gamma, numerator and denominator
# might overflow.
# return 2 * math.exp(
# math.fsum([math.lgamma(0.5*(a+1)) for a in alpha])
# - math.lgamma(math.fsum([0.5*(a+1) for a in alpha]))
# )
# Explicitly cast a to int (from numpy.int64) to work around bug
# <https://github.com/sympy/sympy/issues/13618>.
return 2 * (
prod([gamma(Rational(int(a)+1, 2)) for a in alpha])
/ gamma(sum([Rational(int(a)+1, 2) for a in alpha]))
)
def __init__(self, degree):
self.degree = degree
if degree == 2:
data = {
's2': [
[fr(1, 3), fr(1, 6)],
]
}
else:
assert degree == 3
data = {
's3': [
[-fr(27, 48)],
],
's2': [
[+fr(25, 48), fr(1, 5)],
]
}
self.bary, self.weights = untangle2(data)
self.points = self.bary[:, 1:]
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 dotest(s):
x = Symbol("x")
y = Symbol("y")
l = [
Rational(2),
Float("1.3"),
x,
y,
pow(x, y)*y,
5,
5.5
]
for x in l:
for y in l:
s(x, y)
return True
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 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 __pow__(self, other):
if self.data is None:
raise ValueError("No power without ndarray data.")
numpy = import_module('numpy')
free = self.free
marray = self.data
for metric in free:
marray = numpy.tensordot(
marray,
numpy.tensordot(
metric[0]._tensortype.data,
marray,
(1, 0)
),
(0, 0)
)
pow2 = marray[()]
return pow2 ** (Rational(1, 2) * other)
def test_legendre_poly():
raises(ValueError, lambda: legendre_poly(-1, x))
assert legendre_poly(1, x, polys=True) == Poly(x)
assert legendre_poly(0, x) == 1
assert legendre_poly(1, x) == x
assert legendre_poly(2, x) == Q(3, 2)*x**2 - Q(1, 2)
assert legendre_poly(3, x) == Q(5, 2)*x**3 - Q(3, 2)*x
assert legendre_poly(4, x) == Q(35, 8)*x**4 - Q(30, 8)*x**2 + Q(3, 8)
assert legendre_poly(5, x) == Q(63, 8)*x**5 - Q(70, 8)*x**3 + Q(15, 8)*x
assert legendre_poly(6, x) == Q(
231, 16)*x**6 - Q(315, 16)*x**4 + Q(105, 16)*x**2 - Q(5, 16)
assert legendre_poly(1).dummy_eq(x)
assert legendre_poly(1, polys=True) == Poly(x)
def test_laguerre_poly():
raises(ValueError, lambda: laguerre_poly(-1, x))
assert laguerre_poly(1, x, polys=True) == Poly(-x + 1)
assert laguerre_poly(0, x) == 1
assert laguerre_poly(1, x) == -x + 1
assert laguerre_poly(2, x) == Q(1, 2)*x**2 - Q(4, 2)*x + 1
assert laguerre_poly(3, x) == -Q(1, 6)*x**3 + Q(9, 6)*x**2 - Q(18, 6)*x + 1
assert laguerre_poly(4, x) == Q(
1, 24)*x**4 - Q(16, 24)*x**3 + Q(72, 24)*x**2 - Q(96, 24)*x + 1
assert laguerre_poly(5, x) == -Q(1, 120)*x**5 + Q(25, 120)*x**4 - Q(
200, 120)*x**3 + Q(600, 120)*x**2 - Q(600, 120)*x + 1
assert laguerre_poly(6, x) == Q(1, 720)*x**6 - Q(36, 720)*x**5 + Q(450, 720)*x**4 - Q(2400, 720)*x**3 + Q(5400, 720)*x**2 - Q(4320, 720)*x + 1
assert laguerre_poly(0, x, a) == 1
assert laguerre_poly(1, x, a) == -x + a + 1
assert laguerre_poly(2, x, a) == x**2/2 + (-a - 2)*x + a**2/2 + 3*a/2 + 1
assert laguerre_poly(3, x, a) == -x**3/6 + (a/2 + Q(
3)/2)*x**2 + (-a**2/2 - 5*a/2 - 3)*x + a**3/6 + a**2 + 11*a/6 + 1
assert laguerre_poly(1).dummy_eq(-x + 1)
assert laguerre_poly(1, polys=True) == Poly(-x + 1)
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 _print_Rational(self, e):
"""Print a Rational object.
:param e: The expression.
:rtype : bce.dom.mathml.all.Base
:return: The printed MathML object.
"""
assert isinstance(e, _sympy.Rational)
if e.q == 1:
# Don't do division if the denominator is 1.
return _mathml.NumberComponent(str(e.p))
return _mathml.FractionComponent(
_mathml.NumberComponent(str(e.p)),
_mathml.NumberComponent(str(e.q))
)
def convert_float_string_to_rational(symbol):
"""Convert a string that contains a float number to SymPy's internal Rational type.
:type symbol: str
:param symbol: The string.
:return: Converted value.
"""
# Find the position of decimal dot.
dot_pos = symbol.find(".")
if dot_pos < 0:
raise ValueError("Missing decimal dot.")
# Get the value of its numerator and denominator.
numerator = int(symbol[:dot_pos] + symbol[dot_pos + 1:])
denominator = 10 ** (len(symbol) - 1 - dot_pos)
return _sympy.Rational(numerator, denominator)
def dotest(s):
x = Symbol("x")
y = Symbol("y")
l = [
Rational(2),
Float("1.3"),
x,
y,
pow(x, y)*y,
5,
5.5
]
for x in l:
for y in l:
s(x, y)
return True
def as_coefficients_dict(self):
"""Return a dictionary mapping terms to their Rational coefficient.
Since the dictionary is a defaultdict, inquiries about terms which
were not present will return a coefficient of 0. If an expression is
not an Add it is considered to have a single term.
Examples
========
>>> from sympy.abc import a, x
>>> (3*x + a*x + 4).as_coefficients_dict()
{1: 4, x: 3, a*x: 1}
>>> _[a]
0
>>> (3*a*x).as_coefficients_dict()
{a*x: 3}
"""
c, m = self.as_coeff_Mul()
if not c.is_Rational:
c = S.One
m = self
d = defaultdict(int)
d.update({m: c})
return d
def primitive(self):
"""Return the positive Rational that can be extracted non-recursively
from every term of self (i.e., self is treated like an Add). This is
like the as_coeff_Mul() method but primitive always extracts a positive
Rational (never a negative or a Float).
Examples
========
>>> from sympy.abc import x
>>> (3*(x + 1)**2).primitive()
(3, (x + 1)**2)
>>> a = (6*x + 2); a.primitive()
(2, 3*x + 1)
>>> b = (x/2 + 3); b.primitive()
(1/2, x + 6)
>>> (a*b).primitive() == (1, a*b)
True
"""
if not self:
return S.One, S.Zero
c, r = self.as_coeff_Mul(rational=True)
if c.is_negative:
c, r = -c, -r
return c, r
def test_definition():
# want to test if the system can have several units of the same dimension
dm = Quantity("dm", length, Rational(1, 10))
base = (m, s)
base_dim = (m.dimension, s.dimension)
ms = UnitSystem(base, (c, dm), "MS", "MS system")
assert set(ms._base_units) == set(base)
assert set(ms._units) == set((m, s, c, dm))
#assert ms._units == DimensionSystem._sort_dims(base + (velocity,))
assert ms.name == "MS"
assert ms.descr == "MS system"
assert ms._system._base_dims == DimensionSystem.sort_dims(base_dim)
assert set(ms._system._dims) == set(base_dim + (velocity,))
def test_units():
assert convert_to((5*m/s * day) / km, 1) == 432
assert convert_to(foot / meter, meter) == Rational('0.3048')
# amu is a pure mass so mass/mass gives a number, not an amount (mol)
# TODO: need better simplification routine:
assert str(convert_to(grams/amu, grams).n(2)) == '6.0e+23'
# Light from the sun needs about 8.3 minutes to reach earth
t = (1*au / speed_of_light) / minute
# TODO: need a better way to simplify expressions containing units:
t = convert_to(convert_to(t, meter / minute), meter)
assert t == 49865956897/5995849160
# TODO: fix this, it should give `m` without `Abs`
assert sqrt(m**2) == Abs(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
