# Ticket #10519: trac_10519-asymptotic_multiple-v3.patch

File trac_10519-asymptotic_multiple-v3.patch, 136.0 KB (added by Frédéric Chapoton, 10 years ago)
• ## sage/combinat/all.py

# HG changeset patch
# User Alex Raichev <alex.raichev@gmail.com>
# Date 1344471313 -43200
# Node ID 9fb72a90633451c8903259b6d5cfe38252902521
# Parent  d06cf4b2215d37d3a87a58f65ac53234502dd471
#10519: asymptotics for multivariate generating functions

diff --git a/sage/combinat/all.py b/sage/combinat/all.py
 a from degree_sequences import DegreeSeque from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck from sidon_sets import sidon_sets from amgf import FFPD
• ## new file sage/combinat/amgf.py

diff --git a/sage/combinat/amgf.py b/sage/combinat/amgf.py
new file mode 100644
 - r""" Let $F(x) = \sum_{\nu \in \NN^d} F_{\nu} x^\nu$ be a multivariate power series with complex coefficients that converges in a neighborhood of the origin. Assume that $F = G/H$ for some functions $G$ and $H$ holomorphic in a neighborhood of the origin. Assume also that $H$ is a polynomial. This Python module for use within Sage _ computes asymptotics for the coefficients $F_{r \alpha}$ as $r \to \infty$ with $r \alpha \in \NN^d$ for $\alpha$ in a permissible subset of $d$-tuples of positive reals. More specifically, it computes arbitrary terms of the asymptotic expansion for $F_{r \alpha}$ when the asymptotics are controlled by a multiple point of the alegbraic variety $H = 0$. The algorithms and formulas implemented here come from [RaWi2008a]_ and [RaWi2012]_. .. [AiYu1983] I.A. Aizenberg and A.P. Yuzhakov, "Integral representations and residues in multidimensional complex analysis", Translations of Mathematical Monographs, 58. American Mathematical Society, Providence, RI, 1983. x+283 pp. ISBN: 0-8218-4511-X. .. [Raic2012] Alexander Raichev, "Leinartas's partial fraction decomposition", _. .. [RaWi2008a] Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", Electronic Journal of Combinatorics, Vol. 15 (2008), R89, _. .. [RaWi2012] Alexander Raichev and Mark C. Wilson, "Asymptotics of coefficients of multivariate generating functions: improvements for smooth points", To appear in 2012 in the Online Journal of Analytic Combinatorics, _. AUTHORS: - Alexander Raichev (2008-10-01): Initial version - Alexander Raichev (2010-09-28): Corrected many functions - Alexander Raichev (2010-12-15): Updated documentation - Alexander Raichev (2011-03-09): Fixed a division by zero bug in relative_error() - Alexander Raichev (2011-04-26): Rewrote in object-oreinted style - Alexander Raichev (2011-05-06): Fixed bug in cohomologous_integrand() and fixed _crit_cone_combo() to work in SR - Alexander Raichev (2012-08-06): Major rewrite. Created class FFPD and moved functions around. EXAMPLES:: A smooth point example (Example 5.4 of [RaWi2008a]_):: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: q = 1/2 sage: qq = q.denominator() sage: H = 1 - q*x + q*x*y - x^2*y sage: Hfac = H.factor() sage: G = (1 - q*x)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = list(qq*vector([2, 1 - q])) sage: print alpha [4, 1] sage: I = F.smooth_critical_ideal(alpha) sage: print I Ideal (y^2 - 2*y + 1, x + 1/4*y - 5/4) of Multivariate Polynomial Ring in x, y over Rational Field sage: s = solve(I.gens(), [SR(x) for x in R.gens()], solution_dict=true) sage: print s [{y: 1, x: 1}] sage: p = s[0] sage: asy = F.asymptotics(p, alpha, 1) # long time Creating auxiliary functions... Computing derivatives of auxiallary functions... Computing derivatives of more auxiliary functions... Computing second order differential operator actions... sage: print asy # long time (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)), 1, 1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3))) sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1]) # long time Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [((4, 1), 0.1875000000, [0.1953794675], [-0.04202382689]), ((8, 2), 0.1523437500, [0.1550727862], [-0.01791367323]), ((16, 4), 0.1221771240, [0.1230813519], [-0.007400959228]), ((32, 8), 0.09739671811, [0.09768973377], [-0.003008475766]), ((64, 16), 0.07744253816, [0.07753639308], [-0.001211929722])] A multiple point example (Example 6.5 of [RaWi2012]_):: sage: R.= PolynomialRing(QQ) sage: H = (1 - 2*x - y)**2 * (1 - x - 2*y)**2 sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F (1, [(x + 2*y - 1, 2), (2*x + y - 1, 2)]) sage: I = F.singular_ideal() sage: print I Ideal (x - 1/3, y - 1/3) of Multivariate Polynomial Ring in x, y over Rational Field sage: p = {x: 1/3, y: 1/3} sage: print F.is_convenient_multiple_point(p) (True, 'convenient in variables [x, y]') sage: alpha = (var('a'), var('b')) sage: print F.asymptotic_decomposition(alpha) # long time [(0, []), (-1/9*(2*a^2*y^2 - 5*a*b*x*y + 2*b^2*x^2)*r^2/(x^2*y^2) + 1/9*(5*(a + b)*x*y - 6*a*y^2 - 6*b*x^2)*r/(x^2*y^2) - 1/9*(4*x^2 - 5*x*y + 4*y^2)/(x^2*y^2), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])] sage: print F.asymptotics(p, alpha, 2) # long time (-3*((2*a^2 - 5*a*b + 2*b^2)*r^2 + (a + b)*r + 3)*((1/3)^(-b)*(1/3)^(-a))^r, (1/3)^(-b)*(1/3)^(-a), -3*(2*a^2 - 5*a*b + 2*b^2)*r^2 - 3*(a + b)*r - 9) sage: alpha = [4, 3] sage: asy = F.asymptotics(p, alpha, 2) # long time sage: print asy # long time (3*(10*r^2 - 7*r - 3)*2187^r, 2187, 30*r^2 - 21*r - 9) sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1]) # long time Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [((4, 3), 30.72702332, [0.0000000000], [1.000000000]), ((8, 6), 111.9315678, [69.00000000], [0.3835519207]), ((16, 12), 442.7813138, [387.0000000], [0.1259793763]), ((32, 24), 1799.879232, [1743.000000], [0.03160169385])] """ #***************************************************************************** #       Copyright (C) 2008 Alexander Raichev # #  Distributed under the terms of the GNU General Public License (GPL) #                  http://www.gnu.org/licenses/ #***************************************************************************** from functools import total_ordering # Sage libraries from sage.categories.unique_factorization_domains   import UniqueFactorizationDomains from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.polynomial.polynomial_ring  import is_PolynomialRing from sage.rings.polynomial.multi_polynomial_ring_generic    import is_MPolynomialRing from sage.symbolic.ring import SR from sage.geometry.cone import Cone from sage.calculus.functional   import diff from sage.calculus.functions    import jacobian from sage.calculus.var  import function, var from sage.combinat.cartesian_product    import CartesianProduct from sage.combinat.combinat     import stirling_number1 from sage.combinat.permutation  import Permutation from sage.combinat.tuple    import UnorderedTuples from sage.functions.log     import exp, log from sage.functions.other   import factorial, gamma, sqrt from sage.matrix.constructor    import matrix from sage.misc.misc import add from sage.misc.misc_c   import prod from sage.misc.mrange   import cartesian_product_iterator, mrange from sage.modules.free_module_element   import vector from sage.rings.arith   import binomial from sage.rings.all import CC from sage.rings.fraction_field  import FractionField from sage.rings.integer import Integer from sage.rings.integer_ring    import ZZ from sage.rings.polynomial.polynomial_ring_constructor  import PolynomialRing from sage.rings.rational_field  import QQ from sage.sets.set  import Set from sage.structure.sage_object import SageObject from sage.symbolic.constants    import pi from sage.symbolic.relation import solve @total_ordering class FFPD(object): r""" Represents a fraction with factored polynomial denominator (FFPD) $p/(q_1^{e_1} \cdots q_n^{e_n})$ by storing the parts $p$ and $[(q_1, e_1), \ldots, (q_n, e_n)]$. Here $q_1, \ldots, q_n$ are elements of a 0- or multi-variate factorial polynomial ring $R$ , $q_1, \ldots, q_n$ are distinct irreducible elements of $R$ , $e_1, \ldots, e_n$ are positive integers, and $p$ is a function of the indeterminates of $R$ (a Sage Symbolic Expression). An element $r$ with no polynomial denominator is represented as $[r, (,)]$. AUTHORS: - Alexander Raichev (2012-07-26) """ def __init__(self, numerator=None, denominator_factored=None, quotient=None, reduce_=True): r""" Create a FFPD instance. INPUT: - numerator - (Optional; default=None) An element $p$ of a 0- or 1-variate factorial polynomial ring $R$. - denominator_factored - (Optional; default=None) A list of the form $[(q_1, e_1), \ldots, (q_n, e_n)]$ where the $q_1, \ldots, q_n$ are distinct irreducible elements of $R$ and the $e_i$ are positive integers. - quotient - (Optional; default=None) An element of a field of fractions of a factorial ring. - reduce_ - (Optional; default=True) If True, then represent $p/(q_1^{e_1} \cdots q_n^{e_n})$ in lowest terms. If False, then won't attempt to divide $p$ by any of the $q_i$. OUTPUT: A FFPD instance representing the rational expression $p/(q_1^{e_1} \cdots q_n^{e_n})$. To get a non-None output, one of numerator or quotient must be non-None. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: df = [x, 1], [y, 1], [x*y+1, 1] sage: f = FFPD(x, df) sage: print f (1, [(y, 1), (x*y + 1, 1)]) sage: ff = FFPD(x, df, reduce_=False) sage: print ff (x, [(y, 1), (x, 1), (x*y + 1, 1)]) :: sage: f = FFPD(x + y, [(x + y, 1)]) sage: print f (1, []) :: sage: R. = PolynomialRing(QQ) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print FFPD(quotient=f) (5*x^7 - 5*x^6 + 5/3*x^5 - 5/3*x^4 + 2*x^3 - 2/3*x^2 + 1/3*x - 1/3, [(x - 1, 1), (x, 1), (x^2 + 1/3, 1)]) :: sage: R. = PolynomialRing(QQ) sage: f = 2*y/(5*(x^3 - 1)*(y + 1)) sage: print FFPD(quotient=f) (2/5*y, [(y + 1, 1), (x - 1, 1), (x^2 + x + 1, 1)]) :: sage: R.= PolynomialRing(QQ) sage: p = 1/x^2 sage: q = 3*x**2*y sage: qs = q.factor() sage: f = FFPD(p/qs.unit(), qs) sage: print f (1/(3*x^2), [(y, 1), (x, 2)]) :: sage: R. = PolynomialRing(QQ) sage: f = FFPD(cos(x)*x*y^2, [(x, 2), (y, 1)]) sage: print f (x*y^2*cos(x), [(y, 1), (x, 2)]) :: sage: R. = PolynomialRing(QQ) sage: G = exp(x + y) sage: H = (1 - 2*x - y) * (1 - x - 2*y) sage: a = FFPD(quotient=G/H) sage: print a (e^(x + y)/(2*x^2 + 5*x*y + 2*y^2 - 3*x - 3*y + 1), []) sage: print a._ring None sage: b = FFPD(G, H.factor()) sage: print b (e^(x + y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)]) sage: print b._ring Multivariate Polynomial Ring in x, y over Rational Field Singular throws a 'not implemented' error when trying to factor in a multivariate polynomial ring over an inexact field :: sage: R.= PolynomialRing(CC) sage: f = (x + 1)/(x*y*(x*y + 1)^2) sage: FFPD(quotient=f) Traceback (most recent call last): ... TypeError: Singular error: ? not implemented ? error occurred in or before STDIN line 17: def sage9=factorize(sage8); """ # Attributes are # self._numerator # self._denominator_factored # self._ring if quotient is not None: p = quotient.numerator() q = quotient.denominator() R = q.parent() self._numerator = quotient self._denominator_factored = [] if is_PolynomialRing(R) or is_MPolynomialRing(R): self._ring = R if not R(q).is_unit(): # Factor q try: df = q.factor() except NotImplementedError: # Singular's factor() needs 'proof=False'. df = q.factor(proof=False) self._numerator = p/df.unit() df = sorted([tuple(t) for t in df]) # Sort for consitency. self._denominator_factored = df else: self._ring = None # Done. No reducing needed, as Sage reduced quotient beforehand. return self._numerator = numerator if denominator_factored: self._denominator_factored = sorted([tuple(t) for t in denominator_factored]) self._ring = denominator_factored[0][0].parent() else: self._denominator_factored = [] self._ring = None R = self._ring if R is not None and numerator in R and reduce_: # Reduce fraction if possible. numer = R(self._numerator) df = self._denominator_factored new_df = [] for (q, e) in df: ee = e quo, rem = numer.quo_rem(q) while rem == 0 and ee > 0: ee -= 1 numer = quo quo, rem = numer.quo_rem(q) if ee > 0: new_df.append((q, ee)) self._numerator = numer self._denominator_factored = new_df def numerator(self): r""" Return the numerator of self. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F.numerator() -e^y """ return self._numerator def denominator(self): r""" Return the denominator of self. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F.denominator() 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()]) def denominator_factored(self): r""" Return the factorization in self.ring() of the denominator of self but without the unit part. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F.denominator_factored() [(x - 1, 1), (x*y + x + y - 1, 2)] """ return self._denominator_factored def ring(self): r""" Return the ring of the denominator of self, which is None in the case where self doesn't have a denominator. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F.ring() Multivariate Polynomial Ring in x, y over Rational Field sage: F = FFPD(quotient=G/H) sage: print F (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), []) sage: print F.ring() None """ return self._ring def dimension(self): r""" Return the number of indeterminates of self.ring(). EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F.dimension() 2 """ R = self.ring() if is_PolynomialRing(R) or is_MPolynomialRing(R): return R.ngens() else: return None def list(self): r""" Convert self into a list for printing. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F # indirect doctest (-e^y, [(x - 1, 1), (x*y + x + y - 1, 2)]) """ return (self.numerator(), self.denominator_factored()) def quotient(self): r""" Convert self into a quotient. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2*(1-x) sage: Hfac = H.factor() sage: G = exp(y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F (-e^y, [(x - 1, 1), (x*y + x + y - 1, 2)]) sage: print F.quotient() -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() def __str__(self): r""" Returns a string representation of self EXAMPLES:: """ return str(self.list()) def __eq__(self, other): r""" Two FFPD instances are equal iff they represent the same fraction. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: df = [x, 1], [y, 1], [x*y+1, 1] sage: f = FFPD(x, df) sage: ff = FFPD(x, df, reduce_=False) sage: f == ff True sage: g = FFPD(y, df) sage: g == f False :: sage: R. = PolynomialRing(QQ) sage: G = exp(x + y) sage: H = (1 - 2*x - y) * (1 - x - 2*y) sage: a = FFPD(quotient=G/H) sage: b = FFPD(G, H.factor()) sage: bool(a == b) True """ return self.quotient() == other.quotient() def __ne__(self, other): r""" EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: df = [x, 1], [y, 1], [x*y+1, 1] sage: f = FFPD(x, df) sage: ff = FFPD(x, df, reduce_=False) sage: f != ff False sage: g = FFPD(y, df) sage: g != f # indirect doctest True """ return not (self == other) def __lt__(self, other): r""" FFPD A is less than FFPD B iff (the denominator factorization of A is shorter than that of B) or (the denominator factorization lengths are equal and the denominator of A is less than that of B in their ring) or (the denominator factorization lengths are equal and the denominators are equal and the numerator of A is less than that of B in their ring). EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: df = [x, 1], [y, 1], [x*y+1, 1] sage: f = FFPD(x, df) sage: ff = FFPD(x, df, reduce_=False) sage: g = FFPD(y, df) sage: h = FFPD(exp(x), df) sage: i = FFPD(sin(x + 2), df) sage: print f, ff (1, [(y, 1), (x*y + 1, 1)]) (x, [(y, 1), (x, 1), (x*y + 1, 1)]) sage: print f < ff True sage: print f < g True sage: print g < h True sage: print h < i False """ sn = self.numerator() on = other.numerator() sdf = self.denominator_factored() odf = other.denominator_factored() 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)) def univariate_decomposition(self): r""" Return the usual univariate partial fraction decomposition of self as a FFPDSum instance. Assume that self lies in the field of fractions of a univariate factorial polynomial ring. EXAMPLES:: sage: from sage.combinat.amgf import * One variable:: sage: R. = PolynomialRing(QQ) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print f (15*x^7 - 15*x^6 + 5*x^5 - 5*x^4 + 6*x^3 - 2*x^2 + x - 1)/(3*x^4 - 3*x^3 + x^2 - x) sage: decomp = FFPD(quotient=f).univariate_decomposition() sage: print decomp [(5*x^3, []), (1, [(x - 1, 1)]), (1, [(x, 1)]), (1/3, [(x^2 + 1/3, 1)])] sage: print decomp.sum().quotient() == f True One variable with numerator in symbolic ring:: sage: R. = PolynomialRing(QQ) sage: f = 5*x^3 + 1/x + 1/(x-1) + exp(x)/(3*x^2 + 1) sage: print f e^x/(3*x^2 + 1) + ((5*(x - 1)*x^3 + 2)*x - 1)/((x - 1)*x) sage: decomp = FFPD(quotient=f).univariate_decomposition() sage: print decomp [(e^x/(3*x^2 + 1) + ((5*(x - 1)*x^3 + 2)*x - 1)/((x - 1)*x), [])] One variable over a finite field:: sage: R. = PolynomialRing(GF(2)) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print f (x^6 + x^4 + 1)/(x^3 + x) sage: decomp = FFPD(quotient=f).univariate_decomposition() sage: print decomp [(x^3, []), (1, [(x, 1)]), (x, [(x + 1, 2)])] sage: print decomp.sum().quotient() == f True One variable over an inexact field:: sage: R. = PolynomialRing(CC) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print f (15.0000000000000*x^7 - 15.0000000000000*x^6 + 5.00000000000000*x^5 - 5.00000000000000*x^4 + 6.00000000000000*x^3 - 2.00000000000000*x^2 + x - 1.00000000000000)/(3.00000000000000*x^4 - 3.00000000000000*x^3 + x^2 - x) sage: decomp = FFPD(quotient=f).univariate_decomposition() sage: print decomp [(5.00000000000000*x^3, []), (1.00000000000000, [(x - 1.00000000000000, 1)]), (-0.288675134594813*I, [(x - 0.577350269189626*I, 1)]), (1.00000000000000, [(x, 1)]), (0.288675134594813*I, [(x + 0.577350269189626*I, 1)])] sage: print decomp.sum().quotient() == f # Rounding error coming False NOTE:: Let $f = p/q$ be a rational expression where $p$ and $q$ lie in a univariate factorial polynomial ring $R$. Let $q_1^{e_1} \cdots q_n^{e_n}$ be the unique factorization of $q$ in $R$ into irreducible factors. Then $f$ can be written uniquely as (*)   $p_0 + \sum_{i=1}^{m} p_i/ q_i^{e_i}$, for some $p_j \in R$. I call (*) the *usual partial fraction decomposition* of f. AUTHORS: - Robert Bradshaw (2007-05-31) - Alexander Raichev (2012-06-25) """ if self.dimension() is None or self.dimension() > 1: return FFPDSum([self]) R = self.ring() p = self.numerator() q = self.denominator() if p in R: whole, p = p.quo_rem(q) else: whole = p p = R(1) df = self.denominator_factored() decomp = [FFPD(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) # The inverse exists because the product and a**m # are relatively prime. decomp.append(FFPD(numer, [(a, m)])) return FFPDSum(decomp) def nullstellensatz_certificate(self): r""" Let $[(q_1, e_1), \ldots, (q_n, e_n)]$ be the denominator factorization of self. Return a list of polynomials $h_1, \ldots, h_m$ in self.ring() that satisfies $h_1 q_1 + \cdots + h_m q_n = 1$ if it exists. Otherwise return None. Only works for multivariate self. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: G = sin(x) sage: H = x^2 * (x*y + 1) sage: f = FFPD(G, H.factor()) sage: L = f.nullstellensatz_certificate() sage: print L [y^2, -x*y + 1] sage: df = f.denominator_factored() sage: sum([L[i]*df[i][0]**df[i][1] for i in xrange(len(df))]) == 1 True :: sage: f = 1/(x*y) sage: L = FFPD(quotient=f).nullstellensatz_certificate() sage: L is None True """ R = self.ring() if R is None: 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) else: return None def nullstellensatz_decomposition(self): r""" Return a Nullstellensatz decomposition of self as a FFPDSum instance. Recursive. Only works for multivariate self. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: f = 1/(x*(x*y + 1)) sage: decomp = FFPD(quotient=f).nullstellensatz_decomposition() sage: print decomp [(0, []), (1, [(x, 1)]), (-y, [(x*y + 1, 1)])] sage: decomp.sum().quotient() == f True sage: for r in decomp: ...       L = r.nullstellensatz_certificate() ...       L is None ... True True True :: sage: R. = PolynomialRing(QQ) sage: G = sin(y) sage: H = x*(x*y + 1) sage: f = FFPD(G, H.factor()) sage: decomp = f.nullstellensatz_decomposition() sage: print decomp [(0, []), (sin(y), [(x, 1)]), (-y*sin(y), [(x*y + 1, 1)])] sage: bool(decomp.sum().quotient() == G/H) True sage: for r in decomp: ...       L = r.nullstellensatz_certificate() ...       L is None ... True True True NOTE:: Let $f = p/q$ where $q$ lies in a $d$ -variate polynomial ring $K[X]$ for some field $K$ and $d \ge 1$. Let $q_1^{e_1} \cdots q_n^{e_n}$ be the unique factorization of $q$ in $K[X]$ into irreducible factors and let $V_i$ be the algebraic variety $\{x \in L^d: q_i(x) = 0\}$ of $q_i$ over the algebraic closure $L$ of $K$. By [Raic2012]_, $f$ can be written as (*)   $\sum p_A/\prod_{i \in A} q_i^{e_i}$, where the $p_A$ are products of $p$ and elements in $K[X]$ and the sum is taken over all subsets $A \subseteq \{1, \ldots, m\}$ such that $\cap_{i\in A} T_i \neq \emptyset$. I call (*) a *Nullstellensatz decomposition* of $f$. Nullstellensatz decompositions are not unique. The algorithm used comes from [Raic2012]_. """ L = self.nullstellensatz_certificate() if L is None: # No decomposing possible. return FFPDSum([self]) # Otherwise decompose recursively. decomp = FFPDSum() p = self.numerator() df = self.denominator_factored() m = len(df) iteration1 = FFPDSum([FFPD(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. for r in iteration1: decomp.extend(r.nullstellensatz_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() def algebraic_dependence_certificate(self): r""" Return the ideal $J$ of annihilating polynomials for the set of polynomials [q**e for (q, e) in self.denominator_factored()], which could be the zero ideal. The ideal $J$ lies in a polynomial ring over the field self.ring().base_ring() that has m = len(self.denominator_factored()) indeterminates. Return None if self.ring() is None. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: f = 1/(x^2 * (x*y + 1) * y^3) sage: ff = FFPD(quotient=f) sage: J = ff.algebraic_dependence_certificate() sage: print J Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 - 6*T2^5  + T2^6) of Multivariate Polynomial Ring in T0, T1, T2 over Rational Field sage: g = J.gens()[0] sage: df = ff.denominator_factored() sage: g(*(q**e for q, e in df)) == 0 True :: sage: R. = PolynomialRing(QQ) sage: G = exp(x + y) sage: H = x^2 * (x*y + 1) * y^3 sage: ff = FFPD(G, H.factor()) sage: J = ff.algebraic_dependence_certificate() sage: print J Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 - 6*T2^5 + T2^6) of Multivariate Polynomial Ring in T0, T1, T2 over Rational Field sage: g = J.gens()[0] sage: df = ff.denominator_factored() sage: g(*(q**e for q, e in df)) == 0 True :: sage: f = 1/(x^3 * y^2) sage: J = FFPD(quotient=f).algebraic_dependence_certificate() sage: print J Ideal (0) of Multivariate Polynomial Ring in T0, T1 over Rational Field :: sage: f = sin(1)/(x^3 * y^2) sage: J = FFPD(quotient=f).algebraic_dependence_certificate() sage: print J None """ R = self.ring() if R is None: return None df = self.denominator_factored() if not df: return R.ideal()    # The zero ideal. m = len(df) F = R.base_ring() Xs = list(R.gens()) d = len(Xs) # Expand R by 2*m new variables. S = 'S' while S in [str(x) for x in Xs]: S = S + 'S' Ss = [S + str(i) for i in xrange(m)] T = 'T' while T in [str(x) for x in Xs]: T = T + 'T' Ts = [T + str(i) for i in xrange(m)] Vs = [str(x) for x in Xs] + Ss + Ts RR = PolynomialRing(F, Vs) Xs = RR.gens()[:d] Ss = RR.gens()[d: d + m] 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 = 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') return RRR.ideal(J) def algebraic_dependence_decomposition(self, whole_and_parts=True): r""" Return an algebraic dependence decomposition of self as a FFPDSum instance. Recursive. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: f = 1/(x^2 * (x*y + 1) * y^3) sage: ff = FFPD(quotient=f) sage: decomp = ff.algebraic_dependence_decomposition() sage: print decomp [(0, []), (-x, [(x*y + 1, 1)]), (x^2*y^2 - x*y + 1, [(y, 3), (x, 2)])] sage: print decomp.sum().quotient() == f True sage: for r in decomp: ...       J = r.algebraic_dependence_certificate() ...       J is None or J == J.ring().ideal()  # The zero ideal ... True True True :: sage: R. = PolynomialRing(QQ) sage: G = sin(x) sage: H = x^2 * (x*y + 1) * y^3 sage: f = FFPD(G, H.factor()) sage: decomp = f.algebraic_dependence_decomposition() sage: print decomp [(0, []), (x^4*y^3*sin(x), [(x*y + 1, 1)]), (-(x^5*y^5 - x^4*y^4 + x^3*y^3 - x^2*y^2 + x*y - 1)*sin(x), [(y, 3), (x, 2)])] sage: bool(decomp.sum().quotient() == G/H) True sage: for r in decomp: ...       J = r.algebraic_dependence_certificate() ...       J is None or J == J.ring().ideal() ... True True True NOTE:: Let $f = p/q$ where $q$ lies in a $d$ -variate polynomial ring $K[X]$ for some field $K$. Let $q_1^{e_1} \cdots q_n^{e_n}$ be the unique factorization of $q$ in $K[X]$ into irreducible factors and let $V_i$ be the algebraic variety $\{x\in L^d: q_i(x) = 0\}$ of $q_i$ over the algebraic closure $L$ of $K$. By [Raic2012]_, $f$ can be written as (*)   $\sum p_A/\prod_{i \in A} q_i^{b_i}$, where the $b_i$ are positive integers, each $p_A$ is a products of $p$ and an element in $K[X]$, and the sum is taken over all subsets $A \subseteq \{1, \ldots, m\}$ such that $|A| \le d$ and $\{q_i : i\in A\}$ is algebraically independent. I call (*) an *algebraic dependence decomposition* of $f$. Algebraic dependence decompositions are not unique. The algorithm used comes from [Raic2012]_. """ J = self.algebraic_dependence_certificate() if not J: # No decomposing possible. return FFPDSum([self]) # Otherwise decompose recursively. decomp = FFPDSum() p = self.numerator() df = self.denominator_factored() m = len(df) g = J.gens()[0]     # An annihilating polynomial for df. new_vars = J.ring().gens() # Note that each new_vars[j] corresponds to df[j] such that # g([q**e for q, e in df]) = 0. # Assuming here that g.parent() has negdeglex term order # 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()) numers = map(prod, zip(gg.coefficients(), gg.monomials())) 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. 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]) for j in xrange(m)]) iteration1 = FFPDSum() for r in iteration1_temp: num1 = p*g.parent()(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)) iteration1.append(FFPD(num1, denoms1)) # Now decompose each FFPD of iteration1. for r in iteration1: decomp.extend(r.algebraic_dependence_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() def leinartas_decomposition(self): r""" Return a Leinartas decomposition of self as a FFPDSum instance. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: f = 1/x + 1/y + 1/(x*y + 1) sage: decomp = FFPD(quotient=f).leinartas_decomposition() sage: print decomp [(0, []), (1, [(x*y + 1, 1)]), (x + y, [(y, 1), (x, 1)])] sage: print decomp.sum().quotient() == f True sage: for r in decomp: ...       L = r.nullstellensatz_certificate() ...       print L is None ...       J = r.algebraic_dependence_certificate() ...       print J is None or J == J.ring().ideal() ... True True True True True True :: sage: R. = PolynomialRing(QQ) sage: f = sin(x)/x + 1/y + 1/(x*y + 1) sage: G = f.numerator() sage: H = R(f.denominator()) sage: ff = FFPD(G, H.factor()) sage: decomp = ff.leinartas_decomposition() sage: print decomp [(0, []), (-(x*y^2*sin(x) + x^2*y + x*y + y*sin(x) + x)*y, [(y, 1)]), ((x*y^2*sin(x) + x^2*y + x*y + y*sin(x) + x)*x*y, [(x*y + 1, 1)]), (x*y^2*sin(x) + x^2*y + x*y + y*sin(x) + x, [(y, 1), (x, 1)])] sage: bool(decomp.sum().quotient() == f) True sage: for r in decomp: ...       L = r.nullstellensatz_certificate() ...       print L is None ...       J = r.algebraic_dependence_certificate() ...       print J is None or J == J.ring().ideal() ... True True True True True True True True :: sage: R.= PolynomialRing(GF(2, 'a')) sage: f = 1/(x * y * z * (x*y + z)) sage: decomp = FFPD(quotient=f).leinartas_decomposition() sage: print decomp [(0, []), (1, [(z, 2), (x*y + z, 1)]), (1, [(z, 2), (y, 1), (x, 1)])] sage: print decomp.sum().quotient() == f True NOTE:: Let $f = p/q$ where $q$ lies in a $d$ -variate polynomial ring $K[X]$ for some field $K$. Let $q_1^{e_1} \cdots q_n^{e_n}$ be the unique factorization of $q$ in $K[X]$ into irreducible factors and let $V_i$ be the algebraic variety $\{x\in L^d: q_i(x) = 0\}$ of $q_i$ over the algebraic closure $L$ of $K$. By [Raic2012]_, $f$ can be written as (*)   $\sum p_A/\prod_{i \in A} q_i^{b_i}$, where the $b_i$ are positive integers, each $p_A$ is a product of $p$ and an element of $K[X]$, and the sum is taken over all subsets $A \subseteq \{1, \ldots, m\}$ such that (1) $|A| \le d$, (2) $\cap_{i\in A} T_i \neq \emptyset$, and (3) $\{q_i : i\in A\}$ is algebraically independent. In particular, any rational expression in $d$ variables can be represented as a sum of rational expressions whose denominators each contain at most $d$ distinct irreducible factors. I call (*) a *Leinartas decomposition* of $f$. Leinartas decompositions are not unique. The algorithm used comes from [Raic2012]_. """ d = self.dimension() if d == 1: # Sage's lift() function doesn't work in # univariate polynomial rings. # So nullstellensatz_decomposition() won't work. # So use algebraic_dependence_decomposition(), # which is sufficient. temp = FFPDSum([self]) else: temp = self.nullstellensatz_decomposition() decomp = FFPDSum() for r in temp: decomp.extend(r.algebraic_dependence_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() def cohomology_decomposition(self): r""" Let $p/(q_1^{e_1} \cdots q_n^{e_n})$ be the fraction represented by self and let $K[x_1, \ldots, x_d]$ be the polynomial ring in which the $q_i$ lie. Assume that $n \le d$ and that the gradients of the $q_i$ are linearly independent at all points in the intersection $V_1 \cap \ldots \cap V_n$ of the algebraic varieties $V_i = \{x \in L^d: q_i(x) = 0 \}$, where $L$ is the algebraic closure of the field $K$. Return a FFPDSum $f$ such that the differential form $f dx_1 \wedge \cdots \wedge dx_d$ is de Rham cohomologous to the differential form $p/(q_1^{e_1} \cdots q_n^{e_n}) dx_1 \wedge \cdots \wedge dx_d$ and such that the denominator of each summand of $f$ contains no repeated irreducible factors. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: print FFPD(1, [(x, 1), (y, 2)]).cohomology_decomposition() [(0, [])] :: sage: R.= PolynomialRing(QQ) sage: p = 1 sage: qs = [(x*y - 1, 1), (x**2 + y**2 - 1, 2)] sage: f = FFPD(p, qs) sage: print f.cohomology_decomposition() [(0, []), (4/3*x*y + 4/3, [(x^2 + y^2 - 1, 1)]), (1/3, [(x*y - 1, 1), (x^2 + y^2 - 1, 1)])] NOTES: The algorithm used here comes from the proof of Theorem 17.4 of [AiYu1983]_. AUTHORS: - Alexander Raichev (2008-10-01, 2011-01-15, 2012-07-31) """ R = self.ring() df = self.denominator_factored() n = len(df) if R is None or sum([e for (q, e) in df]) <= n: # No decomposing possible. return FFPDSum([self]) # Otherwise decompose recursively. decomp = FFPDSum() p = self.numerator() qs = [q for (q, e) in df] X = sorted(R.gens()) var_sets_n = Set(X).subsets(n) # Compute Jacobian determinants for qs. dets = [] for v in var_sets_n: # Sort v according to the term order of R. x = sorted(v) jac = jacobian(qs, x) dets.append(R(jac.determinant())) # Get a Nullstellensatz certificate for qs and dets. L = R(1).lift(R.ideal(qs + dets)) # Do first iteration of decomposition. iteration1 = FFPDSum() # Contributions from qs. for i in xrange(n): if L[i] == 0: continue # Cancel one df[i] from denominator. new_df = [list(t) for t in df] if new_df[i][1] > 1: new_df[i][1] -= 1 else: del(new_df[i]) iteration1.append(FFPD(p*L[i], new_df)) # Contributions from dets. # Compute each contribution's cohomologous form using # the least index j such that new_df[j][1] > 1. # Know such an index exists by first 'if' statement at # the top. for j in xrange(n): if df[j][1] > 1: J = j break new_df = [list(t) for t in df] new_df[J][1] -= 1 for k in xrange(var_sets_n.cardinality()): if L[n + k] == 0: continue # 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])] + [SR(qs[j]) for j in xrange(n) if j != J], [SR(xx) for xx in x]) det = jac.determinant() psign = FFPD.permutation_sign(x, X) iteration1.append(FFPD((-1)**J*det/\ (psign*new_df[J][1]), new_df)) # Now decompose each FFPD of iteration1. for r in iteration1: decomp.extend(r.cohomology_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() @staticmethod def permutation_sign(s, u): r""" This function returns the sign of the permutation on 1, ..., len(u) that is induced by the sublist s of u. For internal use by cohomology_decomposition(). INPUT: - s - A sublist of u. - u - A list. OUTPUT: The sign of the permutation obtained by taking indices within u of the list s + sc, where sc is u with the elements of s removed. EXAMPLES:: sage: from sage.combinat.amgf import * sage: u = ['a','b','c','d','e'] sage: s = ['b','d'] sage: FFPD.permutation_sign(s, u) -1 sage: s = ['d','b'] sage: FFPD.permutation_sign(s, u) 1 AUTHORS: - Alexander Raichev (2008-10-01, 2012-07-31) """ # Convert lists to lists of numbers in {1,..., 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)))) P = Permutation(B + C) return P.signature() def asymptotic_decomposition(self, alpha, asy_var=None): r""" Return a FFPDSum that has the same asymptotic expansion as self in the direction alpha but each of whose FFPDs has a denominator factorization of the form $[(q_1, 1), \ldots, (q_n, 1)]$, where n is at most d = self.dimension(). The output results from a Leinartas decomposition followed by a cohomology decomposition. INPUT: - alpha - A d-tuple of positive integers or symbolic variables. - asy_var - (Optional; default=None) A symbolic variable with respect to which to compute asymptotics. If None is given the set asy_var = var('r'). EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - 2*x -y)*(1 - x -2*y)**2 sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = var('a, b') sage: r = var('r') sage: print F.asymptotic_decomposition(alpha, r) # long time [(0, []), (-1/3*(a*y - 2*b*x)*r/(x*y) + 1/3*(2*x - y)/(x*y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])] AUTHORS: - Alexander Raichev (2012-08-01) """ R = self.ring() if R is None: return None d = self.dimension() n = len(self.denominator_factored()) X = [SR(x) for x in R.gens()] # Reduce number of distinct factors in denominator of self # down to at most d. decomp1 = FFPDSum([self]) if n > d: decomp1 = decomp1[0].leinartas_decomposition() # Reduce to no repeated factors in denominator of each element # of decomp1. # Compute the cohomology decomposition for each # 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)]) decomp2 = FFPDSum() for f in decomp1: ff = FFPD(f.numerator()*cauchy_stuff, f.denominator_factored()) decomp2.extend(ff.cohomology_decomposition()) decomp2 = decomp2.combine_like_terms() # Divide out cauchy_stuff from integrands. decomp3 = FFPDSum() for f in decomp2: ff = FFPD((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): r""" Return the first $N$ terms (some of which could be zero) of the asymptotic expansion of the Maclaurin ray coefficients $F_{r\alpha}$ of the function $F$ represented by self as $r\to\infty$, where $r$ = asy_var and alpha is a tuple of positive integers of length d = self.dimension(). Assume that $F$ is holomorphic in a neighborhood of the origin, that the denominator factorization of self is also the unique factorization of the denominator of $F$ in the local analytic ring at $p$ (not just in the polynomial ring self.ring()), that $p$ is a convenient strictly minimal smooth or multiple point with all nonzero coordinates that is critical and nondegenerate for alpha. INPUT: - p - A dictionary with keys that can be coerced to equal self.ring().gens(). - alpha - A tuple of length self.dimension() of positive integers or, if $p$ is a smooth point, possibly of symbolic variables. - N - A positive integer. - numerical - (Optional; default=0) A natural number. If numerical > 0, then return a numerical approximation of $F_{r \alpha}$ with numerical digits of precision. Otherwise return exact values. - asy_var - (Optional; default=None) A  symbolic variable. The variable of the asymptotic expansion. If none is given, var('r') will be assigned. OUTPUT: The tuple (asy, exp_scale, subexp_part). Here asy is the sum of the first $N$ terms (some of which might be 0) of the asymptotic expansion of $F_{r\alpha}$ as $r\to\infty$; exp_scale**r is the exponential factor of asy; subexp_part is the subexponential factor of asy. EXAMPLES:: sage: from sage.combinat.amgf import * A smooth point example:: sage: R.= PolynomialRing(QQ) sage: H = (1 - x - y - x*y)**2 sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac);print(F) (1, [(x*y + x + y - 1, 2)]) sage: alpha = [4, 3] sage: p = {y: 1/3, x: 1/2} sage: asy = F.asymptotics(p, alpha, 2) # long time Creating auxiliary functions... Computing derivatives of auxiallary functions... Computing derivatives of more auxiliary functions... Computing second order differential operator actions... sage: print asy # long time (1/6000*(3600*sqrt(2)*sqrt(3)*sqrt(5)*sqrt(r)/sqrt(pi) + 463*sqrt(2)*sqrt(3)*sqrt(5)/(sqrt(pi)*sqrt(r)))*432^r, 432, 1/6000*(3600*sqrt(5)*r + 463*sqrt(5))*sqrt(2)*sqrt(3)/(sqrt(pi)*sqrt(r))) sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1]) # long time Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [((4, 3), 2.083333333, [2.092576110], [-0.004436533009]), ((8, 6), 2.787374614, [2.790732875], [-0.001204811281]), ((16, 12), 3.826259447, [3.827462310], [-0.0003143703383]), ((32, 24), 5.328112821, [5.328540787], [-0.00008032229296]), ((64, 48), 7.475927885, [7.476079664], [-0.00002030233658])] A multiple point example:: sage: R.= PolynomialRing(QQ) sage: H = (4 - 2*x - y - z)**2*(4 - x - 2*y - z) sage: Hfac = H.factor() sage: G = 16/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F (-16, [(x + 2*y + z - 4, 1), (2*x + y + z - 4, 2)]) sage: alpha = [3, 3, 2] sage: p = {x: 1, y: 1, z: 1} sage: asy = F.asymptotics(p, alpha, 2) # long time Creating auxiliary functions... Computing derivatives of auxiliary functions... Computing derivatives of more auxiliary functions... Computing second-order differential operator actions... sage: print asy # long time (4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(r)), 1, 1/216*(288*sqrt(3)*r + 47*sqrt(3))/(sqrt(pi)*sqrt(r))) sage: print F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1]) # long time Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [((3, 3, 2), 0.9812164307, [1.515572606], [-0.5445854340]), ((6, 6, 4), 1.576181132, [1.992989399], [-0.2644418580]), ((12, 12, 8), 2.485286378, [2.712196351], [-0.09130133851]), ((24, 24, 16), 3.700576827, [3.760447895], [-0.01617884750])] NOTES: The algorithms used here come from [RaWi2008a]_ and [RaWi2012]_. AUTHORS: - Alexander Raichev (2008-10-01, 2010-09-28, 2011-04-27, 2012-08-03) """ R = self.ring() if R is None: return None # Coerce keys of p into R. p = FFPD.coerce_point(R, p) if asy_var is None: asy_var = var('r') d = self.dimension() X = list(R.gens()) alpha = list(alpha) df = self.denominator_factored() # Find greatest i such that X[i] is a convenient coordinate, # that is, such that for all (h, e) in df, we have # (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]: i -= 1 coordinate = i # Decompose self into a sum of simpler pieces. # Where each piece has a denominator with at most d distinct # factors and no repeated factors. decomp = self.asymptotic_decomposition(alpha, asy_var) # Sum asymptotic expansions of terms of decomp. asy_expansions = [] for f in decomp: if f.quotient() == 0: continue critical_cone_p = f.critical_cone(p, coordinate) if [alpha[i] for i in xrange(d)] == [True for i in xrange(d)] and\ alpha not in critical_cone_p[1]: # f does not contribute to asymptotics of self # in direction alpha. continue n = len(f.denominator_factored()) if n == 1: # Smooth point. asy_expansions.append( f.asymptotics_smooth(p, alpha, N, asy_var, coordinate, numerical)) else: # Multiple point with n <= d. asy_expansions.append( f.asymptotics_multiple(p, alpha, N, asy_var, coordinate, numerical)) if asy_expansions: asy = sum([a[0] for a in asy_expansions]) exp_scale = asy_expansions[0][1] # Same for all. subexp_part = sum([a[2] for a in asy_expansions]).\ simplify_full() return (asy, exp_scale, subexp_part) else: return None def asymptotics_smooth(self, p, alpha, N, asy_var, coordinate=None, numerical=0): r""" Does what asymptotics() does but only in the case of a convenient smooth point and assuming that the denominator of self contains no repeated factors. INPUT: - p - A dictionary with keys that can be coerced to equal self.ring().gens(). - alpha - A tuple of length d = self.dimension() of positive integers or, if $p$ is a smooth point, possibly of symbolic variables. - N - A positive integer. - asy_var - (Optional; default=None) A symbolic variable. The variable of the asymptotic expansion. If none is given, var('r') will be assigned. - coordinate- (Optional; default=None) An integer in {0, ..., d-1} indicating a convenient coordinate to base the asymptotic calculations on. If None is assigned, then choose coordinate=d-1. - numerical - (Optional; default=0) A natural number. If numerical > 0, then return a numerical approximation of the Maclaurin ray coefficients of self with numerical digits of precision. Otherwise return exact values. NOTES: The formulas used for computing the asymptotic expansions are Theorems 3.2 and 3.3 [RaWi2008a]_ with the exponent of H equal to 1. Theorem 3.2 is a specialization of Theorem 3.4 of [RaWi2012]_ with $n=1$. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = 1-x-y-x*y sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = [3, 2] sage: p = {y: 1/2*sqrt(13) - 3/2, x: 1/3*sqrt(13) - 2/3} sage: print F.asymptotics_smooth(p, alpha, 2, var('r'), numerical=3) # long time Creating auxiliary functions... Computing derivatives of auxiallary functions... Computing derivatives of more auxiliary functions... Computing second order differential operator actions... ((0.369/sqrt(r) - 0.0186/r^(3/2))*71.2^r, 71.2, 0.369/sqrt(r) - 0.0186/r^(3/2)) :: sage: R. = PolynomialRing(QQ) sage: q = 1/2 sage: qq = q.denominator() sage: H = 1 - q*x + q*x*y - x^2*y sage: Hfac = H.factor() sage: G = (1 - q*x)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = list(qq*vector([2, 1 - q])) sage: print alpha [4, 1] sage: p = {x: 1, y: 1} sage: print F.asymptotics_smooth(p, alpha, 5, var('r')) # long time Creating auxiliary functions... Computing derivatives of auxiallary functions... Computing derivatives of more auxiliary functions... Computing second order differential operator actions... (1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*r^(5/3)), 1, 1/12*2^(2/3)*sqrt(3)*gamma(1/3)/(pi*r^(1/3)) - 1/96*2^(1/3)*sqrt(3)*gamma(2/3)/(pi*r^(5/3))) AUTHORS: - Alexander Raichev (2008-10-01, 2010-09-28, 2012-08-01) """ R = self.ring() if R is None: return None d = self.dimension() I = sqrt(-Integer(1)) # Coerce everything into the Symbolic Ring. X = [SR(x) for x in R.gens()] G = SR(self.numerator()) H = SR(self.denominator()) p = dict([(SR(x), p[x]) for x in R.gens()]) alpha = [SR(a) for a in alpha] # Put given convenient coordinate at end of variable list. if coordinate is not None: x = X.pop(coordinate) X.append(x) a = alpha.pop(coordinate) alpha.append(a) # 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: 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)] ) # Setup. print "Creating auxiliary functions..." # Implicit functions. h = function('h', *tuple(X[:d - 1])) 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 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)]) ht = h.subs(e) Ht = H.subs(e).subs({X[d - 1]: Integer(1)/ht}) 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]) # 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]]}) # 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]), 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 = FFPD.diff_all(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 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 = FFPD.diff_all(H, X, 2*N, ending=[X[d - 1]], sub_final=P) 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 ={} atP.update({U.subs(P): diff(H, X[d - 1]).subs(P)}) end = [X[d - 1]] 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 s in UnorderedTuples(X, l): Uderivs[diff(U, s).subs(P)] = Integer(0) elif k > 0 and k < 2*N: all_zero = True Uderivs =  FFPD.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))]: 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 s in UnorderedTuples(X, l): Uderivs.update({diff(U, s).subs(P): Integer(0)}) else: # Have to compute the rest of the derivatives. Uderivs = FFPD.diff_prod(Hderivs, U, Hcheck, X, range(k + 1, 2*N + 1), end, Uderivs, atP) else: Uderivs = FFPD.diff_prod(Hderivs, U, Hcheck, X, range(1, 2*N + 1), end, Uderivs, atP) atP.update(Uderivs) # In general, this algorithm is not designed to handle the case of a # singular Phit''(Tstar). # However, when d = 2 the algorithm can cope. if d == 2: # Compute v, the order of vanishing at Tstar of Phit. # It is at least 2. v = Integer(2) Phitderiv = diff(Phit, T[0], 2) splat = Phitderiv.subs(Tstar).subs(atP).subs(p).simplify() while splat == 0: v += 1 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)], X[0]).subs(hderivs1)}) atP.update({diff(h, X[0], v).subs(P): hderivs[diff(h, X[0], v)].subs(atP)}) Phitderiv = diff(Phitderiv, T[0]) splat = Phitderiv.subs(Tstar).subs(atP).subs(p).simplify() if d == 2 and v > 2: t = T[0]  # Simplify variable names. 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, # in case v is even. # Otherwise, compute up to orders N - 1 and N - 1 + v, # respectively. # To speed up later computations, # create symbolic functions AA and BB # to stand in for the expressions At and Phitu, respectively. print "Computing derivatives of more auxiliary functions..." AA = function('AA', t) BB = function('BB', t) if v.mod(2) == 0: At_derivs = FFPD.diff_all(At, T, 2*N - 2, sub=hderivs1, sub_final=[Tstar, atP], rekey=AA) Phitu_derivs = FFPD.diff_all(Phitu, T, 2*N - 2 +v, sub=hderivs1, sub_final=[Tstar, atP], zero_order=v + 1, rekey=BB) else: At_derivs = FFPD.diff_all(At, T, N - 1, sub=hderivs1, sub_final=[Tstar, atP], rekey=AA) Phitu_derivs = FFPD.diff_all(Phitu, T, N - 1 + v, sub=hderivs1, sub_final=[Tstar, atP], zero_order=v + 1 , rekey=BB) AABB_derivs = At_derivs AABB_derivs.update(Phitu_derivs) AABB_derivs[AA] = At.subs(Tstar).subs(atP) AABB_derivs[BB] = Phitu.subs(Tstar).subs(atP) print "Computing second order differential operator actions..." DD = FFPD.diff_op_simple(AA, BB, AABB_derivs, t, v, a, N) # Plug above into asymptotic formula. 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) ]) else: 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) ]) # Asymptotics for d >= 2 case. # A singular Phit''(Tstar) will cause a crash in this case. else: 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 = 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 # are sufficiently differentiable functions so that mixed partials # are equal.  Thus only need to compute representative partials. # Choose nondecreasing sequences as representative differentiation- # order sequences. # To speed up later computations, # create symbolic functions AA and BB # to stand in for the expressions At and Phitu, respectively. print "Computing derivatives of more auxiliary functions..." AA = function('AA', *tuple(T)) At_derivs = FFPD.diff_all(At, T, 2*N - 2, sub=hderivs1, sub_final =[Tstar, atP], rekey=AA) BB = function('BB', *tuple(T)) Phitu_derivs = FFPD.diff_all(Phitu, T, 2*N, sub=hderivs1, sub_final =[Tstar, atP], rekey=BB, zero_order=3) AABB_derivs = At_derivs AABB_derivs.update(Phitu_derivs) AABB_derivs[AA] = At.subs(Tstar).subs(atP) AABB_derivs[BB] = Phitu.subs(Tstar).subs(atP) print "Computing second order differential operator actions..." DD = FFPD.diff_op(AA, BB, AABB_derivs, T, a_inv, 1 , N) # Plug above into asymptotic formula. 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) ]) 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) 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) def asymptotics_multiple(self, p, alpha, N, asy_var, coordinate=None, numerical=0): r""" Does what asymptotics() does but only in the case of a convenient multiple point nondegenerate for alpha and assuming that that the number of distinct irreducible factors of the denominator of self is at most self.dimension() and that no factors are repeated. Assume also that p.values() are not symbolic variables. INPUT: - p - A dictionary with keys that can be coerced to equal self.ring().gens(). - alpha - A tuple of length d = self.dimension() of positive integers or, if $p$ is a smooth point, possibly of symbolic variables. - N - A positive integer. - asy_var - (Optional; default=None) A symbolic variable. The variable of the asymptotic expansion. If none is given, var('r') will be assigned. - coordinate- (Optional; default=None) An integer in {0, ..., d-1} indicating a convenient coordinate to base the asymptotic calculations on. If None is assigned, then choose coordinate=d-1. - numerical - (Optional; default=0) A natural number. If numerical > 0, then return a numerical approximation of the Maclaurin ray coefficients of self with numerical digits of precision. Otherwise return exact values. NOTES: The formulas used for computing the asymptotic expansion are Theorem 3.4 and Theorem 3.7 of [RaWi2012]_. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (4 - 2*x - y - z)*(4 - x -2*y - z) sage: Hfac = H.factor() sage: G = 16/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F (16, [(x + 2*y + z - 4, 1), (2*x + y + z - 4, 1)]) sage: p = {x: 1, y: 1, z: 1} sage: alpha = [3, 3, 2] sage: print F.asymptotics_multiple(p, alpha, 2, var('r')) # long time Creating auxiliary functions... Computing derivatives of auxiliary functions... Computing derivatives of more auxiliary functions... Computing second-order differential operator actions... (4/3*sqrt(3)/(sqrt(pi)*sqrt(r)) - 25/216*sqrt(3)/(sqrt(pi)*r^(3/2)), 1, 4/3*sqrt(3)/(sqrt(pi)*sqrt(r)) - 25/216*sqrt(3)/(sqrt(pi)*r^(3/2))) :: sage: R.= PolynomialRing(QQ) 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) sage: print F (1, [(x*y + x - 1, 1), (2*x^2*y*z + x^2*z - 1, 1)]) sage: p = {x: 1/2, z: 4/3, y: 1} sage: alpha = [8, 3, 3] sage: print F.asymptotics_multiple(p, alpha, 2, var('r'), coordinate=1) # long time Creating auxiliary functions... Computing derivatives of auxiliary functions... Computing derivatives of more auxiliary functions... Computing second-order differential operator actions... (1/172872*(24696*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(r)) - 1231*sqrt(3)*sqrt(7)/(sqrt(pi)*r^(3/2)))*108^r, 108, 1/7*sqrt(3)*sqrt(7)/(sqrt(pi)*sqrt(r)) - 1231/172872*sqrt(3)*sqrt(7)/(sqrt(pi)*r^(3/2))) :: sage: R.= PolynomialRing(QQ) sage: H = (1 - 2*x - y) * (1 - x - 2*y) sage: Hfac = H.factor() sage: G = exp(x + y)/Hfac.unit() sage: F = FFPD(G, Hfac) sage: print F (e^(x + y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)]) sage: p = {x: 1/3, y: 1/3} sage: alpha = (var('a'), var('b')) sage: print F.asymptotics_multiple(p, alpha, 2, var('r')) # long time (3*((1/3)^(-b)*(1/3)^(-a))^r*e^(2/3), (1/3)^(-b)*(1/3)^(-a), 3*e^(2/3)) AUTHORS: - Alexander Raichev (2008-10-01, 2010-09-28, 2012-08-02) """ from itertools import product R = self.ring() if R is None: return None # Coerce keys of p into R. p = FFPD.coerce_point(R, p) d = self.dimension() I = sqrt(-Integer(1)) # Coerce everything into the Symbolic Ring. X = [SR(x) for x in R.gens()] G = SR(self.numerator()) 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()]) Sstar = self.crit_cone_combo(p, alpha, coordinate) # Put the given convenient variable at end of variable list. if coordinate is not None: x = X.pop(coordinate) X.append(x) 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())) 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) # 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: 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)] ) # Setup. print "Creating auxiliary functions..." # Create T and S variables. 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)] s = 's' while s in [str(x) for x in X]: s = s + 't' S = [var(s + str(i)) for i in xrange(n - 1)] Sstar = dict([(S[j], Sstar[j]) for j in xrange(n - 1)]) thetastar = dict([(t, Integer(0)) for t in T]) thetastar.update(Sstar) # Create implicit functions. h = [function('h' + str(j), *tuple(X[:d - 1])) for j in xrange(n)] 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) for j in xrange(n)] 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] Ht = [H[j].subs(e).subs({X[d - 1]: Integer(1)/ht[j]}) for j in xrange(n)] hsumt = add([S[j]*ht[j] for j in xrange(n - 1)]) +\ (Integer(1) - add(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)]) # 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])) # 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]), 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 = FFPD.diff_all(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 # 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. print "Computing derivatives of auxiliary functions..." m = min(n, N) end = [X[d-1] for j in xrange(n)] Hprodderivs = FFPD.diff_all(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 ={} 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 s in UnorderedTuples(X, l): Uderivs[diff(U, s).subs(P)] = Integer(0) elif k > 0 and k < 2*N - 2 + m - 1: all_zero = True Uderivs = FFPD.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))]: 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 s in UnorderedTuples(X, l): Uderivs.update({diff(U, s).subs(P): Integer(0)}) else: # Have to compute the rest of the derivatives. Uderivs = FFPD.diff_prod(Hprodderivs, U, Hcheck, X, range(k + 1, 2*N - 2 + m), end, Uderivs, atP) else: Uderivs = FFPD.diff_prod(Hprodderivs, U, Hcheck, X, 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 = 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 # are sufficiently differentiable functions so that mixed partials # are equal.  Thus only need to compute representative partials. # Choose nondecreasing sequences as representative differentiation- # order sequences. # To speed up later computations, create symbolic functions AA and BB # to stand in for the expressions At and Phitu respectively. print "Computing derivatives of more auxiliary functions..." AA = [function('A' + str(j), *tuple(T + S)) for j in xrange(n)] At_derivs = FFPD.diff_all(At, T + S, 2*N - 2, sub=hderivs1, sub_final =[thetastar, atP], rekey=AA) BB = function('BB', *tuple(T + S)) Phitu_derivs = FFPD.diff_all(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) for j in xrange(n): AABB_derivs[AA[j]] = At[j].subs(thetastar).subs(atP) AABB_derivs[BB] = Phitu.subs(thetastar).subs(atP) print "Computing second-order differential operator actions..." DD = FFPD.diff_op(AA, BB, AABB_derivs, T + S, a_inv, n, N) 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)]) 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) 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) @staticmethod def subs_all(f, sub, simplify=False): r""" Return the items of $f$ substituted by the dictionaries of $sub$ in order of their appearance in $sub$. INPUT: - f - An individual or list of symbolic expressions or dictionaries - sub - An individual or list of dictionaries. - simplify - Boolean (default: False). OUTPUT: The items of $f$ substituted by the dictionaries of $sub$ in order of their appearance in $sub$. The subs() command is used. If simplify is True, then simplify() is used after substitution. EXAMPLES:: sage: from sage.combinat.amgf import * sage: var('x, y, z') (x, y, z) sage: a = {x:1} sage: b = {y:2} sage: c = {z:3} sage: FFPD.subs_all(x + y + z, a) y + z + 1 sage: FFPD.subs_all(x + y + z, [c, a]) y + 4 sage: FFPD.subs_all([x + y + z, y^2], b) [x + z + 2, 4] sage: FFPD.subs_all([x + y + z, y^2], [b, c]) [x + 5, 4] :: sage: var('x, y') (x, y) sage: a = {'foo': x**2 + y**2, 'bar': x - y} sage: b = {x: 1 , y: 2} sage: FFPD.subs_all(a, b) {'foo': 5, 'bar': -1} AUTHORS: - Alexander Raichev (2009-05-05) """ singleton = False if not isinstance(f, list): f = [f] singleton = True if not isinstance(sub, list): sub = [sub] g = [] for ff in f: for D in sub: if isinstance(ff, dict): 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) ] else: return g[Integer(0) ].simplify() elif singleton and not simplify: return g[Integer(0) ] elif not singleton and simplify: G = [] for gg in g: if isinstance(gg, dict): G.append(gg) else: G.append(gg.simplify()) return G else: return g @staticmethod def diff_all(f, V, n, ending=[], sub=None, sub_final=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 variables in $V$. The default is to key the dictionary by all nondecreasing sequences in $V$ of length 1 up to length $n$. For internal use. INPUT: - f - An individual or list of $\mathcal{C}^{n+1}$ functions. - V - A list of variables occurring in $f$. - n - A natural number. - ending - A list of variables in $V$. - sub - An individual or list of dictionaries. - sub_final - An individual or list of dictionaries. - rekey - A callable symbolic function in $V$ or list thereof. - zero_order - A natural number. OUTPUT: The dictionary ${s_1:deriv_1,..., s_r:deriv_r}$. Here $s_1,\ldots, s_r$ is a listing of all nondecreasing sequences of length 1 up to length $n$ over the alphabet $V$, where $w > v$ in $X$ iff $str(w) > str(v)$, and $deriv_j$ is the derivative of $f$ with respect to the derivative sequence $s_j$ and simplified with respect to the substitutions in $sub$ and evaluated at $sub_final$. Moreover, all derivatives with respect to sequences of length less than $zero_order$ (derivatives of order less than $zero_order$ ) will be made zero. If $rekey$ is nonempty, then $s_1,\ldots, s_r$ will be replaced by the symbolic derivatives of the functions in $rekey$. If $ending$ is nonempty, then every derivative sequence $s_j$ will be suffixed by $ending$. EXAMPLES:: I'd like to print the entire dictionaries, but that doesn't yield consistent output order for doctesting. Order of keys changes.:: sage: from sage.combinat.amgf import * sage: f = function('f', x) sage: dd = FFPD.diff_all(f, [x], 3) sage: dd[(x, x, x)] D[0, 0, 0](f)(x) :: sage: d1 = {diff(f, x): 4*x^3} sage: dd = FFPD.diff_all(f,[x], 3, sub=d1) sage: dd[(x, x, x)] 24*x :: sage: dd = FFPD.diff_all(f,[x], 3, sub=d1, rekey=f) sage: dd[diff(f, x, 3)] 24*x :: sage: a = {x:1} sage: dd = FFPD.diff_all(f,[x], 3, sub=d1, rekey=f, sub_final=a) sage: dd[diff(f, x, 3)] 24 :: sage: X = var('x, y, z') sage: f = function('f',*X) sage: dd = FFPD.diff_all(f, X, 2, ending=[y, y, y]) sage: dd[(z, y, y, y)] D[1, 1, 1, 2](f)(x, y, z) :: sage: g = function('g',*X) sage: dd = FFPD.diff_all([f, g], X, 2) sage: dd[(0, y, z)] D[1, 2](f)(x, y, z) :: sage: dd[(1, z, z)] D[2, 2](g)(x, y, z) :: sage: f = exp(x*y*z) sage: ff = function('ff',*X) sage: dd = FFPD.diff_all(f, X, 2, rekey=ff) sage: dd[diff(ff, x, z)] x*y^2*z*e^(x*y*z) + y*e^(x*y*z) AUTHORS: - Alexander Raichev (2008-10-01, 2009-04-01, 2010-02-01) """ singleton=False if not isinstance(f, list): f = [f] singleton=True # Build the dictionary of derivatives iteratively from a list # of nondecreasing derivative-order sequences. derivs = {} r = len(f) if ending: seeds = [ending] start = Integer(1) else: seeds = [[v] for v in V] start = Integer(2) if singleton: for s in seeds: derivs[tuple(s)] = FFPD.subs_all(diff(f[0], s), sub) for l in xrange(start, n + 1): for t in UnorderedTuples(V, l): s = tuple(t + ending) derivs[s] = FFPD.subs_all(diff(derivs[s[1:]], s[0]), sub) else: # Make the dictionary keys of the form (j, sequence of variables), # where j in range(r). for s in seeds: value = FFPD.subs_all([diff(f[j], s) for j in xrange(r)], sub) 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 = FFPD.subs_all([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: # Zero out all the derivatives of order < zero_order if singleton: for k in derivs.keys(): if len(k) < zero_order: derivs[k] = Integer(0) else: # Ignore the first of element of k, which is an index. for k in derivs.keys(): if len(k) - 1  < zero_order: derivs[k] = Integer(0) if sub_final: # Substitute sub_final into the values of derivs. for k in derivs.keys(): derivs[k] = FFPD.subs_all(derivs[k], sub_final) if rekey: # Rekey the derivs dictionary by the value of rekey. F = rekey if singleton: # F must be a singleton. 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()] ) return derivs @staticmethod def diff_op(A, B, AB_derivs, V, M, r, N): 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 of $A$ and $B$ evaluated at $p$ that are necessary for the computation. For internal use by the functions asymptotics_smooth() and asymptotics_multiple(). INPUT: - A - A single or length r list of symbolic functions in the 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. 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. - M - A symmetric $l \times l$ matrix, where $l$ is the length of V. - r, N - Natural numbers. OUTPUT: A dictionary whose keys are natural number tuples of the form $(j, k, l)$, where $l \le 2k$, $j \le r-1$, and $j+k \le N-1$, and whose values are $DD^(l+k)(A[j] B^l)$ evaluated at a point $p$, where $DD$ is the linear second-order differential operator $-\sum_{i=0}^{l-1} \sum_{j=0}^{l-1} M[i][j] \partial^2 /(\partial V[j] \partial V[i])$. EXAMPLES:: sage: from sage.combinat.amgf import * sage: T = var('x, y') sage: A = function('A',*tuple(T)) sage: B = function('B',*tuple(T)) sage: AB_derivs = {} sage: M = matrix([[1, 2],[2, 1]]) sage: DD = FFPD.diff_op(A, B, AB_derivs, T, M, 1, 2) sage: DD.keys() [(0, 1, 2), (0, 1, 1), (0, 1, 0), (0, 0, 0)] sage: len(DD[(0, 1, 2)]) 246 AUTHORS: - Alexander Raichev (2008-10-01, 2010-01-12) """ if not isinstance(A, list): A = [A] # First, compute the necessary product derivatives of A and B. 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) # 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): # 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 : 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)]) diffo = Integer(0) for t in P: if product_derivs[(j, k, l) + FFPD.diff_seq(V, t)] !=\ Integer(0): MM = Integer(1) for (a, b) in t: MM = MM * M[a][b] if a != b: MM = Integer(2) *MM diffo = diffo + MM*product_derivs[(j, k, l) +\ FFPD.diff_seq(V, t)] DD[(j, k, l)] = (-Integer(1) )**(k+l)*diffo return DD @staticmethod def diff_seq(V, s): r""" Given a list s of tuples of natural numbers, return the list of elements of V with indices the elements of the elements of s. This function is for internal use by the function diff_op(). INPUT: - V - A list. - s - A list of tuples of natural numbers in the interval range(len(V)). OUTPUT: The tuple tuple([V[tt] for tt in sorted(t)]), where t is the list of elements of the elements of s. EXAMPLES:: sage: from sage.combinat.amgf import * sage: V = list(var('x, t, z')) sage: FFPD.diff_seq(V,([0, 1],[0, 2, 1],[0, 0])) (x, x, x, x, t, t, z) AUTHORS: - Alexander Raichev (2009-05-19) """ t = [] for ss in s: t.extend(ss) return tuple([V[tt] for tt in sorted(t)]) @staticmethod def diff_op_simple(A, B, AB_derivs, x, v, a, N): r""" Return $DD^(e k + v l)(A B^l)$ evaluated at a point $p$ for various natural numbers $e, k, l$ that depend on $v$ and $N$. Here $DD$ is a specific linear differential operator that depends on $a$ and $v$ , $A$ and $B$ are symbolic functions, and $AB_derivs$ contains all the derivatives of $A$ and $B$ evaluated at $p$ that are necessary for the computation. For internal use by the function asymptotics_smooth(). INPUT: - 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 N if v is odd and the (symbolic) derivatives of B 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 - Symbolic variable. - a - A complex number. - v, N - Natural numbers. OUTPUT: A dictionary whose keys are natural number pairs of the form $(k, l)$, where $k < N$ and $l \le 2k$ and whose values are $DD^(e k + v l)(A B^l)$ evaluated at a point $p$. Here $e=2$ if $v$ is even, $e=1$ if $v$ is odd, and $DD$ is the linear differential operator $(a^{-1/v} d/dt)$ if $v$ is even and $(|a|^{-1/v} i \text{sgn}(a) d/dt)$ if $v$ is odd. EXAMPLES:: sage: from sage.combinat.amgf import * sage: A = function('A', x) sage: B = function('B', x) sage: AB_derivs = {} sage: FFPD.diff_op_simple(A, B, AB_derivs, x, 3, 2, 2) {(1, 0): 1/2*I*2^(2/3)*D[0](A)(x), (0, 0): A(x), (1, 1): 1/4*(A(x)*D[0, 0, 0, 0](B)(x) + B(x)*D[0, 0, 0, 0](A)(x) + 4*D[0](A)(x)*D[0, 0, 0](B)(x) + 4*D[0](B)(x)*D[0, 0, 0](A)(x) + 6*D[0, 0](A)(x)*D[0, 0](B)(x))*2^(2/3)} AUTHORS: - Alexander Raichev (2010-01-15) """ I = sqrt(-Integer(1)) DD = {} 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) 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) return DD @staticmethod def diff_prod(f_derivs, u, g, X, interval, end, uderivs, atc): r""" Take various derivatives of the equation $f = ug$, evaluate them at a point $c$, and solve for the derivatives of $u$. For internal use by the function asymptotics_multiple(). INPUT: - f_derivs - A dictionary whose keys are all tuples of the form s + end, where s is a sequence of variables from X whose length lies in interval, and whose values are the derivatives of a function $f$ evaluated at $c$. - u - A callable symbolic function. - g - An expression or callable symbolic function. - X - A list of symbolic variables. - interval - A list of positive integers. Call the first and last values $n$ and $nn$, respectively. - end - A possibly empty list of repetitions of the variable z, where z is the last element of X. - uderivs - A dictionary whose keys are the symbolic derivatives of order 0 to order $n-1$ of u evaluated at $c$ and whose values are the corresponding derivatives evaluated at $c$. - atc - A dictionary whose keys are the keys of $c$ and all the symbolic derivatives of order 0 to order $nn$ of g evaluated $c$ and whose values are the corresponding derivatives evaluated at $c$. OUTPUT: A dictionary whose keys are the derivatives of u up to order $nn$ and whose values are those derivatives evaluated at $c$. EXAMPLES:: I'd like to print out the entire dictionary, but that does not give consistent output for doctesting. Order of keys changes :: sage: from sage.combinat.amgf import * sage: u = function('u', x) sage: g = function('g', x) sage: fd = {(x,):1,(x, x):1} sage: ud = {u(x=2): 1} sage: atc = {x: 2, g(x=2): 3, diff(g, x)(x=2): 5} sage: atc[diff(g, x, x)(x=2)] = 7 sage: dd = FFPD.diff_prod(fd, u, g, [x], [1, 2], [], ud, atc) sage: dd[diff(u, x, 2)(x=2)] 22/9 NOTES: This function works by differentiating the equation $f = ug$ with respect to the variable sequence s + end, for all tuples s of X of lengths in interval, evaluating at the point $c$ , and solving for the remaining derivatives of u. This function assumes that u never appears in the differentiations of $f = ug$ after evaluating at $c$. AUTHORS: - Alexander Raichev (2009-05-14, 2010-01-21) """ for l in interval: D = {} rhs = [] lhs = [] 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)) # 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' +\ ''.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(FFPD.subs_all(D, sol[Integer(0) ])) return uderivs def crit_cone_combo(self, p, alpha, coordinate=None): r""" Return an auxiliary point associated to the multiple point p of the factors self. For internal use by asymptotics_multiple(). INPUT: - p - A dictionary with keys that can be coerced to equal self.ring().gens(). - alpha - A list of rationals. OUTPUT: A solution of the matrix equation $y \Gamma = \alpha'$ for $y$ , where $\Gamma$ is the matrix given by [FFPD.direction(v) for v in self.log_grads(p)] and $\alpha'$ is FFPD.direction(alpha) EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: p = exp(x) sage: df = [(1 - 2*x - y, 1), (1 - x - 2*y, 1)] sage: f = FFPD(p, df) sage: p = {x: 1/3, y: 1/3} sage: alpha = (var('a'), var('b')) sage: print f.crit_cone_combo(p, alpha) [1/3*(2*a - b)/b, -2/3*(a - 2*b)/b] NOTES: Use this function only when $\Gamma$ is well-defined and there is a unique solution to the matrix equation $y \Gamma = \alpha'$. Fails otherwise. AUTHORS: - Alexander Raichev (2008-10-01, 2008-11-25, 2009-03-04, 2010-09-08, 2010-12-02, 2012-08-02) """ # Assuming here that each log_grads(f) has nonzero final component. # Then 'direction' will not throw a division by zero error. R = self.ring() if R is None: return None # Coerce keys of p into R. p = FFPD.coerce_point(R, p) d = self.dimension() n = len(self.denominator_factored()) Gamma = matrix([FFPD.direction(v, coordinate) for v in self.log_grads(p)]) beta = FFPD.direction(alpha, coordinate) # 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 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] @staticmethod def direction(v, coordinate=None): r""" Returns [vv/v[coordinate] for vv in v] where coordinate is the last index of v if not specified otherwise. INPUT: - v - A vector. - coordinate - (Optional; default=None) An index for v. EXAMPLES:: sage: from sage.combinat.amgf import * sage: FFPD.direction([2, 3, 5]) (2/5, 3/5, 1) sage: FFPD.direction([2, 3, 5], 0) (1, 3/2, 5/2) AUTHORS: - Alexander Raichev (2010-08-25) """ if coordinate is None: coordinate = len(v) - 1 return tuple([vv/v[coordinate] for vv in v]) def grads(self, p): r""" Return a list of the gradients of the polynomials [q for (q, e) in self.denominator_factored()] evalutated at p. INPUT: - p - (Optional: default=None) A dictionary whose keys are the generators of self.ring(). EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: p = exp(x) sage: df = [(x**3 + 3*y^2, 5), (x*y, 2), (y, 1)] sage: f = FFPD(p, df) sage: print f (e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)]) sage: print R.gens() (x, y) sage: p = None sage: print f.grads(p) [(0, 1), (y, x), (3*x^2, 6*y)] :: sage: p = {x: sqrt(2), y: var('a')} sage: print f.grads(p) [(0, 1), (a, sqrt(2)), (6, 6*a)] AUTHORS: - Alexander Raichev (2009-03-06) """ R = self.ring() if R is None: return # Coerce keys of p into R. p = FFPD.coerce_point(R, p) X = R.gens() d = self.dimension() H = [h for (h, e) in self.denominator_factored()] n = len(H) return [tuple([diff(H[i], X[j]).subs(p) for j in xrange(d)]) for i in xrange(n)] def log_grads(self, p): r""" Return a list of the logarithmic gradients of the polynomials [q for (q, e) in self.denominator_factored()] evalutated at p. INPUT: - p - (Optional: default=None) A dictionary whose keys are the generators of self.ring(). NOTE: The logarithmic gradient of a function $f$ at point $p$ is the vector $(x_1 \partial_1 f(x), \ldots, x_d \partial_d f(x) )$ evaluated at $p$. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: p = exp(x) sage: df = [(x**3 + 3*y^2, 5), (x*y, 2), (y, 1)] sage: f = FFPD(p, df) sage: print f (e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)]) sage: print R.gens() (x, y) sage: p = None sage: print f.log_grads(p) [(0, y), (x*y, x*y), (3*x^3, 6*y^2)] :: sage: p = {x: sqrt(2), y: var('a')} sage: print f.log_grads(p) [(0, a), (sqrt(2)*a, sqrt(2)*a), (6*sqrt(2), 6*a^2)] AUTHORS: - Alexander Raichev (2009-03-06) """ R = self.ring() if R is None: return None # Coerce keys of p into R. p = FFPD.coerce_point(R, p) X = R.gens() 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)]) for i in xrange(n)] def critical_cone(self, p, coordinate=None): r""" Return the critical cone of the convenient multiple point p. INPUT: - p - A dictionary with keys that can be coerced to equal self.ring().gens() and values in a field. - coordinate - (Optional; default=None) A natural number. OUTPUT: A list of vectors that generate the critical cone of p and the cone itself, which is None if the values of p don't lie in QQ. Divide logarithmic gradients by their component coordinate entries. If coordinate=None, then search from d-1 down to 0 for the first index j such that for all i self.log_grads()[i][j] != 0 and set coordinate=j. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: G = 1 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) sage: p = {x: 1/2, y: 1, z: 4/3} sage: print F.critical_cone(p) ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N) AUTHORS: - Alexander Raichev (2010-08-25, 2012-08-02) """ R = self.ring() if R is None: return # Coerce keys of p into R. p = FFPD.coerce_point(R, p) X = R.gens() d = self.dimension() lg = self.log_grads(p) n = len(lg) if coordinate not in xrange(d): # Search from d-1 down to 0 for a coordinate j such that # for all i we have lg[i][j] != 0. # One is guaranteed to exist in the case of a convenient multiple # point. for j in reversed(xrange(d)): if 0 not in [lg[i][j] for i in xrange(n)]: coordinate = j break Gamma = [FFPD.direction(v, coordinate) for v in lg] try: cone = Cone(Gamma) except TypeError: cone = None return (Gamma, cone) def is_convenient_multiple_point(self, p): r""" Return True if p is a convenient multiple point of self and False otherwise. Also return a short comment. INPUT: - p - A dictionary with keys that can be coerced to equal self.ring().gens(). OUTPUT: A pair (verdict, comment). In case p is a convenient multiple point, verdict=True and comment ='No problem'. In case p is not, verdict=False and comment is string explaining why p fails to be a convenient multiple point. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x*(1 + y))*(1 - z*x**2*(1 + 2*y)) sage: df = H.factor() 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} sage: print F.is_convenient_multiple_point(p1) (True, 'convenient in variables [x, y]') sage: print F.is_convenient_multiple_point(p2) (False, 'not a singular point') NOTES: See [RaWi2012]_ for more details. AUTHORS: - Alexander Raichev (2011-04-18, 2012-08-02) """ R = self.ring() if R is None: return # Coerce keys of p into R. p = FFPD.coerce_point(R, p) H = [h for (h, e) in self.denominator_factored()] n = len(H) d = self.dimension() # Test 1: Are the factors in H zero at p? if [h.subs(p) for h in H] != [0 for h in H]: # Failed test 1.  Move on to next point. return (False, 'not a singular point') # Test 2: Are the factors in H smooth at p? grads = self.grads(p) for v in grads: if v == [0 for i in xrange(d)]: return (False, 'not smooth point of factors') # Test 3: Do the factors in H intersect transversely at p? if n <= d: M = matrix(grads) if M.rank() != n: return (False, 'not a transverse intersection') else: # Check all sub-multisets of grads of size d. for S in Subsets(grads, d, submultiset=True): M = matrix(S) if M.rank() != d: return (False, 'not a transverse intersection') # Test 4: Is p convenient? M = matrix(self.log_grads(p)) convenient_coordinates = [] for j in xrange(d): if 0 not in M.columns()[j]: convenient_coordinates.append(j) if not convenient_coordinates: return (False, 'multiple point but not convenient') # Tests all passed X = R.gens() return (True, 'convenient in variables %s' %\ [X[i] for i in convenient_coordinates]) def singular_ideal(self): r""" Let $R$ be the ring of self and $H$ its denominator. Let $Hred$ be the reduction (square-free part) of $H$. Return the ideal in $R$ generated by $Hred$ and its partial derivatives. If the coefficient field of $R$ is algebraically closed, then the output is the ideal of the singular locus (which is a variety) of the variety of $H$. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = (1 - x*(1 + y))**3*(1 - z*x**2*(1 + 2*y)) sage: df = H.factor() 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) of Multivariate Polynomial Ring in x, y, z over Rational Field AUTHORS: - Alexander Raichev (2008-10-01, 2008-11-20, 2010-12-03, 2011-04-18, 2012-08-03) """ R = self.ring() if R is None: return Hred = prod([h for (h, e) in self.denominator_factored()]) J = R.ideal([Hred] + Hred.gradient()) return R.ideal(J.groebner_basis()) def smooth_critical_ideal(self, alpha): r""" Let $R$ be the ring of self and $H$ its denominator. Return the ideal in $R$ of smooth critical points of the variety of $H$ for the direction alpha. If the variety $V$ of $H$ has no smooth points, then return the ideal in $R$ of $V$. INPUT: - alpha - A d-tuple of positive integers and/or symbolic entries, where d = self.ring().ngens(). EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: H = (1-x-y-x*y)^2 sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = var('a1, a2') sage: F.smooth_critical_ideal(alpha) Ideal (y^2 + ((-2*a1)/(-a2))*y - 1, x + (a2/(-a1))*y + (a1 - a2)/(-a1)) of Multivariate Polynomial Ring in x, y over Fraction Field of Multivariate Polynomial Ring in a1, a2 over Rational Field :: sage: R. = PolynomialRing(QQ) sage: H = (1-x-y-x*y)^2 sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = [7/3, var('a')] sage: F.smooth_critical_ideal(alpha) Ideal (y^2 + (-14/(-3*a))*y - 1, x + (-3/7*a)*y + 3/7*a - 1) of Multivariate Polynomial Ring in x, y over Fraction Field of Univariate Polynomial Ring in a over Rational Field NOTES: See [RaWi2012]_ for more details. AUTHORS: - Alexander Raichev (2008-10-01, 2008-11-20, 2009-03-09, 2010-12-02, 2011-04-18, 2012-08-03) """ R = self.ring() if R is None: return Hred = prod([h for (h, e) in self.denominator_factored()]) K = R.base_ring() d = self.dimension() # Expand K by the variables of alpha if there are any. indets = [] for a in alpha: if a not in K and a in SR: indets.append(a) indets = list(Set(indets))   # Delete duplicates in indets. if indets: L = FractionField(PolynomialRing(K, indets)) S = R.change_ring(L) # Coerce alpha into L. alpha = [L(a) for a in alpha] else: S = R # 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)]) return S.ideal(J.groebner_basis()) def maclaurin_coefficients(self, multi_indices, numerical=0): r""" Returns the Maclaurin coefficients of self that have multi-indices alpha, 2*alpha, ..., r*alpha. INPUT: - multi_indices - A list of tuples of positive integers, where each tuple has length self.ring().ngens(). - numerical - (Optional; default=0) A natural number. If positive, return numerical approximations of coefficients with numerical digits of accuracy. OUTPUT: A dictionary of the form (nu, Maclaurin coefficient of index nu of self). EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: H = (4 - 2*x - y - z) * (4 - x - 2*y - z) sage: Hfac = H.factor() sage: G = 16/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = vector([3, 3, 2]) sage: interval = [1, 2, 4] sage: S = [r*alpha for r in interval] sage: print F.maclaurin_coefficients(S, numerical=10) {(6, 6, 4): 0.7005249476, (12, 12, 8): 0.5847732654, (3, 3, 2): 0.7849731445} NOTES: Uses iterated univariate Maclaurin expansions. Slow. AUTHORS: - Alexander Raichev (2011-04-08, 2012-08-03) """ R = self.ring() if R is None: return d = self.dimension() # Create biggest multi-index needed. alpha = [] for i in xrange(d): alpha.append(max((nu[i] for nu in multi_indices))) # Compute Maclaurin expansion of self up to index alpha. # Use iterated univariate expansions. # Slow! 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 = R(f) # Collect coefficients. coeffs = {} X = R.gens() for nu in multi_indices: monomial = prod([X[i]**nu[i] for i in xrange(d)]) coeffs[tuple(nu)] = F.monomial_coefficient(monomial) if numerical: coeffs[tuple(nu)] = coeffs[tuple(nu)].n(digits=numerical) return coeffs def relative_error(self, approx, alpha, interval, exp_scale=Integer(1), digits=10): r""" Returns the relative error between the values of the Maclaurin coefficients of self with multi-indices r alpha for r in interval and the values of the functions (of the variable r) in approx. INPUT: - approx - An individual or list of symbolic expressions in one variable. - alpha - A list of positive integers of length self.ring().ngens() - interval - A list of positive integers. - exp_scale - (Optional; default=1) A number. OUTPUT: A list whose entries are of the form [r*alpha, a_r, b_r, err_r] for r in interval. Here r*alpha is a tuple; a_r is the r*alpha (multi-index) coefficient of the Maclaurin series for self divided by exp_scale**r; b_r is a list of the values of the functions in approx evaluated at r and divided by exp_scale**m; err_r is the list of relative errors (a_r - f)/a_r for f in b_r. All outputs are decimal approximations. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: H = 1 - x - y - x*y sage: Hfac = H.factor() sage: G = 1/Hfac.unit() sage: F = FFPD(G, Hfac) sage: alpha = [1, 1] sage: r = var('r') sage: a1 = (0.573/sqrt(r))*5.83^r sage: a2 = (0.573/sqrt(r) - 0.0674/r^(3/2))*5.83^r sage: es = 5.83 sage: F.relative_error([a1, a2], alpha, [1, 2, 4, 8], es) # long time Calculating errors table in the form exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors... [((1, 1), 0.5145797599, [0.5730000000, 0.5056000000], [-0.1135300000, 0.01745066667]), ((2, 2), 0.3824778089, [0.4051721856, 0.3813426871], [-0.05933514614, 0.002967810973]), ((4, 4), 0.2778630595, [0.2865000000, 0.2780750000], [-0.03108344267, -0.0007627515584]), ((8, 8), 0.1991088276, [0.2025860928, 0.1996074055], [-0.01746414394, -0.002504047242])] AUTHORS: - Alexander Raichev (2009-05-18, 2011-04-18, 2012-08-03) """ if not isinstance(approx, list): approx = [approx] av = approx[0].variables()[0] print "Calculating errors table in the form" print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..." # Get Maclaurin coefficients of self. multi_indices = [r*vector(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) 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]]) else: stats_row.append([(mac[beta] - a)/mac[beta] for a in mac_approx[beta]]) stats.append(tuple(stats_row)) return stats @staticmethod def coerce_point(R, p): r""" Coerce the keys of the dictionary p into the ring R. Assume that it is possible. EXAMPLES:: AUTHORS: - Alexander Raichev (2009-05-18, 2011-04-18, 2012-08-03) """ result = p if p is not None and p.keys() and p.keys()[0].parent() != R: try: result = dict([(x, p[SR(x)]) for x in R.gens()]) except TypeError: pass return result class FFPDSum(list): r""" A list representing the sum of FFPD objects with distinct denominator factorizations. AUTHORS: - Alexander Raichev (2012-06-25) """ def __str__(self): r""" Returns a string representation of self EXAMPLES:: """ return str([r.list() for r in self]) def __eq__(self, other): r""" Returns True if self is equal to other EXAMPLES:: """ return sorted(self) == sorted(other) def __ne__(self, other): r""" Returns True if self is not equal to other EXAMPLES:: """ return not (self == other) def ring(self): r""" Return the polynomial ring of the denominators of self. If self does not have any denominators, then return None. EXAMPLES:: """ for r in self: R = r.ring() if R is not None: return R return None def whole_and_parts(self): r""" Rewrite self as a FFPDSum of a (possibly zero) polynomial FFPD followed by reduced rational expression FFPDs. Only useful for multivariate decompositions. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ, 'x, y') sage: f = x**2 + 3*y + 1/x + 1/y sage: f = FFPD(quotient=f) sage: print f (x^3*y + 3*x*y^2 + x + y, [(y, 1), (x, 1)]) sage: print FFPDSum([f]).whole_and_parts() [(x^2 + 3*y, []), (x + y, [(y, 1), (x, 1)])] :: sage: R. = PolynomialRing(QQ) sage: f = cos(x)**2 + 3*y + 1/x + 1/y sage: print f 1/x + 1/y + cos(x)^2 + 3*y sage: G = f.numerator() sage: H = R(f.denominator()) sage: f = FFPD(G, H.factor()) sage: print f (x*y*cos(x)^2 + 3*x*y^2 + x + y, [(y, 1), (x, 1)]) sage: print FFPDSum([f]).whole_and_parts() [(0, []), (x*y*cos(x)^2 + 3*x*y^2 + x + y, [(y, 1), (x, 1)])] """ whole = 0 parts = [] R = self.ring() for r in self: # Since r has already passed through FFPD.__init__()'s reducing # procedure, r is already in lowest terms. # Check if can write r as a mixed fraction: whole + fraction. p = r.numerator() q = r.denominator() if q == 1: # r is already whole whole += p else: try: # Coerce p into R and divide p by q p = R(p) a, b = p.quo_rem(q) except TypeError: # p is not in R and so can't divide p by q a = 0 b = p whole += a parts.append(FFPD(b, r.denominator_factored(), reduce_=False)) return FFPDSum([FFPD(whole, ())] + parts) def combine_like_terms(self): r""" Combine terms in self with the same denominator. Only useful for multivariate decompositions. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R.= PolynomialRing(QQ) sage: f = FFPD(quotient=1/(x * y * (x*y + 1))) sage: g = FFPD(quotient=x/(x * y * (x*y + 1))) sage: s = FFPDSum([f, g, f]) sage: t = s.combine_like_terms() sage: print s [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (1, [(y, 1), (x*y + 1, 1)]), (1, [(y, 1), (x, 1), (x*y + 1, 1)])] sage: print t [(1, [(y, 1), (x*y + 1, 1)]), (2, [(y, 1), (x, 1), (x*y + 1, 1)])] :: sage: R.= PolynomialRing(QQ) sage: H = x * y * (x*y + 1) sage: f = FFPD(1, H.factor()) sage: g = FFPD(exp(x + y), H.factor()) sage: s = FFPDSum([f, g]) sage: print s [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (e^(x + y), [(y, 1), (x, 1), (x*y + 1, 1)])] sage: t = s.combine_like_terms() sage: print t [(e^(x + y) + 1, [(y, 1), (x, 1), (x*y + 1, 1)])] """ if not self: return self # Combine like terms. FFPDs = sorted(self) new_FFPDs = [] temp = FFPDs[0] for f in FFPDs[1:]: if  temp.denominator_factored() == f.denominator_factored(): # Add f to temp. num = temp.numerator() + f.numerator() temp = FFPD(num, temp.denominator_factored()) else: # Append temp to new_FFPDs and update temp. new_FFPDs.append(temp) temp = f new_FFPDs.append(temp) return FFPDSum(new_FFPDs) def sum(self): r""" Return the sum of the FFPDs in self as a FFPD. EXAMPLES:: sage: from sage.combinat.amgf import * sage: R. = PolynomialRing(QQ) sage: df = (x, 1), (y, 1), (x*y + 1, 1) sage: f = FFPD(2, df) sage: g = FFPD(2*x*y, df) sage: print FFPDSum([f, g]) [(2, [(y, 1), (x, 1), (x*y + 1, 1)]), (2, [(x*y + 1, 1)])] sage: print FFPDSum([f, g]).sum() (2, [(y, 1), (x, 1)]) :: sage: R. = PolynomialRing(QQ) sage: f = FFPD(cos(x), [(x, 2)]) sage: g = FFPD(cos(y), [(x, 1), (y, 2)]) sage: print FFPDSum([f, g]) [(cos(x), [(x, 2)]), (cos(y), [(y, 2), (x, 1)])] sage: print FFPDSum([f, g]).sum() (y^2*cos(x) + x*cos(y), [(y, 2), (x, 2)]) """ if not self: return self # Compute the sum's numerator and denominator. R = self.ring() summy = sum((f.quotient() for f in self)) numer = summy.numerator() denom = R(summy.denominator()) # 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. if denom == 1: # Done return FFPD(numer, df, reduce_=False) factors = [] for f in self: factors.extend([q for (q, e) in f.denominator_factored()]) # Eliminate repeats from factors and sort. factors = sorted(list(set(factors))) # The irreducible factors of denom lie in factors. # Use this fact to build df. for q in factors: e = 0 quo, rem = denom.quo_rem(q) while rem == 0: e += 1 denom = quo quo, rem = denom.quo_rem(q) if e > 0: df.append((q, e)) return FFPD(numer, df, reduce_=False)