summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorDaniel Krenn <devel@danielkrenn.at>2014-12-28 10:33:16 +0100
committerDaniel Krenn <devel@danielkrenn.at>2014-12-28 10:33:16 +0100
commit53d60ea14b049384d1f4becba615e37f05530f43 (patch)
tree325196383a60e17110e2226ca4c3ce575f5adf60
parentcoercion from other FFPD rings implemented (diff)
parentmoved code like in "9914678 - moved static methods to parent" to allow merging (diff)
Merge branch 'u/chapoton/10519' into t/10519-on-6.5.beta4
* u/chapoton/10519: moved code like in "9914678 - moved static methods to parent" to allow merging trac #10519 pep8 fully compliant Conflicts: src/sage/combinat/asymptotics_multivariate_generating_functions.py
-rw-r--r--src/sage/combinat/asymptotics_multivariate_generating_functions.py524
1 files changed, 269 insertions, 255 deletions
diff --git a/src/sage/combinat/asymptotics_multivariate_generating_functions.py b/src/sage/combinat/asymptotics_multivariate_generating_functions.py
index 94d310e..fd1b119 100644
--- a/src/sage/combinat/asymptotics_multivariate_generating_functions.py
+++ b/src/sage/combinat/asymptotics_multivariate_generating_functions.py
@@ -102,7 +102,7 @@ Another smooth point example (Example 5.4 of [RaWi2008a]_)::
Computing derivatives of auxiallary functions...
Computing derivatives of more auxiliary functions...
Computing second order differential operator actions...
- sage: asy
+ sage: asy
(1/12*sqrt(3)*2^(2/3)*gamma(1/3)/(pi*r^(1/3)), 1,
1/12*sqrt(3)*2^(2/3)*gamma(1/3)/(pi*r^(1/3)))
sage: F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1])
@@ -184,7 +184,7 @@ from sage.misc.mrange import cartesian_product_iterator, mrange
from sage.modules.free_module_element import vector
from sage.rings.arith import binomial, xgcd
from sage.rings.all import CC
-from sage.rings.fraction_field import FractionField
+from sage.rings.fraction_field import FractionField
from sage.rings.integer import Integer
from sage.rings.rational_field import QQ
from sage.sets.set import Set
@@ -378,7 +378,7 @@ class FFPDElement(sage.structure.element.RingElement):
x^3*y^2 + 2*x^3*y + x^2*y^2 + x^3 - 2*x^2*y - x*y^2 - 3*x^2 - 2*x*y
- y^2 + 3*x + 2*y - 1
"""
- return prod([q**e for q, e in self.denominator_factored()])
+ return prod(q ** e for q, e in self.denominator_factored())
def __cmp__(self, other):
@@ -511,7 +511,7 @@ class FFPDElement(sage.structure.element.RingElement):
-e^y/(x^3*y^2 + 2*x^3*y + x^2*y^2 + x^3 - 2*x^2*y - x*y^2 - 3*x^2 -
2*x*y - y^2 + 3*x + 2*y - 1)
"""
- return self.numerator()/self.denominator()
+ return self.numerator() / self.denominator()
def __repr__(self):
@@ -668,9 +668,9 @@ class FFPDElement(sage.structure.element.RingElement):
sd = self.denominator()
od = other.denominator()
- return bool(len(sdf) < len(odf) or\
- (len(sdf) == len(odf) and sd < od) or\
- (len(sdf) == len(odf) and sd == od and sn < on))
+ return bool(len(sdf) < len(odf) or
+ (len(sdf) == len(odf) and sd < od) or
+ (len(sdf) == len(odf) and sd == od and sn < on))
def univariate_decomposition(self):
@@ -771,12 +771,13 @@ class FFPDElement(sage.structure.element.RingElement):
whole, p = p.quo_rem(q)
else:
whole = p
- p = R(1)
+ p = R.one()
df = self.denominator_factored()
decomp = [self.parent()(whole, [])]
for (a, m) in df:
- numer = p * prod([b**n for (b, n) in df if b != a]).\
- inverse_mod(a**m) % (a**m)
+ numer = p * prod([b ** n
+ for (b, n) in df
+ if b != a]).inverse_mod(a ** m) % (a ** m)
# The inverse exists because the product and a**m
# are relatively prime.
decomp.append(self.parent()(numer, [(a, m)]))
@@ -820,9 +821,9 @@ class FFPDElement(sage.structure.element.RingElement):
return None
df = self.denominator_factored()
- J = R.ideal([q**e for q, e in df])
- if R(1) in J:
- return R(1).lift(J)
+ J = R.ideal([q ** e for q, e in df])
+ if R.one() in J:
+ return R.one().lift(J)
return None
def nullstellensatz_decomposition(self):
@@ -895,8 +896,8 @@ class FFPDElement(sage.structure.element.RingElement):
p = self.numerator()
df = self.denominator_factored()
m = len(df)
- iteration1 = FFPDSum([self.parent()(p*L[i],[df[j]
- for j in xrange(m) if j != i])
+ iteration1 = FFPDSum([self.parent()(p * L[i],
+ [df[j] for j in xrange(m) if j != i])
for i in xrange(m) if L[i] != 0])
# Now decompose each FFPD of iteration1.
@@ -977,7 +978,7 @@ class FFPDElement(sage.structure.element.RingElement):
Xs = list(R.gens())
d = len(Xs)
- # Expand R by 2*m new variables.
+ # Expand R by 2 * m new variables.
S = 'S'
while S in [str(x) for x in Xs]:
S = S + 'S'
@@ -994,13 +995,13 @@ class FFPDElement(sage.structure.element.RingElement):
Ts = RR.gens()[d + m: d + 2 * m]
# Compute the appropriate elimination ideal.
- J = RR.ideal([ Ss[j] - RR(df[j][0]) for j in xrange(m)] +\
- [ Ss[j]**df[j][1] - Ts[j] for j in xrange(m)])
+ J = RR.ideal([Ss[j] - RR(df[j][0]) for j in xrange(m)] +
+ [Ss[j] ** df[j][1] - Ts[j] for j in xrange(m)])
J = J.elimination_ideal(Xs + Ss)
# Coerce J into the polynomial ring in the indeteminates Ts[m:].
# I choose the negdeglex order because i find it useful in my work.
- RRR = PolynomialRing(F, [str(t) for t in Ts], order ='negdeglex')
+ RRR = PolynomialRing(F, [str(t) for t in Ts], order='negdeglex')
return RRR.ideal(J)
def algebraic_dependence_decomposition(self, whole_and_parts=True):
@@ -1090,25 +1091,25 @@ class FFPDElement(sage.structure.element.RingElement):
# so that g.lt() is indeed the monomial we want below.
# Use g to rewrite r into a sum of FFPDs,
# each with < m distinct denominator factors.
- gg = (g.lt() - g)/(g.lc())
+ gg = (g.lt() - g) / (g.lc())
numers = map(prod, zip(gg.coefficients(), gg.monomials()))
- e = list(g.lt().exponents())[0:m]
+ e = list(g.lt().exponents())[0: m]
denoms = [(new_vars[j], e[0][j] + 1) for j in xrange(m)]
# Write r in terms of new_vars,
# cancel factors in the denominator, and combine like terms.
FFPD = FractionWithFactoredDenominatorRing(J.ring())
- iteration1_temp = FFPDSum([FFPD(a, denoms) for a in numers]).\
- combine_like_terms()
+ iteration1_temp = FFPDSum([FFPD(a, denoms)
+ for a in numers]).combine_like_terms()
# Substitute in df.
- qpowsub = dict([(new_vars[j], df[j][0]**df[j][1])
+ qpowsub = dict([(new_vars[j], df[j][0] ** df[j][1])
for j in xrange(m)])
iteration1 = FFPDSum()
for r in iteration1_temp:
- num1 = p*J.ring()(r.numerator()).subs(qpowsub)
+ num1 = p * J.ring()(r.numerator()).subs(qpowsub)
denoms1 = []
for q, e in r.denominator_factored():
j = new_vars.index(q)
- denoms1.append((df[j][0], df[j][1]*e))
+ denoms1.append((df[j][0], df[j][1] * e))
iteration1.append(self.parent()(num1, denoms1))
# Now decompose each FFPD of iteration1.
for r in iteration1:
@@ -1334,7 +1335,7 @@ class FFPDElement(sage.structure.element.RingElement):
new_df[i][1] -= 1
else:
del(new_df[i])
- iteration1.append(self.parent()(p*L[i], new_df))
+ iteration1.append(self.parent()(p * L[i], new_df))
# Contributions from dets.
# Compute each contribution's cohomologous form using
@@ -1353,13 +1354,13 @@ class FFPDElement(sage.structure.element.RingElement):
# Sort variables according to the term order of R.
x = sorted(var_sets_n[k])
# Compute Jacobian in the Symbolic Ring.
- jac = jacobian([SR(p*L[n + k])] +
+ jac = jacobian([SR(p * L[n + k])] +
[SR(qs[j]) for j in xrange(n) if j != J],
[SR(xx) for xx in x])
det = jac.determinant()
psign = FractionWithFactoredDenominatorRing._permutation_sign(x, X)
- iteration1.append(self.parent()((-1)**J*det/\
- (psign*new_df[J][1]),
+ iteration1.append(self.parent()((-1) ** J * det /
+ (psign * new_df[J][1]),
new_df))
# Now decompose each FFPD of iteration1.
@@ -1434,10 +1435,11 @@ class FFPDElement(sage.structure.element.RingElement):
# Cauchy differential form generated by each element of decomp.
if asy_var is None:
asy_var = var('r')
- cauchy_stuff = prod([X[j]**(-alpha[j]*asy_var - 1) for j in xrange(d)])
+ cauchy_stuff = prod([X[j] ** (-alpha[j] * asy_var - 1)
+ for j in xrange(d)])
decomp2 = FFPDSum()
for f in decomp1:
- ff = self.parent()(f.numerator()*cauchy_stuff,
+ ff = self.parent()(f.numerator() * cauchy_stuff,
f.denominator_factored())
decomp2.extend(ff.cohomology_decomposition())
decomp2 = decomp2.combine_like_terms()
@@ -1445,14 +1447,15 @@ class FFPDElement(sage.structure.element.RingElement):
# Divide out cauchy_stuff from integrands.
decomp3 = FFPDSum()
for f in decomp2:
- ff = self.parent()((f.numerator()/cauchy_stuff).\
- simplify_full().collect(asy_var),
+ ff = self.parent()((f.numerator() /
+ cauchy_stuff).simplify_full().collect(asy_var),
f.denominator_factored())
decomp3.append(ff)
return decomp3
- def asymptotics(self, p, alpha, N, asy_var=None, numerical=0, verbose=False):
+ def asymptotics(self, p, alpha, N, asy_var=None, numerical=0,
+ verbose=False):
r"""
Return the first `N` terms (some of which could be zero)
of the asymptotic expansion of the Maclaurin ray coefficients
@@ -1587,7 +1590,7 @@ class FFPDElement(sage.structure.element.RingElement):
# (X[i]*diff(h, X[i])).subs(p) != 0.
# Assuming such an i exists.
i = d - 1
- while 0 in [(X[i]*diff(h, X[i])).subs(p) for (h, e) in df]:
+ while 0 in [(X[i] * diff(h, X[i])).subs(p) for (h, e) in df]:
i -= 1
coordinate = i
@@ -1714,22 +1717,22 @@ class FFPDElement(sage.structure.element.RingElement):
# I'll just past the code from the multiple point case.
if d == 1:
det = jacobian(H, X).subs(p).determinant().abs()
- exp_scale = prod([(p[X[i]]**(-alpha[i])).subs(p)
- for i in xrange(d)] )
- subexp_part = -G.subs(p)/(det*prod(p.values()))
+ exp_scale = prod([(p[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)])
+ subexp_part = -G.subs(p) / (det * prod(p.values()))
if numerical:
exp_scale = exp_scale.n(digits=numerical)
subexp_part = subexp_part.n(digits=numerical)
- return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
+ return (exp_scale ** asy_var * subexp_part, exp_scale, subexp_part)
# If p is a tuple of rationals, then compute with it directly.
# Otherwise, compute symbolically and plug in p at the end.
- if vector(p.values()) in QQ**d:
+ if vector(p.values()) in QQ ** d:
P = p
else:
sP = [var('p' + str(j)) for j in xrange(d)]
- P = dict( [(X[j], sP[j]) for j in xrange(d)] )
- p = dict( [(sP[j], p[X[j]]) for j in xrange(d)] )
+ P = dict([(X[j], sP[j]) for j in xrange(d)])
+ p = dict([(sP[j], p[X[j]]) for j in xrange(d)])
# Setup.
if verbose:
@@ -1739,83 +1742,84 @@ class FFPDElement(sage.structure.element.RingElement):
U = function('U', *tuple(X))
# All other functions are defined in terms of h, U, and
# explicit functions.
- Gcheck = -G/U * (h/X[d - 1])
- A = Gcheck.subs({X[d - 1]: Integer(1)/h})/h
+ Gcheck = -G / U * (h / X[d - 1])
+ A = Gcheck.subs({X[d - 1]: Integer(1) / h}) / h
t = 't'
while t in [str(x) for x in X]:
t = t + 't'
T = [var(t + str(i)) for i in xrange(d - 1)]
- e = dict([(X[i], P[X[i]]*exp(I*T[i])) for i in xrange(d - 1)])
+ e = dict([(X[i], P[X[i]] * exp(I * T[i])) for i in xrange(d - 1)])
ht = h.subs(e)
At = A.subs(e)
- Phit = -log(P[X[d - 1]]*ht) + \
- I * add([alpha[i]/alpha[d - 1]*T[i] for i in xrange(d - 1)])
- Tstar = dict([(t, Integer(0)) for t in T])
+ Phit = (-log(P[X[d - 1]] * ht) +
+ I * sum([alpha[i] / alpha[d - 1] * T[i]
+ for i in xrange(d - 1)]))
+ Tstar = {t: Integer(0) for t in T}
# Store h and U and all their derivatives evaluated at P.
atP = P.copy()
- atP.update({h.subs(P): Integer(1)/P[X[d - 1]]})
+ atP.update({h.subs(P): Integer(1) / P[X[d - 1]]})
- # Compute the derivatives of h up to order 2*N, evaluate at P,
+ # Compute the derivatives of h up to order 2 * N, evaluate at P,
# and store in atP.
# Keep a copy of unevaluated h derivatives for use in the case
# d = 2 and v > 2 below.
hderivs1 = {} # First derivatives of h.
for i in xrange(d - 1):
- s = solve( diff(H.subs({X[d - 1]: Integer(1)/h}), X[i]),
+ s = solve(diff(H.subs({X[d - 1]: Integer(1) / h}), X[i]),
diff(h, X[i]))[0].rhs().simplify()
hderivs1.update({diff(h, X[i]): s})
atP.update({diff(h, X[i]).subs(P): s.subs(P).subs(atP)})
hderivs = FractionWithFactoredDenominatorRing._diff_all(
- h, X[0: d - 1], 2*N, sub=hderivs1, rekey=h)
+ h, X[0: d - 1], 2 * N, sub=hderivs1, rekey=h)
for k in hderivs.keys():
- atP.update({k.subs(P):hderivs[k].subs(atP)})
+ atP.update({k.subs(P): hderivs[k].subs(atP)})
- # Compute the derivatives of U up to order 2*N and evaluate at P.
+ # Compute the derivatives of U up to order 2 * N and evaluate at P.
# To do this, differentiate H = U*Hcheck over and over, evaluate at P,
# and solve for the derivatives of U at P.
# Need the derivatives of H with short keys to pass on
# to diff_prod later.
Hderivs = FractionWithFactoredDenominatorRing._diff_all(
- H, X, 2*N, ending=[X[d - 1]], sub_final=P)
+ H, X, 2 * N, ending=[X[d - 1]], sub_final=P)
if verbose:
print("Computing derivatives of auxiallary functions...")
# For convenience in checking if all the nontrivial derivatives of U
# at p are zero a few line below, store the value of U(p) in atP
# instead of in Uderivs.
- Uderivs ={}
+ Uderivs = {}
atP.update({U.subs(P): diff(H, X[d - 1]).subs(P)})
end = [X[d - 1]]
- Hcheck = X[d - 1] - Integer(1)/h
+ Hcheck = X[d - 1] - Integer(1) / h
k = H.polynomial(CC).degree() - 1
if k == 0:
# Then we can conclude that all higher derivatives of U are zero.
- for l in xrange(1, 2*N + 1):
+ for l in xrange(1, 2 * N + 1):
for s in UnorderedTuples(X, l):
Uderivs[diff(U, s).subs(P)] = Integer(0)
- elif k > 0 and k < 2*N:
+ elif k > 0 and k < 2 * N:
all_zero = True
- Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
+ Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
Hderivs, U, Hcheck, X,
range(1, k + 1), end, Uderivs, atP)
# Check for a nonzero U derivative.
- if Uderivs.values() != [Integer(0) for i in xrange(len(Uderivs))]:
+ if any(u for u in Uderivs.values()):
all_zero = False
if all_zero:
# Then, using a proposition at the end of [RaWi2012], we can
# conclude that all higher derivatives of U are zero.
- for l in xrange(k + 1, 2*N +1):
+ for l in xrange(k + 1, 2 * N + 1):
for s in UnorderedTuples(X, l):
Uderivs.update({diff(U, s).subs(P): Integer(0)})
else:
# Have to compute the rest of the derivatives.
Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
Hderivs, U, Hcheck, X,
- range(k + 1, 2*N + 1), end, Uderivs,
+ range(k + 1, 2 * N + 1), end, Uderivs,
atP)
else:
Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
Hderivs, U, Hcheck, X,
- range(1, 2*N + 1), end, Uderivs, atP)
+ range(1, 2 * N + 1), end, Uderivs, atP)
atP.update(Uderivs)
# In general, this algorithm is not designed to handle the case of a
@@ -1829,7 +1833,7 @@ class FFPDElement(sage.structure.element.RingElement):
splat = Phitderiv.subs(Tstar).subs(atP).subs(p).simplify()
while splat == 0:
v += 1
- if v > 2*N:
+ if v > 2 * N:
# Then need to compute more derivatives of h for atP.
hderivs.update({diff(h, X[0], v):
diff(hderivs[diff(h, X[0], v - 1)],
@@ -1841,8 +1845,8 @@ class FFPDElement(sage.structure.element.RingElement):
if d == 2 and v > 2:
t = T[0] # Simplify variable names.
- a = splat/factorial(v)
- Phitu = Phit - a*t**v
+ a = splat / factorial(v)
+ Phitu = Phit - a * t ** v
# Compute all partial derivatives of At and Phitu
# up to orders 2*(N - 1) and 2*(N - 1) + v, respectively,
@@ -1858,11 +1862,11 @@ class FFPDElement(sage.structure.element.RingElement):
BB = function('BB', t)
if v.mod(2) == 0:
At_derivs = FractionWithFactoredDenominatorRing._diff_all(
- At, T, 2*N - 2,
+ At, T, 2 * N - 2,
sub=hderivs1, sub_final=[Tstar, atP],
rekey=AA)
Phitu_derivs = FractionWithFactoredDenominatorRing._diff_all(
- Phitu, T, 2*N - 2 +v,
+ Phitu, T, 2 * N - 2 +v,
sub=hderivs1,
sub_final=[Tstar, atP],
zero_order=v + 1, rekey=BB)
@@ -1888,25 +1892,26 @@ class FFPDElement(sage.structure.element.RingElement):
L = []
if v.mod(2) == 0:
for k in xrange(N):
- L.append(add([
- (-1)**l * gamma((2*k + v*l + 1)/v)/\
- (factorial(l) * factorial(2*k + v*l))*\
- DD[(k, l)] for l in xrange(0, 2*k + 1) ]))
- chunk = a**(-1/v)/(pi*v)*add([
- alpha[d - 1 ]**(-(2*k + 1)/v)*\
- L[k]*asy_var**(-(2*k + 1)/v) for k in xrange(N) ])
+ L.append(sum([(-1) ** l * gamma((2 * k + v * l + 1) / v) /
+ (factorial(l) * factorial(2 * k + v * l)) *
+ DD[(k, l)] for l in xrange(0, 2 * k + 1)]))
+ chunk = (a ** (-1 / v) / (pi * v) *
+ sum([alpha[d - 1] ** (-(2 * k + 1) / v) *
+ L[k] * asy_var ** (-(2 * k + 1) / v)
+ for k in xrange(N)]))
else:
- zeta = exp(I*pi/(2*v))
+ zeta = exp(I * pi / (2 * v))
for k in xrange(N):
- L.append(add([
- (-1)**l*gamma((k + v*l + 1)/v)/\
- (factorial(l)*factorial(k + v*l))*\
- (zeta**(k + v*l + 1) +\
- (-1)**(k + v*l)*zeta**(-(k + v*l + 1)))*\
- DD[(k, l)] for l in xrange(0, k + 1) ]))
- chunk = abs(a)**(-1/v)/(2*pi*v)*add([
- alpha[d - 1]**(-(k + 1)/v)*\
- L[k] *asy_var**(-(k + 1)/v) for k in xrange(N) ])
+ L.append(sum([(-1) ** l * gamma((k + v * l + 1) / v) /
+ (factorial(l) * factorial(k + v * l)) *
+ (zeta ** (k + v * l + 1) +
+ (-1) ** (k + v * l) *
+ zeta ** (-(k + v * l + 1))) *
+ DD[(k, l)] for l in xrange(0, k + 1)]))
+ chunk = (abs(a) ** (-1 / v) / (2 * pi * v) *
+ sum([alpha[d - 1] ** (-(k + 1) / v) *
+ L[k] * asy_var ** (-(k + 1) / v)
+ for k in xrange(N)]))
# Asymptotics for d >= 2 case.
# A singular Phit''(Tstar) will cause a crash in this case.
@@ -1914,11 +1919,11 @@ class FFPDElement(sage.structure.element.RingElement):
Phit1 = jacobian(Phit, T).subs(hderivs1)
a = jacobian(Phit1, T).subs(hderivs1).subs(Tstar).subs(atP)
a_inv = a.inverse()
- Phitu = Phit - (Integer(1)/Integer(2))*matrix([T])*\
- a*matrix([T]).transpose()
+ Phitu = (Phit - (1 / QQ(2)) * matrix([T]) *
+ a * matrix([T]).transpose())
Phitu = Phitu[0][0]
# Compute all partial derivatives of At and Phitu up to
- # orders 2*N-2 and 2*N, respectively.
+ # orders 2 * N-2 and 2 * N, respectively.
# Take advantage of the fact that At and Phitu
# are sufficiently differentiable functions so that mixed partials
# are equal. Thus only need to compute representative partials.
@@ -1931,11 +1936,11 @@ class FFPDElement(sage.structure.element.RingElement):
print("Computing derivatives of more auxiliary functions...")
AA = function('AA', *tuple(T))
At_derivs = FractionWithFactoredDenominatorRing._diff_all(
- At, T, 2*N - 2, sub=hderivs1,
+ At, T, 2 * N - 2, sub=hderivs1,
sub_final =[Tstar, atP], rekey=AA)
BB = function('BB', *tuple(T))
Phitu_derivs = FractionWithFactoredDenominatorRing._diff_all(
- Phitu, T, 2*N, sub=hderivs1,
+ Phitu, T, 2 * N, sub=hderivs1,
sub_final =[Tstar, atP], rekey=BB,
zero_order=3)
AABB_derivs = At_derivs
@@ -1948,32 +1953,33 @@ class FFPDElement(sage.structure.element.RingElement):
AA, BB, AABB_derivs, T, a_inv, 1 , N)
# Plug above into asymptotic formula.
- L =[]
+ L = []
for k in xrange(N):
- L.append(add([
- DD[(0, k, l)]/((-1)**k*2**(l+k)*\
- factorial(l)*factorial(l+k))
- for l in xrange(0, 2*k + 1) ]))
- chunk = add([ (2*pi)**((1 - d)/Integer(2))*\
- a.determinant()**(-Integer(1)/Integer(2))*\
- alpha[d - 1]**((Integer(1) - d)/Integer(2) - k)*L[k]*\
- asy_var**((Integer(1) - d)/Integer(2) - k)
- for k in xrange(N) ])
+ L.append(sum([DD[(0, k, l)] / ((-1) ** k * 2 ** (l + k) *
+ factorial(l) * factorial(l + k))
+ for l in xrange(0, 2 * k + 1)]))
+ chunk = sum([(2 * pi) ** ((1 - d) / Integer(2)) *
+ a.determinant() ** (-Integer(1) / Integer(2)) *
+ alpha[d - 1] ** ((Integer(1) - d) / Integer(2) - k) *
+ L[k] *
+ asy_var ** ((Integer(1) - d) / Integer(2) - k)
+ for k in xrange(N)])
chunk = chunk.subs(p).simplify()
coeffs = chunk.coefficients(asy_var)
coeffs.reverse()
coeffs = coeffs[:N]
if numerical:
- subexp_part = add([co[0].subs(p).n(digits=numerical)*\
- asy_var**co[1] for co in coeffs])
- exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
+ subexp_part = sum([co[0].subs(p).n(digits=numerical) *
+ asy_var ** co[1] for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
for i in xrange(d)]).n(digits=numerical)
else:
- subexp_part = add([co[0].subs(p)*asy_var**co[1] for co in coeffs])
- exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
- for i in xrange(d)])
- return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
+ subexp_part = sum([co[0].subs(p) * asy_var ** co[1]
+ for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)])
+ return (exp_scale ** asy_var * subexp_part, exp_scale, subexp_part)
def asymptotics_multiple(self, p, alpha, N, asy_var, coordinate=None,
numerical=0, verbose=False):
@@ -2081,7 +2087,7 @@ class FFPDElement(sage.structure.element.RingElement):
H = [SR(h) for (h, e) in self.denominator_factored()]
Hprod = prod(H)
n = len(H)
- P = dict([(SR(x), p[x]) for x in R.gens()])
+ P = {SR(x): p[x] for x in R.gens()}
Sstar = self._crit_cone_combo(p, alpha, coordinate)
# Put the given convenient variable at end of variable list.
@@ -2091,25 +2097,24 @@ class FFPDElement(sage.structure.element.RingElement):
a = alpha.pop(coordinate)
alpha.append(a)
-
# Case n = d.
if n == d:
det = jacobian(H, X).subs(P).determinant().abs()
- exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(P)
- for i in xrange(d)] )
- subexp_part = G.subs(P)/(det*prod(P.values()))
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(P)
+ for i in xrange(d)])
+ subexp_part = G.subs(P) / (det * prod(P.values()))
if numerical:
exp_scale = exp_scale.n(digits=numerical)
subexp_part = subexp_part.n(digits=numerical)
- return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
+ return (exp_scale ** asy_var * subexp_part, exp_scale, subexp_part)
# Case n < d.
# If P is a tuple of rationals, then compute with it directly.
# Otherwise, compute symbolically and plug in P at the end.
- if vector(P.values()) not in QQ**d:
+ if vector(P.values()) not in QQ ** d:
sP = [var('p' + str(j)) for j in xrange(d)]
- P = dict( [(X[j], sP[j]) for j in xrange(d)] )
- p = dict( [(sP[j], p[X[j]]) for j in xrange(d)] )
+ P = dict([(X[j], sP[j]) for j in xrange(d)])
+ p = dict([(sP[j], p[X[j]]) for j in xrange(d)])
# Setup.
if verbose:
@@ -2131,35 +2136,36 @@ class FFPDElement(sage.structure.element.RingElement):
U = function('U', *tuple(X))
# All other functions are defined in terms of h, U, and
# explicit functions.
- Hcheck = prod([X[d - 1] - Integer(1)/h[j] for j in xrange(n)])
- Gcheck = -G/U * prod([-h[j]/X[d - 1] for j in xrange(n)])
- A = [(-1)**(n - 1)*X[d - 1]**(-n + j)*\
- diff(Gcheck.subs({X[d - 1]: Integer(1)/X[d - 1]}), X[d - 1], j)
+ Hcheck = prod([X[d - 1] - Integer(1) / h[j] for j in xrange(n)])
+ Gcheck = -G / U * prod([-h[j] / X[d - 1] for j in xrange(n)])
+ A = [(-1) ** (n - 1) * X[d - 1] ** (-n + j) *
+ diff(Gcheck.subs({X[d - 1]: Integer(1) / X[d - 1]}), X[d - 1], j)
for j in xrange(n)]
- e = dict([(X[i], P[X[i]]*exp(I*T[i])) for i in xrange(d - 1)])
+ e = dict([(X[i], P[X[i]] * exp(I * T[i])) for i in xrange(d - 1)])
ht = [hh.subs(e) for hh in h]
- hsumt = add([S[j]*ht[j] for j in xrange(n - 1)]) +\
- (Integer(1) - add(S))*ht[n - 1]
+ hsumt = (sum([S[j] * ht[j] for j in xrange(n - 1)]) +
+ (Integer(1) - sum(S)) * ht[n - 1])
At = [AA.subs(e).subs({X[d - 1]: hsumt}) for AA in A]
- Phit = -log(P[X[d - 1]]*hsumt) +\
- I*add([alpha[i]/alpha[d - 1]*T[i] for i in xrange(d - 1)])
+ Phit = (-log(P[X[d - 1]] * hsumt) +
+ I * sum([alpha[i] / alpha[d - 1] * T[i]
+ for i in xrange(d - 1)]))
# atP Stores h and U and all their derivatives evaluated at C.
atP = P.copy()
- atP.update(dict([(hh.subs(P), Integer(1)/P[X[d - 1]]) for hh in h]))
+ atP.update(dict([(hh.subs(P), Integer(1) / P[X[d - 1]]) for hh in h]))
- # Compute the derivatives of h up to order 2*N and evaluate at P.
+ # Compute the derivatives of h up to order 2 * N and evaluate at P.
hderivs1 = {} # First derivatives of h.
for (i, j) in mrange([d - 1, n]):
- s = solve(diff(H[j].subs({X[d - 1]: Integer(1)/h[j]}), X[i]),
+ s = solve(diff(H[j].subs({X[d - 1]: Integer(1) / h[j]}), X[i]),
diff(h[j], X[i]))[0].rhs().simplify()
hderivs1.update({diff(h[j], X[i]): s})
atP.update({diff(h[j], X[i]).subs(P): s.subs(P).subs(atP)})
hderivs = FractionWithFactoredDenominatorRing._diff_all(
- h, X[0:d - 1], 2*N, sub=hderivs1, rekey=h)
+ h, X[0:d - 1], 2 * N, sub=hderivs1, rekey=h)
for k in hderivs.keys():
atP.update({k.subs(P): hderivs[k].subs(atP)})
- # Compute the derivatives of U up to order 2*N - 2 + min{n, N} - 1 and
+ # Compute the derivatives of U up to order 2 * N - 2 + min{n, N} - 1 and
# evaluate at P.
# To do this, differentiate H = U*Hcheck over and over, evaluate at P,
# and solve for the derivatives of U at P.
@@ -2170,49 +2176,49 @@ class FFPDElement(sage.structure.element.RingElement):
m = min(n, N)
end = [X[d-1] for j in xrange(n)]
Hprodderivs = FractionWithFactoredDenominatorRing._diff_all(
- Hprod, X, 2*N - 2 + n, ending=end,
+ Hprod, X, 2 * N - 2 + n, ending=end,
sub_final=P)
atP.update({U.subs(P): diff(Hprod, X[d - 1], n).subs(P)/factorial(n)})
- Uderivs ={}
+ Uderivs = {}
k = Hprod.polynomial(CC).degree() - n
if k == 0:
# Then we can conclude that all higher derivatives of U are zero.
- for l in xrange(1, 2*N - 2 + m):
+ for l in xrange(1, 2 * N - 2 + m):
for s in UnorderedTuples(X, l):
Uderivs[diff(U, s).subs(P)] = Integer(0)
- elif k > 0 and k < 2*N - 2 + m - 1:
+ elif k > 0 and k < 2 * N - 2 + m - 1:
all_zero = True
Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
Hprodderivs, U, Hcheck, X,
range(1, k + 1), end, Uderivs, atP)
# Check for a nonzero U derivative.
- if Uderivs.values() != [Integer(0) for i in xrange(len(Uderivs))]:
+ if any(u for u in Uderivs.values()):
all_zero = False
if all_zero:
# Then all higher derivatives of U are zero.
- for l in xrange(k + 1, 2*N - 2 + m):
+ for l in xrange(k + 1, 2 * N - 2 + m):
for s in UnorderedTuples(X, l):
Uderivs.update({diff(U, s).subs(P): Integer(0)})
else:
# Have to compute the rest of the derivatives.
Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
Hprodderivs, U, Hcheck, X,
- range(k + 1, 2*N - 2 + m), end,
+ range(k + 1, 2 * N - 2 + m), end,
Uderivs, atP)
else:
Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
Hprodderivs, U, Hcheck, X,
- range(1, 2*N - 2 + m), end, Uderivs, atP)
+ range(1, 2 * N - 2 + m), end, Uderivs, atP)
atP.update(Uderivs)
Phit1 = jacobian(Phit, T + S).subs(hderivs1)
a = jacobian(Phit1, T + S).subs(hderivs1).subs(thetastar).subs(atP)
a_inv = a.inverse()
- Phitu = Phit - (Integer(1)/Integer(2))*matrix([T + S])*a*\
- matrix([T + S]).transpose()
+ Phitu = (Phit - (1 / Integer(2)) * matrix([T + S]) * a *
+ matrix([T + S]).transpose())
Phitu = Phitu[0][0]
- # Compute all partial derivatives of At and Phitu up to orders 2*N - 2
- # and 2*N, respectively. Take advantage of the fact that At and Phitu
+ # Compute all partial derivatives of At and Phitu up to orders 2 * N - 2
+ # and 2 * N, respectively. Take advantage of the fact that At and Phitu
# are sufficiently differentiable functions so that mixed partials
# are equal. Thus only need to compute representative partials.
# Choose nondecreasing sequences as representative differentiation-
@@ -2223,12 +2229,12 @@ class FFPDElement(sage.structure.element.RingElement):
print("Computing derivatives of more auxiliary functions...")
AA = [function('A' + str(j), *tuple(T + S)) for j in xrange(n)]
At_derivs = FractionWithFactoredDenominatorRing._diff_all(
- At, T + S, 2*N - 2, sub=hderivs1,
- sub_final =[thetastar, atP], rekey=AA)
+ At, T + S, 2 * N - 2, sub=hderivs1,
+ sub_final=[thetastar, atP], rekey=AA)
BB = function('BB', *tuple(T + S))
Phitu_derivs = FractionWithFactoredDenominatorRing._diff_all(
- Phitu, T + S, 2*N, sub=hderivs1,
- sub_final =[thetastar, atP], rekey=BB,
+ Phitu, T + S, 2 * N, sub=hderivs1,
+ sub_final=[thetastar, atP], rekey=BB,
zero_order=3)
AABB_derivs = At_derivs
AABB_derivs.update(Phitu_derivs)
@@ -2243,33 +2249,38 @@ class FFPDElement(sage.structure.element.RingElement):
L = {}
for (j, k) in product(xrange(min(n, N)), xrange(max(0, N - 1 - n), N)):
if j + k <= N - 1:
- L[(j, k)] = add([DD[(j, k, l)]/((-1)**k*2**(k + l)*\
- factorial(l)*factorial(k + l))
- for l in xrange(2*k + 1)])
- det = a.determinant()**(-Integer(1)/Integer(2))*\
- (2*pi)**((n - d)/Integer(2))
- chunk = det*add([
- (alpha[d - 1]*asy_var)**((n - d)/Integer(2) - q)*\
- add([L[(j, k)]*binomial(n - 1, j)*\
- stirling_number1(n - j, n + k - q)*(-1)**(q - j - k)
- for (j, k) in product(xrange(min(n - 1, q) + 1),
- xrange(max(0, q - n), q + 1))
- if j + k <= q])
- for q in xrange(N)])
+ L[(j, k)] = sum([DD[(j, k, l)] / ((-1) ** k * 2 ** (k + l) *
+ factorial(l) *
+ factorial(k + l))
+ for l in xrange(2 * k + 1)])
+ det = (a.determinant() ** (-1 / Integer(2)) *
+ (2 * pi) ** ((n - d) / Integer(2)))
+ chunk = det * sum([(alpha[d - 1] * asy_var) ** ((n - d) /
+ Integer(2) - q) *
+ sum([L[(j, k)] * binomial(n - 1, j) *
+ stirling_number1(n - j, n + k - q) *
+ (-1) ** (q - j - k)
+ for (j, k) in product(xrange(min(n - 1, q) + 1),
+ xrange(max(0, q - n),
+ q + 1))
+ if j + k <= q])
+ for q in xrange(N)])
chunk = chunk.subs(P).simplify()
coeffs = chunk.coefficients(asy_var)
coeffs.reverse()
coeffs = coeffs[:N]
if numerical:
- subexp_part = add([co[0].subs(p).n(digits=numerical)*asy_var**co[1]
+ subexp_part = sum([co[0].subs(p).n(digits=numerical) *
+ asy_var ** co[1]
for co in coeffs])
- exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
for i in xrange(d)]).n(digits=numerical)
else:
- subexp_part = add([co[0].subs(p)*asy_var**co[1] for co in coeffs])
- exp_scale = prod([(P[X[i]]**(-alpha[i])).subs(p)
+ subexp_part = sum([co[0].subs(p) * asy_var ** co[1]
+ for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
for i in xrange(d)])
- return (exp_scale**asy_var*subexp_part, exp_scale, subexp_part)
+ return (exp_scale ** asy_var * subexp_part, exp_scale, subexp_part)
def _crit_cone_combo(self, p, alpha, coordinate=None):
@@ -2331,8 +2342,8 @@ class FFPDElement(sage.structure.element.RingElement):
# solve_left() fails when working in SR :-(.
# So use solve() instead.
# Gamma.solve_left(vector(beta))
- V = [var('sss'+str(i)) for i in range(n)]
- M = matrix(V)*Gamma
+ V = [var('sss' + str(i)) for i in range(n)]
+ M = matrix(V) * Gamma
eqns = [M[0][j] == beta[j] for j in range(d)]
s = solve(eqns, V, solution_dict=True)[0] # Assume a unique solution.
return [s[v] for v in V]
@@ -2354,7 +2365,7 @@ class FFPDElement(sage.structure.element.RingElement):
sage: R.<x,y> = PolynomialRing(QQ)
sage: FFPD = FractionWithFactoredDenominatorRing(R)
sage: p = exp(x)
- sage: df = [(x**3 + 3*y^2, 5), (x*y, 2), (y, 1)]
+ sage: df = [(x^3 + 3*y^2, 5), (x*y, 2), (y, 1)]
sage: f = FFPD(p, df)
sage: f
(e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)])
@@ -2401,7 +2412,7 @@ class FFPDElement(sage.structure.element.RingElement):
sage: R.<x,y> = PolynomialRing(QQ)
sage: FFPD = FractionWithFactoredDenominatorRing(R)
sage: p = exp(x)
- sage: df = [(x**3 + 3*y^2, 5), (x*y, 2), (y, 1)]
+ sage: df = [(x^3 + 3*y^2, 5), (x*y, 2), (y, 1)]
sage: f = FFPD(p, df)
sage: f
(e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)])
@@ -2426,7 +2437,7 @@ class FFPDElement(sage.structure.element.RingElement):
d = self.dimension()
H = [h for (h, e) in self.denominator_factored()]
n = len(H)
- return [tuple([(X[j]*diff(H[i], X[j])).subs(p) for j in xrange(d)])
+ return [tuple([(X[j] * diff(H[i], X[j])).subs(p) for j in xrange(d)])
for i in xrange(n)]
def critical_cone(self, p, coordinate=None):
@@ -2451,10 +2462,10 @@ class FFPDElement(sage.structure.element.RingElement):
EXAMPLES::
sage: from sage.combinat.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
- sage: R.<x,y,z>= PolynomialRing(QQ)
+ sage: R.<x,y,z> = PolynomialRing(QQ)
sage: FFPD = FractionWithFactoredDenominatorRing(R)
sage: G = 1
- sage: H = (1 - x*(1 + y))*(1 - z*x**2*(1 + 2*y))
+ sage: H = (1 - x*(1 + y)) * (1 - z*x**2*(1 + 2*y))
sage: Hfac = H.factor()
sage: G = 1/Hfac.unit()
sage: F = FFPD(G, Hfac)
@@ -2512,11 +2523,11 @@ class FFPDElement(sage.structure.element.RingElement):
EXAMPLES::
sage: from sage.combinat.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
- sage: R.<x,y,z>= PolynomialRing(QQ)
+ sage: R.<x,y,z> = PolynomialRing(QQ)
sage: FFPD = FractionWithFactoredDenominatorRing(R)
- sage: H = (1 - x*(1 + y))*(1 - z*x**2*(1 + 2*y))
+ sage: H = (1 - x*(1 + y)) * (1 - z*x**2*(1 + 2*y))
sage: df = H.factor()
- sage: G = 1/df.unit()
+ sage: G = 1 / df.unit()
sage: F = FFPD(G, df)
sage: p1 = {x: 1/2, y: 1, z: 4/3}
sage: p2 = {x: 1, y: 2, z: 1/2}
@@ -2570,8 +2581,8 @@ class FFPDElement(sage.structure.element.RingElement):
# Tests all passed
X = R.gens()
- return ( True, 'convenient in variables {}'.format(
- [X[i] for i in convenient_coordinates]) )
+ convenientX = [X[i] for i in convenient_coordinates]
+ return (True, 'convenient in variables {}'.format(convenientX))
def singular_ideal(self):
r"""
@@ -2588,11 +2599,11 @@ class FFPDElement(sage.structure.element.RingElement):
EXAMPLES::
sage: from sage.combinat.asymptotics_multivariate_generating_functions import *
- sage: R.<x,y,z>= PolynomialRing(QQ)
+ sage: R.<x,y,z> = PolynomialRing(QQ)
sage: FFPD = FractionWithFactoredDenominatorRing(R)
- sage: H = (1 - x*(1 + y))**3*(1 - z*x**2*(1 + 2*y))
+ sage: H = (1 - x*(1 + y))^3 * (1 - z*x**2*(1 + 2*y))
sage: df = H.factor()
- sage: G = 1/df.unit()
+ sage: G = 1 / df.unit()
sage: F = FFPD(G, df)
sage: F.singular_ideal()
Ideal (x*y + x - 1, y^2 - 2*y*z + 2*y - z + 1, x*z + y - 2*z + 1)
@@ -2628,7 +2639,7 @@ class FFPDElement(sage.structure.element.RingElement):
sage: from sage.combinat.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
sage: R.<x,y> = PolynomialRing(QQ)
sage: FFPD = FractionWithFactoredDenominatorRing(R)
- sage: H = (1-x-y-x*y)^2
+ sage: H = (1 - x - y - x*y)^2
sage: Hfac = H.factor()
sage: G = 1/Hfac.unit()
sage: F = FFPD(G, Hfac)
@@ -2673,10 +2684,10 @@ class FFPDElement(sage.structure.element.RingElement):
# Find smooth, critical points for alpha.
X = S.gens()
Hred = S(Hred)
- J = S.ideal([Hred] +\
- [alpha[d - 1]*X[i]*diff(Hred, X[i]) -\
- alpha[i]*X[d - 1]*diff(Hred, X[d - 1])
- for i in xrange(d - 1)])
+ J = S.ideal([Hred] +
+ [alpha[d - 1] * X[i] * diff(Hred, X[i]) -
+ alpha[i] * X[d - 1] * diff(Hred, X[d - 1])
+ for i in xrange(d - 1)])
return S.ideal(J.groebner_basis())
def maclaurin_coefficients(self, multi_indices, numerical=0):
@@ -2708,7 +2719,7 @@ class FFPDElement(sage.structure.element.RingElement):
sage: FFPD = FractionWithFactoredDenominatorRing(R)
sage: H = 2 - 3*x
sage: Hfac = H.factor()
- sage: G = 1/Hfac.unit()
+ sage: G = 1 / Hfac.unit()
sage: F = FFPD(G, Hfac)
sage: F
(-1/3, [(x - 2/3, 1)])
@@ -2721,7 +2732,7 @@ class FFPDElement(sage.structure.element.RingElement):
sage: FFPD = FractionWithFactoredDenominatorRing(R)
sage: H = (4 - 2*x - y - z) * (4 - x - 2*y - z)
sage: Hfac = H.factor()
- sage: G = 16/Hfac.unit()
+ sage: G = 16 / Hfac.unit()
sage: F = FFPD(G, Hfac)
sage: alpha = vector([3, 3, 2])
sage: interval = [1, 2, 4]
@@ -2755,7 +2766,7 @@ class FFPDElement(sage.structure.element.RingElement):
# Create biggest multi-index needed.
alpha = []
for i in xrange(d):
- alpha.append(max((nu[i] for nu in multi_indices)))
+ alpha.append(max((nu[i] for nu in multi_indices)))
# Compute Maclaurin expansion of self up to index alpha.
# Use iterated univariate expansions.
@@ -2763,13 +2774,13 @@ class FFPDElement(sage.structure.element.RingElement):
f = SR(self.quotient())
X = [SR(x) for x in R.gens()]
for i in xrange(d):
- f = f.taylor(X[i], 0, alpha[i])
+ f = f.taylor(X[i], 0, alpha[i])
F = R(f)
# Collect coefficients.
X = R.gens()
for nu in multi_indices:
- monomial = prod([X[i]**nu[i] for i in xrange(d)])
+ monomial = prod([X[i] ** nu[i] for i in xrange(d)])
val = F.monomial_coefficient(monomial)
if numerical:
val = val.n(digits=numerical)
@@ -2813,7 +2824,7 @@ class FFPDElement(sage.structure.element.RingElement):
sage: FFPD = FractionWithFactoredDenominatorRing(R)
sage: H = 1 - x - y - x*y
sage: Hfac = H.factor()
- sage: G = 1/Hfac.unit()
+ sage: G = 1 / Hfac.unit()
sage: F = FFPD(G, Hfac)
sage: alpha = [1, 1]
sage: r = var('r')
@@ -2836,22 +2847,24 @@ class FFPDElement(sage.structure.element.RingElement):
#print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..."
# Get Maclaurin coefficients of self.
- multi_indices = [r*vector(alpha) for r in interval]
+ alpha = vector(alpha)
+ multi_indices = [r * alpha for r in interval]
mac = self.maclaurin_coefficients(multi_indices, numerical=digits)
#mac = self.old_maclaurin_coefficients(alpha, max(interval))
mac_approx = {}
stats = []
for r in interval:
- beta = tuple(r*vector(alpha))
- mac[beta] = (mac[beta]/exp_scale**r).n(digits=digits)
- mac_approx[beta] = [(f.subs({av:r})/exp_scale**r).n(digits=digits)
+ exp_s_r = exp_scale ** r
+ beta = tuple(r * alpha)
+ mac[beta] = (mac[beta] / exp_s_r).n(digits=digits)
+ mac_approx[beta] = [(f.subs({av: r}) / exp_s_r).n(digits=digits)
for f in approx]
stats_row = [beta, mac[beta], mac_approx[beta]]
if mac[beta] == 0:
- stats_row.extend([None for a in mac_approx[beta]])
+ stats_row.extend([None for a in mac_approx[beta]])
else:
- stats_row.append([(mac[beta] - a)/mac[beta]
- for a in mac_approx[beta]])
+ stats_row.append([(mac[beta] - a) / mac[beta]
+ for a in mac_approx[beta]])
stats.append(tuple(stats_row))
return stats
@@ -3303,16 +3316,16 @@ class FractionWithFactoredDenominatorRing(
EXAMPLES::
sage: from sage.combinat.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
- sage: u = ['a','b','c','d','e']
- sage: s = ['b','d']
+ sage: u = ['a', 'b', 'c', 'd', 'e']
+ sage: s = ['b', 'd']
sage: FFDR._permutation_sign(s, u)
-1
- sage: s = ['d','b']
+ sage: s = ['d', 'b']
sage: FFDR._permutation_sign(s, u)
1
"""
# Convert lists to lists of numbers in {1,..., len(u)}
- A = [i + 1 for i in xrange(len(u))]
+ A = [i + 1 for i in xrange(len(u))]
B = [u.index(x) + 1 for x in s]
C = sorted(list(Set(A).difference(Set(B))))
@@ -3376,18 +3389,18 @@ class FractionWithFactoredDenominatorRing(
for ff in f:
for D in sub:
if isinstance(ff, dict):
- ff = dict( [(k, ff[k].subs(D)) for k in ff.keys()] )
+ ff = dict([(k, ff[k].subs(D)) for k in ff.keys()])
else:
ff = ff.subs(D)
g.append(ff)
if singleton and simplify:
- if isinstance(g[Integer(0) ], dict):
- return g[Integer(0) ]
- return g[Integer(0) ].simplify()
+ if isinstance(g[Integer(0)], dict):
+ return g[Integer(0)]
+ return g[Integer(0)].simplify()
if singleton and not simplify:
- return g[Integer(0) ]
+ return g[Integer(0)]
if not singleton and simplify:
G = []
@@ -3403,7 +3416,7 @@ class FractionWithFactoredDenominatorRing(
@staticmethod
def _diff_all(f, V, n, ending=[], sub=None, sub_final=None,
- zero_order=0, rekey=None):
+ zero_order=0, rekey=None):
r"""
Return a dictionary of representative mixed partial
derivatives of `f` from order 1 up to order `n` with respect to the
@@ -3454,16 +3467,16 @@ class FractionWithFactoredDenominatorRing(
D[0, 0, 0](f)(x)
sage: d1 = {diff(f, x): 4*x^3}
- sage: dd = FFDR._diff_all(f,[x], 3, sub=d1)
+ sage: dd = FFDR._diff_all(f, [x], 3, sub=d1)
sage: dd[(x, x, x)]
24*x
- sage: dd = FFDR._diff_all(f,[x], 3, sub=d1, rekey=f)
+ sage: dd = FFDR._diff_all(f, [x], 3, sub=d1, rekey=f)
sage: dd[diff(f, x, 3)]
24*x
sage: a = {x:1}
- sage: dd = FFDR._diff_all(f,[x], 3, sub=d1, rekey=f, sub_final=a)
+ sage: dd = FFDR._diff_all(f, [x], 3, sub=d1, rekey=f, sub_final=a)
sage: dd[diff(f, x, 3)]
24
@@ -3491,10 +3504,10 @@ class FractionWithFactoredDenominatorRing(
sage: dd[diff(ff, x, z)]
x*y^2*z*e^(x*y*z) + y*e^(x*y*z)
"""
- singleton=False
+ singleton = False
if not isinstance(f, list):
f = [f]
- singleton=True
+ singleton = True
# Build the dictionary of derivatives iteratively from a list
# of nondecreasing derivative-order sequences.
@@ -3521,14 +3534,13 @@ class FractionWithFactoredDenominatorRing(
for s in seeds:
value = FractionWithFactoredDenominatorRing.subs_all(
[diff(f[j], s) for j in xrange(r)], sub)
- derivs.update(dict([(tuple([j]+s), value[j])
+ derivs.update(dict([(tuple([j] + s), value[j])
for j in xrange(r)]))
for l in xrange(start, n + 1):
for t in UnorderedTuples(V, l):
s = tuple(t + ending)
value = FractionWithFactoredDenominatorRing.subs_all(
- [diff(derivs[(j,) + s[1:]],
- s[0]) for j in xrange(r)], sub)
+ [diff(derivs[(j,) + s[1:]], s[0]) for j in xrange(r)], sub)
derivs.update(dict([((j,) + s, value[j])
for j in xrange(r)]))
if zero_order:
@@ -3540,7 +3552,7 @@ class FractionWithFactoredDenominatorRing(
else:
# Ignore the first of element of k, which is an index.
for k in derivs.keys():
- if len(k) - 1 < zero_order:
+ if len(k) - 1 < zero_order:
derivs[k] = Integer(0)
if sub_final:
# Substitute sub_final into the values of derivs.
@@ -3552,12 +3564,12 @@ class FractionWithFactoredDenominatorRing(
F = rekey
if singleton:
# F must be a singleton.
- derivs = dict( [(diff(F, list(k)), derivs[k])
- for k in derivs.keys()] )
+ derivs = dict([(diff(F, list(k)), derivs[k])
+ for k in derivs.keys()])
else:
# F must be a list.
- derivs = dict( [(diff(F[k[0]], list(k)[1:]), derivs[k])
- for k in derivs.keys()] )
+ derivs = dict([(diff(F[k[0]], list(k)[1:]), derivs[k])
+ for k in derivs.keys()])
return derivs
@@ -3566,6 +3578,7 @@ class FractionWithFactoredDenominatorRing(
r"""
Return the derivatives `DD^{(l+k)}(A[j] B^l)` evaluated at a point
`p` for various natural numbers `j, k, l` which depend on `r` and `N`.
+
Here `DD` is a specific second-order linear differential operator
that depends on `M` , `A` is a list of symbolic functions,
`B` is symbolic function, and ``AB_derivs`` contains all the derivatives
@@ -3582,8 +3595,8 @@ class FractionWithFactoredDenominatorRing(
variables ``V``
- ``B`` -- a symbolic function in the variables ``V``.
- ``AB_derivs`` -- a dictionary whose keys are the (symbolic)
- derivatives of ``A[0], ..., A[r-1]`` up to order ``2*N-2`` and
- the (symbolic) derivatives of ``B`` up to order ``2*N``;
+ derivatives of ``A[0], ..., A[r-1]`` up to order ``2 * N-2`` and
+ the (symbolic) derivatives of ``B`` up to order ``2 * N``;
the values of the dictionary are complex numbers that are
the keys evaluated at a common point `p`
- ``V`` -- the variables of the ``A[j]`` and ``B``
@@ -3621,37 +3634,36 @@ class FractionWithFactoredDenominatorRing(
product_derivs = {}
for (j, k) in mrange([r, N]):
if j + k < N:
- for l in xrange(2*k + 1):
- for s in UnorderedTuples(V, 2*(k + l)):
- product_derivs[tuple([j, k, l] + s)] = \
- diff(A[j]*B**l, s).subs(AB_derivs)
+ for l in xrange(2 * k + 1):
+ for s in UnorderedTuples(V, 2 * (k + l)):
+ DF = diff(A[j] * B ** l, s).subs(AB_derivs)
+ product_derivs[tuple([j, k, l] + s)] = DF
# Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary.
DD = {}
rows = M.nrows()
for (j, k) in mrange([r, N]):
if j + k < N:
- for l in xrange(2*k + 1):
+ for l in xrange(2 * k + 1):
# Take advantage of the symmetry of M by ignoring
# the upper-diagonal entries of M and multiplying by
# appropriate powers of 2.
- if k + l == 0 :
+ if k + l == 0:
DD[(j, k, l)] = product_derivs[(j, k, l)]
continue
S = [(a, b) for (a, b) in mrange([rows, rows]) if b <= a]
- P = cartesian_product_iterator([S for i in range(k+l)])
+ P = cartesian_product_iterator([S for i in range(k + l)])
diffo = Integer(0)
for t in P:
- if product_derivs[(j, k, l) + FractionWithFactoredDenominatorRing._diff_seq(V, t)] !=\
- Integer(0):
+ if product_derivs[(j, k, l) + FractionWithFactoredDenominatorRing._diff_seq(V, t)]:
MM = Integer(1)
for (a, b) in t:
- MM = MM * M[a][b]
+ MM *= M[a][b]
if a != b:
- MM = Integer(2) *MM
- diffo = diffo + MM*product_derivs[(j, k, l) +\
- FractionWithFactoredDenominatorRing._diff_seq(V, t)]
- DD[(j, k, l)] = (-Integer(1) )**(k+l)*diffo
+ MM *= Integer(2)
+ diffo += MM * product_derivs[(j, k, l) +
+ FractionWithFactoredDenominatorRing._diff_seq(V, t)]
+ DD[(j, k, l)] = (-Integer(1)) ** (k + l) * diffo
return DD
@@ -3705,9 +3717,9 @@ class FractionWithFactoredDenominatorRing(
- ``A, B`` -- Symbolic functions in the variable ``x``
- ``AB_derivs`` - a dictionary whose keys are the (symbolic)
- derivatives of ``A`` up to order ``2*N`` if ``v`` is even or
+ derivatives of ``A`` up to order ``2 * N`` if ``v`` is even or
``N`` if ``v`` is odd and the (symbolic) derivatives of ``B``
- up to order ``2*N + v`` if ``v`` is even or ``N + v``
+ up to order ``2 * N + v`` if ``v`` is even or ``N + v``
if ``v`` is odd; the values of the dictionary are complex numbers
that are the keys evaluated at a common point `p`
- ``x`` -- a symbolic variable
@@ -3739,17 +3751,19 @@ class FractionWithFactoredDenominatorRing(
"""
I = sqrt(-Integer(1))
DD = {}
- if v.mod(Integer(2)) == Integer(0) :
+ if v.mod(Integer(2)) == Integer(0):
for k in xrange(N):
- for l in xrange(2*k + 1):
- DD[(k, l)] = (a**(-Integer(1)/v))**(2*k + v*l)*\
- diff(A*B**l, x, 2*k + v*l).subs(AB_derivs)
+ for l in xrange(2 * k + 1):
+ DD[(k, l)] = ((a ** (-Integer(1) / v)) ** (2 * k + v * l) *
+ diff(A * B ** l, x,
+ 2 * k + v * l).subs(AB_derivs))
else:
for k in xrange(N):
for l in xrange(k + 1):
- DD[(k, l)] = (abs(a)**(-Integer(1)/v)*I*\
- a/abs(a))**(k+v*l)*\
- diff(A*B**l, x, k + v*l).subs(AB_derivs)
+ DD[(k, l)] = ((abs(a) ** (-Integer(1) / v) * I *
+ a / abs(a)) ** (k + v * l) *
+ diff(A * B ** l, x,
+ k + v * l).subs(AB_derivs))
return DD
@@ -3817,18 +3831,18 @@ class FractionWithFactoredDenominatorRing(
for t in UnorderedTuples(X, l):
s = t + end
lhs.append(f_derivs[tuple(s)])
- rhs.append(diff(u*g, s).subs(atc).subs(uderivs))
+ rhs.append(diff(u * g, s).subs(atc).subs(uderivs))
# Since Sage's solve command can't take derivatives as variable
# names, make new variables based on t to stand in for
# diff(u, t) and store them in D.
- D[diff(u, t).subs(atc)] = var('zing' +\
+ D[diff(u, t).subs(atc)] = var('zing' +
''.join([str(x) for x in t]))
eqns = [lhs[i] == rhs[i].subs(uderivs).subs(D)
for i in xrange(len(lhs))]
variables = D.values()
- sol = solve(eqns,*variables, solution_dict=True)
- uderivs.update(FractionWithFactoredDenominatorRing.subs_all(
- D, sol[Integer(0) ]))
+ sol = solve(eqns, *variables, solution_dict=True)
+ uderivs.update(
+ FractionWithFactoredDenominatorRing.subs_all(D, sol[Integer(0)]))
return uderivs
@@ -3853,7 +3867,7 @@ class FractionWithFactoredDenominatorRing(
"""
if coordinate is None:
coordinate = len(v) - 1
- return tuple([vv/v[coordinate] for vv in v])
+ return tuple([vv / v[coordinate] for vv in v])
@staticmethod
@@ -4076,7 +4090,7 @@ class FFPDSum(list):
new_FFPDs = []
temp = FFPDs[0]
for f in FFPDs[1:]:
- if temp.denominator_factored() == f.denominator_factored():
+ if temp.denominator_factored() == f.denominator_factored():
# Add f to temp.
num = temp.numerator() + f.numerator()
temp = f.parent()(num, temp.denominator_factored())
@@ -4123,7 +4137,7 @@ class FFPDSum(list):
# Compute the sum's denominator factorization.
# Could use the factor() command, but it's probably faster to use
# the irreducible factors of the denominators of self.
- df = [] # The denominator factorization for the sum.
+ df = [] # The denominator factorization for the sum.
if denom == 1:
# Done
return FractionWithFactoredDenominatorRing(numer.parent())(numer, df, reduce_=False)