summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorTravis Scrimshaw <tscrimsh at umn.edu>2016-02-13 23:28:56 -0600
committerTravis Scrimshaw <tscrimsh at umn.edu>2016-02-13 23:28:56 -0600
commit64c081ad91046c6ee678cdd0ad625d96e03b5ccc (patch)
tree11a17b6ba3bb204b1ed43311c3bd7f22e5d22a3b
parentUpdated Sage version to 7.1.beta3 (diff)
parentTrac #10519: Revert "mark class as experimental" (diff)
Merge branch 'public/combinat/10519' of trac.sagemath.org:sage into public/combinat/10519
-rw-r--r--src/doc/en/reference/asymptotic/index.rst1
-rw-r--r--src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py4497
2 files changed, 4498 insertions, 0 deletions
diff --git a/src/doc/en/reference/asymptotic/index.rst b/src/doc/en/reference/asymptotic/index.rst
index 86ef7bb..1dcf205 100644
--- a/src/doc/en/reference/asymptotic/index.rst
+++ b/src/doc/en/reference/asymptotic/index.rst
@@ -64,5 +64,6 @@ Asymptotic Expansions --- Table of Contents
sage/rings/asymptotic/growth_group_cartesian
sage/rings/asymptotic/term_monoid
sage/rings/asymptotic/misc
+ sage/rings/asymptotic/asymptotics_multivariate_generating_functions
.. include:: ../footer.txt
diff --git a/src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py b/src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py
new file mode 100644
index 0000000..3e24f39
--- /dev/null
+++ b/src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py
@@ -0,0 +1,4497 @@
+r"""
+Asymptotics of Multivariate Generating Series
+
+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 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 strictly minimal multiple point of the alegbraic variety `H = 0`.
+
+The algorithms and formulas implemented here come from [RaWi2008a]_
+and [RaWi2012]_. For a general reference take a look in the book [PeWi2013].
+
+.. [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*.
+ :arxiv:`1206.4740`.
+
+.. [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 :arxiv:`0803.2914`.
+
+.. [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. :arxiv:`1009.5715`.
+
+.. [PeWi2013] Robin Pemantle and Mark C. Wilson. *Analytic
+ Combinatorics in Several Variables*, Cambridge University Press,
+ 2013.
+
+
+Introductory Examples
+=====================
+
+::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+
+A univariate smooth point example::
+
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: H = (x - 1/2)^3
+ sage: Hfac = H.factor()
+ sage: G = -1/(x + 3)/Hfac.unit()
+ sage: F = FFPD(G, Hfac)
+ sage: F
+ (-1/(x + 3), [(x - 1/2, 3)])
+ sage: alpha = [1]
+ sage: decomp = F.asymptotic_decomposition(alpha)
+ sage: decomp
+ (0, []) +
+ (-1/2*(x^2 + 6*x + 9)*r^2/(x^5 + 9*x^4 + 27*x^3 + 27*x^2)
+ - 1/2*(5*x^2 + 24*x + 27)*r/(x^5 + 9*x^4 + 27*x^3 + 27*x^2)
+ - 3*(x^2 + 3*x + 3)/(x^5 + 9*x^4 + 27*x^3 + 27*x^2),
+ [(x - 1/2, 1)])
+ sage: F1 = decomp[1]
+ sage: p = {x: 1/2}
+ sage: asy = F1.asymptotics(p, alpha, 3)
+ sage: asy
+ (8/343*(49*r^2 + 161*r + 114)*2^r, 2, 8/7*r^2 + 184/49*r + 912/343)
+ sage: F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1])
+ [((1,), 7.555555556, [7.556851312], [-0.0001714971672]),
+ ((2,), 14.74074074, [14.74052478], [0.00001465051901]),
+ ((4,), 35.96502058, [35.96501458], [1.667911934e-7]),
+ ((8,), 105.8425656, [105.8425656], [4.399565380e-11]),
+ ((16,), 355.3119534, [355.3119534], [0.0000000000])]
+
+Another smooth point example (Example 5.4 of [RaWi2008a]_)::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ 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: alpha
+ [4, 1]
+ sage: I = F.smooth_critical_ideal(alpha)
+ sage: 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([SR(z) for z in I.gens()],
+ ....: [SR(z) for z in R.gens()], solution_dict=true)
+ sage: s
+ [{y: 1, x: 1}]
+ sage: p = s[0]
+ sage: asy = F.asymptotics(p, alpha, 1, verbose=True)
+ Creating auxiliary functions...
+ Computing derivatives of auxiallary functions...
+ Computing derivatives of more auxiliary functions...
+ Computing second order differential operator actions...
+ sage: asy
+ (1/24*2^(2/3)*(sqrt(3) + 4/(sqrt(3) + I) + I)*gamma(1/3)/(pi*r^(1/3)),
+ 1,
+ 1/24*2^(2/3)*(sqrt(3) + 4/(sqrt(3) + I) + I)*gamma(1/3)/(pi*r^(1/3)))
+ sage: r = SR('r'); tuple((a*r^(1/3)).full_simplify() / r^(1/3) for a in asy) # make nicer coefficients
+ (1/12*sqrt(3)*2^(2/3)*gamma(1/3)/(pi*r^(1/3)),
+ 1,
+ 1/12*sqrt(3)*2^(2/3)*gamma(1/3)/(pi*r^(1/3)))
+ sage: F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1])
+ [((4, 1), 0.1875000000, [0.1953794675...], [-0.042023826...]),
+ ((8, 2), 0.1523437500, [0.1550727862...], [-0.017913673...]),
+ ((16, 4), 0.1221771240, [0.1230813519...], [-0.0074009592...]),
+ ((32, 8), 0.09739671811, [0.09768973377...], [-0.0030084757...]),
+ ((64, 16), 0.07744253816, [0.07753639308...], [-0.0012119297...])]
+
+A multiple point example (Example 6.5 of [RaWi2012]_)::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F
+ (1, [(x + 2*y - 1, 2), (2*x + y - 1, 2)])
+ sage: I = F.singular_ideal()
+ sage: 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: F.is_convenient_multiple_point(p)
+ (True, 'convenient in variables [x, y]')
+ sage: alpha = (var('a'), var('b'))
+ sage: decomp = F.asymptotic_decomposition(alpha); decomp
+ (0, []) +
+ (-1/9*(2*b^2*x^2 - 5*a*b*x*y + 2*a^2*y^2)*r^2/(x^2*y^2)
+ - 1/9*(6*b*x^2 - 5*(a + b)*x*y + 6*a*y^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: F1 = decomp[1]
+ sage: F1.asymptotics(p, alpha, 2)
+ (-3*((2*a^2 - 5*a*b + 2*b^2)*r^2 + (a + b)*r + 3)*((1/3)^(-a)*(1/3)^(-b))^r,
+ (1/3)^(-a)*(1/3)^(-b), -3*(2*a^2 - 5*a*b + 2*b^2)*r^2 - 3*(a + b)*r - 9)
+ sage: alpha = [4, 3]
+ sage: decomp = F.asymptotic_decomposition(alpha)
+ sage: F1 = decomp[1]
+ sage: asy = F1.asymptotics(p, alpha, 2)
+ sage: asy
+ (3*(10*r^2 - 7*r - 3)*2187^r, 2187, 30*r^2 - 21*r - 9)
+ sage: F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1])
+ [((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])]
+
+TESTS::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (1 - 2*x - y) * (1 - x - 2*y)
+ sage: G = 1
+ sage: Hfac = H.factor()
+ sage: G = G / Hfac.unit()
+ sage: F = FFPD(G, Hfac); F
+ (1, [(x + 2*y - 1, 1), (2*x + y - 1, 1)])
+ sage: p = {x: 1, y: 1}
+ sage: alpha = [1, 1]
+ sage: F.asymptotics(p, alpha, 1)
+ (1/3, 1, 1/3)
+
+::
+
+ sage: R.<x,y,t> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (1 - y) * (1 + x^2) * (1 - t*(1 + x^2 + x*y^2))
+ sage: G = (1 + x) * (1 + x^2 - x*y^2)
+ sage: Hfac = H.factor()
+ sage: G = G / Hfac.unit()
+ sage: F = FFPD(G, Hfac); F
+ (-x^2*y^2 + x^3 - x*y^2 + x^2 + x + 1,
+ [(y - 1, 1), (x^2 + 1, 1), (x*y^2*t + x^2*t + t - 1, 1)])
+ sage: p = {x: 1, y: 1, t: 1/3}
+ sage: alpha = [1, 1, 1]
+ sage: F.asymptotics_multiple(p, alpha, 1, var('r')) # not tested - see #19989
+
+
+Various
+=======
+
+AUTHORS:
+
+- Alexander Raichev (2008)
+- Daniel Krenn (2014, 2016)
+
+
+Classes and Methods
+===================
+"""
+#*****************************************************************************
+# Copyright (C) 2008 Alexander Raichev <tortoise.said@gmail.com>
+# Copyright (C) 2014, 2016 Daniel Krenn <dev@danielkrenn.at>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+# http://www.gnu.org/licenses/
+# *****************************************************************************
+
+
+# libraries
+from functools import total_ordering
+import sage
+from sage.misc.misc_c import prod
+from sage.rings.integer import Integer
+
+
+@total_ordering
+class FractionWithFactoredDenominator(sage.structure.element.RingElement):
+ r"""
+ This element represents a fraction with a factored polynomial
+ denominator. See also its parent
+ :class:`FractionWithFactoredDenominatorRing` for details.
+
+ 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` (e.g., a Sage symbolic expression).
+ An element `r` with no polynomial denominator is represented as ``(r, [])``.
+
+ INPUT:
+
+ - ``numerator`` -- an element `p`. This can be of any ring from which
+ parent's base has coercion in.
+ - ``denominator_factored`` -- 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
+ - ``reduce`` -- (optional) if ``True``, then represent
+ `p/(q_1^{e_1} \cdots q_n^{e_n})` in lowest terms, otherwise
+ this won't attempt to divide `p` by any of the `q_i`.
+
+ OUTPUT:
+
+ An element representing the rational expression
+ `p/(q_1^{e_1} \cdots q_n^{e_n})`.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: df = [x, 1], [y, 1], [x*y+1, 1]
+ sage: f = FFPD(x, df)
+ sage: f
+ (1, [(y, 1), (x*y + 1, 1)])
+ sage: ff = FFPD(x, df, reduce=False)
+ sage: ff
+ (x, [(y, 1), (x, 1), (x*y + 1, 1)])
+
+ sage: f = FFPD(x + y, [(x + y, 1)])
+ sage: f
+ (1, [])
+
+ ::
+
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
+ sage: FFPD(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.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: f = 2*y/(5*(x^3 - 1)*(y + 1))
+ sage: FFPD(f)
+ (2/5*y, [(y + 1, 1), (x - 1, 1), (x^2 + x + 1, 1)])
+
+ sage: p = 1/x^2
+ sage: q = 3*x**2*y
+ sage: qs = q.factor()
+ sage: f = FFPD(p/qs.unit(), qs)
+ sage: f
+ (1/3/x^2, [(y, 1), (x, 2)])
+
+ sage: f = FFPD(cos(x)*x*y^2, [(x, 2), (y, 1)])
+ sage: f
+ (x*y^2*cos(x), [(y, 1), (x, 2)])
+
+ sage: G = exp(x + y)
+ sage: H = (1 - 2*x - y) * (1 - x - 2*y)
+ sage: a = FFPD(G/H)
+ sage: a
+ (e^(x + y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])
+ sage: a.denominator_ring
+ Multivariate Polynomial Ring in x, y over Rational Field
+ sage: b = FFPD(G, H.factor())
+ sage: b
+ (e^(x + y), [(x + 2*y - 1, 1), (2*x + y - 1, 1)])
+ sage: b.denominator_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.<x,y> = PolynomialRing(CC)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = (x + 1)/(x*y*(x*y + 1)^2)
+ sage: FFPD(f)
+ Traceback (most recent call last):
+ ...
+ TypeError: Singular error:
+ ? not implemented
+ ? error occurred in or before STDIN line ...:
+ `def sage...=factorize(sage...);`
+
+ AUTHORS:
+
+ - Alexander Raichev (2012-07-26)
+
+ - Daniel Krenn (2014-12-01)
+ """
+ def __init__(self, parent, numerator, denominator_factored,
+ reduce=True):
+ r"""
+ Create a FractionWithFactoredDenominator instance.
+
+ See :class:`FractionWithFactoredDenominator` for details.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: df = [x, 1], [y, 1], [x*y+1, 1]
+ sage: f = FFPD(x, df)
+ sage: TestSuite(f).run()
+ """
+ super(FractionWithFactoredDenominator, self).__init__(parent)
+
+ from sage.rings.semirings.non_negative_integer_semiring import NN
+ self._numerator = parent.numerator_ring(numerator)
+ self._denominator_factored = list(
+ (parent.denominator_ring(d), NN(n))
+ for d, n in denominator_factored)
+
+ R = self.denominator_ring
+ if 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``.
+
+ OUTPUT:
+
+ The numerator.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F.numerator()
+ -e^y
+ """
+ return self._numerator
+
+
+ def denominator(self):
+ r"""
+ Return the denominator of ``self``.
+
+ OUTPUT:
+
+ The denominator (i.e., the product of the factored denominator).
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: 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 __cmp__(self, other):
+ r"""
+ Compares two elements.
+
+ INPUT:
+
+ - ``other`` -- element to compare with ``self``.
+
+ OUTPUT:
+
+ A comparison value.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x*y, [(x-1, 1), (y-2, 2)])
+ sage: g = FFPD(x, [(x-1, 1), (y-2, 2)])
+ sage: f.__cmp__(f)
+ 0
+ sage: f.__cmp__(g)
+ 1
+ """
+ return cmp(self.numerator() * other.denominator(),
+ other.numerator() * self.denominator())
+
+
+ def denominator_factored(self):
+ r"""
+ Return the factorization in ``self.denominator_ring`` of the denominator of
+ ``self`` but without the unit part.
+
+ OUTPUT:
+
+ The factored denominator as a list of tuple ``(f, m)``, where `f` is
+ a factor and `m` its multiplicity.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F.denominator_factored()
+ [(x - 1, 1), (x*y + x + y - 1, 2)]
+ """
+ return self._denominator_factored
+
+
+ @property
+ def denominator_ring(self):
+ r"""
+ Return the ring of the denominator.
+
+ OUTPUT:
+
+ A ring.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F.denominator_ring
+ Multivariate Polynomial Ring in x, y over Rational Field
+ sage: F = FFPD(G/H)
+ sage: F
+ (e^y, [(x - 1, 1), (x*y + x + y - 1, 2)])
+ sage: F.denominator_ring
+ Multivariate Polynomial Ring in x, y over Rational Field
+ """
+ return self.parent().denominator_ring
+
+
+ @property
+ def numerator_ring(self):
+ r"""
+ Return the ring of the numerator.
+
+ OUTPUT:
+
+ A ring.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F.numerator_ring
+ Symbolic Ring
+ sage: F = FFPD(G/H)
+ sage: F
+ (e^y, [(x - 1, 1), (x*y + x + y - 1, 2)])
+ sage: F.numerator_ring
+ Symbolic Ring
+ """
+ return self.parent().numerator_ring
+
+
+ def dimension(self):
+ r"""
+ Return the number of indeterminates of ``self.denominator_ring``.
+
+ OUTPUT:
+
+ An integer.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F.dimension()
+ 2
+ """
+ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
+ from sage.rings.polynomial.multi_polynomial_ring_generic import is_MPolynomialRing
+ R = self.denominator_ring
+ if is_PolynomialRing(R) or is_MPolynomialRing(R):
+ return R.ngens()
+ raise NotImplementedError('Only polynomial rings are supported as base.')
+
+
+ def quotient(self):
+ r"""
+ Convert ``self`` into a quotient.
+
+ OUTPUT:
+
+ An element.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F
+ (-e^y, [(x - 1, 1), (x*y + x + y - 1, 2)])
+ sage: 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 __repr__(self):
+ r"""
+ Return a string representation of ``self``.
+
+ OUTPUT:
+
+ A string.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x + y, [(y, 1), (x, 1)])
+ sage: f
+ (x + y, [(y, 1), (x, 1)])
+ """
+ return repr((self.numerator(), self.denominator_factored()))
+
+
+ def __eq__(self, other):
+ r"""
+ Tests for equality of the given elements (with taking care of
+ different parents by using the coercion model).
+
+ INPUT:
+
+ - ``other`` -- object to compare with ``self``.
+
+ OUTPUT:
+
+ ``True`` or ``False``.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x, [])
+ sage: f == x
+ True
+ sage: x == f
+ True
+ sage: f == 4
+ False
+ sage: g = FFPD(R(3), [])
+ sage: g == R(3)
+ True
+ sage: g == QQ(3)
+ True
+ sage: g == ZZ(3)
+ True
+ sage: 3 == g
+ True
+ """
+ from sage.structure.element import have_same_parent
+ if have_same_parent(self, other):
+ return self._eq_(other)
+
+ from sage.structure.element import get_coercion_model
+ import operator
+ try:
+ return get_coercion_model().bin_op(self, other, operator.eq)
+ except TypeError:
+ return False
+
+
+ def _eq_(self, other):
+ r"""
+ Two FFPD instances are equal iff they represent the same
+ fraction.
+
+ INPUT:
+
+ - ``other`` -- object to compare with ``self``.
+
+ OUTPUT:
+
+ ``True`` or ``False``.
+
+ It can be assumed that ``self`` and ``other`` have the same parent.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ 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.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: G = exp(x + y)
+ sage: H = (1 - 2*x - y) * (1 - x - 2*y)
+ sage: a = FFPD(G/H)
+ sage: b = FFPD(G, H.factor())
+ sage: bool(a == b)
+ True
+ """
+ return self.quotient() == other.quotient()
+
+
+ def __ne__(self, other):
+ r"""
+ Tests for nonequality of the given elements.
+
+ INPUT:
+
+ - ``other`` -- object to compare with ``self``.
+
+ OUTPUT:
+
+ ``True`` or ``False``.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ 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
+ 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``)
+ of (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).
+
+ INPUT:
+
+ - ``other`` -- object to compare with ``self``.
+
+ OUTPUT:
+
+ ``True`` or ``False``.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: df = [x, 1], [y, 1], [x*y+1, 1]
+ sage: f = FFPD(x, df); f
+ (1, [(y, 1), (x*y + 1, 1)])
+ sage: ff = FFPD(x, df, reduce=False); ff
+ (x, [(y, 1), (x, 1), (x*y + 1, 1)])
+ sage: g = FFPD(y, df)
+ sage: h = FFPD(exp(x), df)
+ sage: i = FFPD(sin(x + 2), df)
+ sage: f < ff
+ True
+ sage: f < g
+ True
+ sage: g < h
+ True
+ sage: 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
+ if this element.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ Assume that the numerator of ``self`` lies in the same
+ univariate factorial polynomial ring as the factors of the denominator.
+
+ 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:
+
+ .. MATH::
+
+ (*) \quad p_0 + \sum_{i=1}^{m} \frac{p_i}{q_i^{e_i}},
+
+ for some `p_j \in R`.
+ We call `(*)` the *usual partial fraction decomposition* of `f`.
+
+ .. NOTE::
+
+ This partial fraction decomposition can be computed using
+ - :meth:`sage.symbolic.expression.Expression.partial_fraction` or
+ - :meth:`sage.categories.quotient_fields.QuotientFields.ElementMethods.partial_fraction_decomposition`
+ as well. However, here we use the already obtained/cached
+ factorization of the denominator. This gives a speed up for
+ non-small instances.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+
+ One variable::
+
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
+ sage: 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(f).univariate_decomposition()
+ sage: decomp
+ (5*x^3, []) +
+ (1, [(x - 1, 1)]) +
+ (1, [(x, 1)]) +
+ (1/3, [(x^2 + 1/3, 1)])
+ sage: decomp.sum().quotient() == f
+ True
+
+ One variable with numerator in symbolic ring::
+
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: f = 5*x^3 + 1/x + 1/(x-1) + exp(x)/(3*x^2 + 1)
+ sage: f
+ (5*x^5 - 5*x^4 + 2*x - 1)/(x^2 - x) + e^x/(3*x^2 + 1)
+ sage: decomp = FFPD(f).univariate_decomposition()
+ sage: decomp
+ (0, []) +
+ (15/4*x^7 - 15/4*x^6 + 5/4*x^5 - 5/4*x^4 + 3/2*x^3 + 1/4*x^2*e^x -
+ 3/4*x^2 - 1/4*x*e^x + 1/2*x - 1/4, [(x - 1, 1)]) +
+ (-15*x^7 + 15*x^6 - 5*x^5 + 5*x^4 - 6*x^3 -
+ x^2*e^x + 3*x^2 + x*e^x - 2*x + 1, [(x, 1)]) +
+ (1/4*(15*x^7 - 15*x^6 + 5*x^5 - 5*x^4 + 6*x^3 + x^2*e^x -
+ 3*x^2 - x*e^x + 2*x - 1)*(3*x - 1), [(x^2 + 1/3, 1)])
+
+ One variable over a finite field::
+
+ sage: R.<x> = PolynomialRing(GF(2))
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
+ sage: f
+ (x^6 + x^4 + 1)/(x^3 + x)
+ sage: decomp = FFPD(f).univariate_decomposition()
+ sage: decomp
+ (x^3, []) + (1, [(x, 1)]) + (x, [(x + 1, 2)])
+ sage: decomp.sum().quotient() == f
+ True
+
+ One variable over an inexact field::
+
+ sage: R.<x> = PolynomialRing(CC)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1)
+ sage: 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(f).univariate_decomposition()
+ sage: 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: decomp.sum().quotient() == f # Rounding error coming
+ False
+
+ TESTS::
+
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: f = exp(x) / (x^2-x)
+ sage: f
+ e^x/(x^2 - x)
+ sage: FFPD(f).univariate_decomposition()
+ (0, []) + (e^x, [(x - 1, 1)]) + (-e^x, [(x, 1)])
+
+ AUTHORS:
+
+ - Robert Bradshaw (2007-05-31)
+ - Alexander Raichev (2012-06-25)
+ - Daniel Krenn (2014-12-01)
+ """
+ if self.dimension() > 1:
+ return FractionWithFactoredDenominatorSum([self])
+
+ R = self.denominator_ring
+ p = self.numerator()
+ q = self.denominator()
+ try:
+ whole, p = R(p).quo_rem(q)
+ mn = R.one()
+ except TypeError, ValueError:
+ whole = R(0)
+ mn = p
+ p = R.one()
+ df = self.denominator_factored()
+ decomp = [self.parent()(whole, [])]
+ denominator = prod(b**n for b, n in df)
+ for a, m in df:
+ am = a**m
+ q, r = denominator.quo_rem(am)
+ assert r==0
+ numer = p * q.inverse_mod(am) % am
+ # The inverse exists because the product and a**m
+ # are relatively prime.
+ decomp.append(self.parent()(mn * numer, [(a, m)]))
+ return FractionWithFactoredDenominatorSum(decomp)
+
+
+ def nullstellensatz_certificate(self):
+ r"""
+ Return a Nullstellensatz certificate of ``self`` if it exists.
+
+ OUTPUT:
+
+ A list of polynomials or ``None``.
+
+ 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.denominator_ring``
+ that satisfies `h_1 q_1 + \cdots + h_m q_n = 1` if it exists.
+ Otherwise return ``None``.
+ Only works for multivariate base rings.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: G = sin(x)
+ sage: H = x^2 * (x*y + 1)
+ sage: f = FFPD(G, H.factor())
+ sage: L = f.nullstellensatz_certificate()
+ sage: 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(f).nullstellensatz_certificate()
+ sage: L is None
+ True
+ """
+ R = self.denominator_ring
+ df = self.denominator_factored()
+ J = R.ideal([q ** e for q, e in df])
+ if R.one() in J:
+ return R.one().lift(J)
+ return None
+
+
+ def nullstellensatz_decomposition(self):
+ r"""
+ Return a Nullstellensatz decomposition of ``self``.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ Let `f = p/q` where `q` lies in a `d` -variate polynomial ring
+ `K[X]` for some field `K` and `d \geq 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 \mid q_i(x) = 0\}`
+ of `q_i` over the algebraic closure `L` of `K`.
+ By [Raic2012]_, `f` can be written as
+
+ .. MATH::
+
+ (*) \quad \sum_A \frac{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`.
+
+ We call `(*)` a *Nullstellensatz decomposition* of `f`.
+ Nullstellensatz decompositions are not unique.
+
+ The algorithm used comes from [Raic2012]_.
+
+ .. NOTE::
+
+ Recursive. Only works for multivariate ``self``.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import *
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 1/(x*(x*y + 1))
+ sage: decomp = FFPD(f).nullstellensatz_decomposition()
+ sage: decomp
+ (0, []) + (1, [(x, 1)]) + (-y, [(x*y + 1, 1)])
+ sage: decomp.sum().quotient() == f
+ True
+ sage: [r.nullstellensatz_certificate() is None for r in decomp]
+ [True, True, True]
+
+ ::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: G = sin(y)
+ sage: H = x*(x*y + 1)
+ sage: f = FFPD(G, H.factor())
+ sage: decomp = f.nullstellensatz_decomposition()
+ sage: decomp
+ (0, []) + (sin(y), [(x, 1)]) + (-y*sin(y), [(x*y + 1, 1)])
+ sage: bool(decomp.sum().quotient() == G/H)
+ True
+ sage: [r.nullstellensatz_certificate() is None for r in decomp]
+ [True, True, True]
+ """
+ L = self.nullstellensatz_certificate()
+ if L is None:
+ # No decomposing possible.
+ return FractionWithFactoredDenominatorSum([self])
+
+ # Otherwise decompose recursively.
+ decomp = FractionWithFactoredDenominatorSum()
+ p = self.numerator()
+ df = self.denominator_factored()
+ m = len(df)
+ iteration1 = FractionWithFactoredDenominatorSum(
+ [self.parent()(p * L[i], [df[j] for j in xrange(m) if j != i])
+ for i in xrange(m) if L[i] != 0])
+
+ # Now decompose each FFPD of iteration1.
+ 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 algebraic dependence certificate of ``self``.
+
+ OUTPUT:
+
+ An ideal.
+
+ 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.denominator_ring.base_ring()`` that has
+ ``m = len(self.denominator_factored())`` indeterminates.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 1/(x^2 * (x*y + 1) * y^3)
+ sage: ff = FFPD(f)
+ sage: J = ff.algebraic_dependence_certificate(); 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.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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(); 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(f).algebraic_dependence_certificate()
+ sage: J
+ Ideal (0) of Multivariate Polynomial Ring in T0, T1 over
+ Rational Field
+
+ ::
+
+ sage: f = sin(1)/(x^3 * y^2)
+ sage: J = FFPD(f).algebraic_dependence_certificate()
+ sage: J
+ Ideal (0) of Multivariate Polynomial Ring in T0, T1 over Rational Field
+ """
+ from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+
+ R = self.denominator_ring
+ 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``.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ 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
+
+ .. MATH::
+
+ (*) \quad \sum_A \frac{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| \leq d` and
+ `\{q_i \mid i \in A\}` is algebraically independent.
+
+ We call `(*)` an *algebraic dependence decomposition* of `f`.
+ Algebraic dependence decompositions are not unique.
+
+ The algorithm used comes from [Raic2012]_.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 1/(x^2 * (x*y + 1) * y^3)
+ sage: ff = FFPD(f)
+ sage: decomp = ff.algebraic_dependence_decomposition()
+ sage: decomp
+ (0, []) + (-x, [(x*y + 1, 1)]) +
+ (x^2*y^2 - x*y + 1, [(y, 3), (x, 2)])
+ sage: 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.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: 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
+ """
+ J = self.algebraic_dependence_certificate()
+ if not J:
+ # No decomposing possible.
+ return FractionWithFactoredDenominatorSum([self])
+
+ # Otherwise decompose recursively.
+ decomp = FractionWithFactoredDenominatorSum()
+ 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.
+ FFPD = FractionWithFactoredDenominatorRing(J.ring())
+ iteration1_temp = FractionWithFactoredDenominatorSum(
+ [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 = FractionWithFactoredDenominatorSum()
+ for r in iteration1_temp:
+ num1 = p * J.ring()(r.numerator()).subs(qpowsub)
+ denoms1 = []
+ for q, e in r.denominator_factored():
+ j = new_vars.index(q)
+ denoms1.append((df[j][0], df[j][1] * e))
+ iteration1.append(self.parent()(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``.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ 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 \mid q_i(x) = 0\}` of `q_i` over the algebraic closure
+ `L` of `K`. By [Raic2012]_, `f` can be written as
+
+ .. MATH::
+
+ (*) \quad \sum_A \frac{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 \mid 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.
+
+ We call `(*)` a *Leinartas decomposition* of `f`.
+ Leinartas decompositions are not unique.
+
+ The algorithm used comes from [Raic2012]_.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = (x^2 + 1)/((x + 2)*(x - 1)*(x^2 + x + 1))
+ sage: decomp = FFPD(f).leinartas_decomposition()
+ sage: decomp
+ (0, []) + (2/9, [(x - 1, 1)]) +
+ (-5/9, [(x + 2, 1)]) + (1/3*x, [(x^2 + x + 1, 1)])
+ sage: decomp.sum().quotient() == f
+ True
+
+ ::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 1/x + 1/y + 1/(x*y + 1)
+ sage: decomp = FFPD(f).leinartas_decomposition()
+ sage: decomp
+ (0, []) + (1, [(x*y + 1, 1)]) + (x + y, [(y, 1), (x, 1)])
+ sage: 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.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: 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.<x,y,z>= PolynomialRing(GF(2, 'a'))
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 1/(x * y * z * (x*y + z))
+ sage: decomp = FFPD(f).leinartas_decomposition()
+ sage: decomp
+ (0, []) + (1, [(z, 2), (x*y + z, 1)]) +
+ (1, [(z, 2), (y, 1), (x, 1)])
+ sage: decomp.sum().quotient() == f
+ True
+ """
+ if self.dimension() == 1:
+ # Sage's lift() function doesn't work in
+ # univariate polynomial rings.
+ # So nullstellensatz_decomposition() won't work.
+ # Can use algebraic_dependence_decomposition(),
+ # which is sufficient.
+ # temp = FractionWithFactoredDenominatorSum([self])
+ # Alternatively can use univariate_decomposition(),
+ # which is more efficient.
+ return self.univariate_decomposition()
+ temp = self.nullstellensatz_decomposition()
+ decomp = FractionWithFactoredDenominatorSum()
+ 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"""
+ Return the cohomology decomposition of ``self``.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ 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 \leq 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 \mid q_i(x) = 0 \}`, where `L` is the algebraic
+ closure of the field `K`.
+ Return a :class:`FractionWithFactoredDenominatorSum`
+ `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.
+
+ The algorithm used here comes from the proof of Theorem 17.4 of
+ [AiYu1983]_.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = 1/(x^2 + x + 1)^3
+ sage: decomp = FFPD(f).cohomology_decomposition()
+ sage: decomp
+ (0, []) + (2/3, [(x^2 + x + 1, 1)])
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: FFPD(1, [(x, 1), (y, 2)]).cohomology_decomposition()
+ (0, [])
+
+ sage: p = 1
+ sage: qs = [(x*y - 1, 1), (x**2 + y**2 - 1, 2)]
+ sage: f = FFPD(p, qs)
+ sage: 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)])
+ """
+ from sage.calculus.functions import jacobian
+ from sage.arith.all import xgcd
+ from sage.sets.set import Set
+ from sage.symbolic.ring import SR
+
+ R = self.denominator_ring
+ df = self.denominator_factored()
+ n = len(df)
+ if sum([e for (q, e) in df]) <= n:
+ # No decomposing possible.
+ return FractionWithFactoredDenominatorSum([self])
+
+ # Otherwise decompose recursively.
+ decomp = FractionWithFactoredDenominatorSum()
+ 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.
+ if self.dimension() == 1:
+ # Sage's lift() function doesn't work in
+ # univariate polynomial rings.
+ # So use xgcd(), which does the same thing in this case.
+ # Note that by assumption qs and dets have length 1.
+ L = xgcd(qs[0], dets[0])[1:]
+ else:
+ L = R(1).lift(R.ideal(qs + dets))
+
+ # Do first iteration of decomposition.
+ iteration1 = FractionWithFactoredDenominatorSum()
+ # 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(self.parent()(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 = FractionWithFactoredDenominatorRing._permutation_sign(x, X)
+ iteration1.append(self.parent()((-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()
+
+
+ def asymptotic_decomposition(self, alpha, asy_var=None):
+ r"""
+ Return a sum
+ that has the same asymptotic expansion
+ as ``self`` in the direction ``alpha`` but each summand has a
+ denominator factorization of the form `[(q_1, 1), \ldots, (q_n, 1)]`,
+ where ``n`` is at most ``d = self.dimension()``.
+
+ INPUT:
+
+ - ``alpha`` -- a `d`-tuple of positive integers or symbolic variables
+ - ``asy_var`` -- (default: ``None``) a symbolic variable with
+ respect to which to compute asymptotics;
+ if ``None`` is given, we set ``asy_var = var('r')``
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ The output results from a Leinartas decomposition followed by a
+ cohomology decomposition.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: f = (x^2 + 1)/((x - 1)^3*(x + 2))
+ sage: F = FFPD(f)
+ sage: alpha = [var('a')]
+ sage: F.asymptotic_decomposition(alpha)
+ (0, []) +
+ (1/54*(5*a^2*x^2 + 2*a^2*x + 11*a^2)*r^2/x^2
+ - 1/54*(5*a*x^2 - 2*a*x - 33*a)*r/x^2 + 11/27/x^2, [(x - 1, 1)]) +
+ (-5/27, [(x + 2, 1)])
+
+ ::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F.asymptotic_decomposition(alpha)
+ (0, []) +
+ (1/3*(2*b*x - a*y)*r/(x*y) + 1/3*(2*x - y)/(x*y),
+ [(x + 2*y - 1, 1), (2*x + y - 1, 1)])
+ """
+ from sage.calculus.var import var
+ from sage.symbolic.ring import SR
+
+ R = self.denominator_ring
+ 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 = FractionWithFactoredDenominatorSum([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 = FractionWithFactoredDenominatorSum()
+ for f in decomp1:
+ ff = self.parent()(f.numerator() * cauchy_stuff,
+ f.denominator_factored())
+ decomp2.extend(ff.cohomology_decomposition())
+ decomp2 = decomp2._combine_like_terms_()
+
+ # Divide out cauchy_stuff from integrands.
+ decomp3 = FractionWithFactoredDenominatorSum()
+ for f in decomp2:
+ ff = self.parent()((f.numerator() /
+ cauchy_stuff).simplify_full().collect(asy_var),
+ f.denominator_factored())
+ decomp3.append(ff)
+
+ return decomp3
+
+
+ def asymptotics(self, p, alpha, N, asy_var=None, numerical=0,
+ verbose=False):
+ r"""
+ Return the asymptotics in the given direction.
+
+ INPUT:
+
+ - ``p`` -- a dictionary with keys that can be coerced to equal
+ ``self.denominator_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
+ - ``asy_var`` -- (default: ``None``) a symbolic variable for the
+ asymptotic expansion; if ``none`` is given, then
+ ``var('r')`` will be assigned
+ - ``numerical`` -- (default: 0) a natural number;
+ if ``numerical`` is greater than 0, then return a numerical
+ approximation of `F_{r \alpha}` with ``numerical`` digits of
+ precision; otherwise return exact values
+ - ``verbose`` -- (default: ``False``) print the current state of
+ the algorithm
+
+ 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``.
+
+ This function
+ returns 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` is ``asy_var`` and ``alpha`` is a tuple of
+ positive integers of length `d` which is ``self.dimension()``.
+ Assume that
+
+ - `F` is holomorphic in a neighborhood of the origin;
+ - the unique factorization of the denominator `H` of `F` in the local
+ algebraic ring at `p` equals its unique factorization in the local
+ analytic ring at `p`;
+ - the unique factorization of `H` in the local algebraic ring at `p`
+ has at most ``d`` irreducible factors, none of which are repeated
+ (one can reduce to this case via :meth:`asymptotic_decomposition()`);
+ - `p` is a convenient strictly minimal smooth or multiple point
+ with all nonzero coordinates that is critical and nondegenerate
+ for ``alpha``.
+
+ The algorithms used here come from [RaWi2008a]_ and [RaWi2012]_.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+
+ A smooth point example::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: decomp = F.asymptotic_decomposition(alpha); decomp
+ (0, []) + (-3/2*r*(y + 1)/y - 1/2*(y + 1)/y, [(x*y + x + y - 1, 1)])
+ sage: F1 = decomp[1]
+ sage: p = {y: 1/3, x: 1/2}
+ sage: asy = F1.asymptotics(p, alpha, 2, verbose=True)
+ Creating auxiliary functions...
+ Computing derivatives of auxiallary functions...
+ Computing derivatives of more auxiliary functions...
+ Computing second order differential operator actions...
+ sage: asy
+ (1/6000*(3600*sqrt(5)*sqrt(3)*sqrt(2)*sqrt(r)/sqrt(pi)
+ + 463*sqrt(5)*sqrt(3)*sqrt(2)/(sqrt(pi)*sqrt(r)))*432^r,
+ 432,
+ 3/5*sqrt(5)*sqrt(3)*sqrt(2)*sqrt(r)/sqrt(pi)
+ + 463/6000*sqrt(5)*sqrt(3)*sqrt(2)/(sqrt(pi)*sqrt(r)))
+ sage: F.relative_error(asy[0], alpha, [1, 2, 4, 8, 16], asy[1])
+ [((4, 3), 2.083333333, [2.092576110], [-0.0044365330...]),
+ ((8, 6), 2.787374614, [2.790732875], [-0.0012048112...]),
+ ((16, 12), 3.826259447, [3.827462310], [-0.0003143703...]),
+ ((32, 24), 5.328112821, [5.328540787], [-0.0000803222...]),
+ ((64, 48), 7.475927885, [7.476079664], [-0.0000203023...])]
+
+ A multiple point example::
+
+ sage: R.<x,y,z>= PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: F
+ (-16, [(x + 2*y + z - 4, 1), (2*x + y + z - 4, 2)])
+ sage: alpha = [3, 3, 2]
+ sage: decomp = F.asymptotic_decomposition(alpha); decomp
+ (0, []) +
+ (16*r*(4*y - 3*z)/(y*z) + 16*(2*y - z)/(y*z),
+ [(x + 2*y + z - 4, 1), (2*x + y + z - 4, 1)])
+ sage: F1 = decomp[1]
+ sage: p = {x: 1, y: 1, z: 1}
+ sage: asy = F1.asymptotics(p, alpha, 2, verbose=True) # long time
+ Creating auxiliary functions...
+ Computing derivatives of auxiliary functions...
+ Computing derivatives of more auxiliary functions...
+ Computing second-order differential operator actions...
+ sage: asy # long time
+ (4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(r)),
+ 1, 4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(r)))
+ sage: F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1]) # long time
+ [((3, 3, 2), 0.9812164307, [1.515572606], [-0.54458543...]),
+ ((6, 6, 4), 1.576181132, [1.992989399], [-0.26444185...]),
+ ((12, 12, 8), 2.485286378, [2.712196351], [-0.091301338...]),
+ ((24, 24, 16), 3.700576827, [3.760447895], [-0.016178847...])]
+ """
+ from sage.calculus.functional import diff
+ from sage.calculus.var import var
+
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._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()
+ n = len(df) # Number of smooth factors
+
+ # 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
+
+ if n == 1:
+ # Smooth point.
+ return self.asymptotics_smooth(p, alpha, N, asy_var, coordinate,
+ numerical, verbose=verbose)
+
+ # Multiple point.
+ return self.asymptotics_multiple(p, alpha, N, asy_var, coordinate,
+ numerical, verbose=verbose)
+
+
+ def asymptotics_smooth(self, p, alpha, N, asy_var, coordinate=None,
+ numerical=0, verbose=False):
+ r"""
+ Same as :meth:`asymptotics()`, but only in the case of a
+ convenient smooth point.
+
+ INPUT:
+
+ - ``p`` -- a dictionary with keys that can be coerced to equal
+ ``self.denominator_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, \ldots, 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 is greater than 0, then return a numerical approximation
+ of the Maclaurin ray coefficients of ``self`` with ``numerical``
+ digits of precision; otherwise return exact values
+ - ``verbose`` -- (default: ``False``) print the current state of
+ the algorithm
+
+ OUTPUT:
+
+ The asymptotic expansion.
+
+ 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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = 2 - 3*x
+ sage: Hfac = H.factor()
+ sage: G = 1/Hfac.unit()
+ sage: F = FFPD(G, Hfac)
+ sage: F
+ (-1/3, [(x - 2/3, 1)])
+ sage: alpha = [2]
+ sage: p = {x: 2/3}
+ sage: asy = F.asymptotics_smooth(p, alpha, 3, asy_var=var('r'))
+ sage: asy
+ (1/2*(9/4)^r, 9/4, 1/2)
+
+ ::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ 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: F.asymptotics_smooth(p, alpha, 2, var('r'), numerical=3, verbose=True)
+ Creating auxiliary functions...
+ Computing derivatives of auxiallary functions...
+ Computing derivatives of more auxiliary functions...
+ Computing second order differential operator actions...
+ (71.2^r*(0.369/sqrt(r) - 0.018.../r^(3/2)), 71.2, 0.369/sqrt(r) - 0.018.../r^(3/2))
+
+ 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: alpha
+ [4, 1]
+ sage: p = {x: 1, y: 1}
+ sage: F.asymptotics_smooth(p, alpha, 5, var('r'), verbose=True) # not tested (140 seconds)
+ Creating auxiliary functions...
+ Computing derivatives of auxiallary functions...
+ Computing derivatives of more auxiliary functions...
+ Computing second order differential operator actions...
+ (1/12*sqrt(3)*2^(2/3)*gamma(1/3)/(pi*r^(1/3))
+ - 1/96*sqrt(3)*2^(1/3)*gamma(2/3)/(pi*r^(5/3)),
+ 1,
+ 1/12*sqrt(3)*2^(2/3)*gamma(1/3)/(pi*r^(1/3))
+ - 1/96*sqrt(3)*2^(1/3)*gamma(2/3)/(pi*r^(5/3)))
+ """
+ from sage.calculus.functional import diff
+ from sage.calculus.functions import jacobian
+ from sage.calculus.var import function
+ from sage.calculus.var import var
+ from sage.combinat.tuple import UnorderedTuples
+ from sage.functions.other import factorial
+ from sage.functions.other import gamma
+ from sage.functions.other import sqrt
+ from sage.functions.log import exp
+ from sage.functions.log import log
+ from sage.matrix.constructor import matrix
+ from sage.modules.free_module_element import vector
+ from sage.symbolic.constants import pi
+ from sage.symbolic.relation import solve
+ from sage.symbolic.ring import SR
+ from sage.rings.all import CC
+ from sage.rings.rational_field import QQ
+
+ R = self.denominator_ring
+ 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)
+
+ # Deal with the simple univariate case first.
+ # Same as the multiple point case with n == d.
+ # but with a negative sign.
+ # I'll just past the code from the multiple point case.
+ if d == 1:
+ det = jacobian(H, X).subs(p).determinant().abs()
+ exp_scale = prod([(p[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)])
+ subexp_part = -G.subs(p) / (det * prod(p.values()))
+ 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)
+
+ # 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.
+ if verbose:
+ 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)
+ At = A.subs(e)
+ Phit = (-log(P[X[d - 1]] * ht) +
+ I * sum([alpha[i] / alpha[d - 1] * T[i]
+ for i in xrange(d - 1)]))
+ Tstar = {t: Integer(0) for t in T}
+ # Store h and U and all their derivatives evaluated at P.
+ atP = P.copy()
+ atP.update({h.subs(P): Integer(1) / P[X[d - 1]]})
+
+ # 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 = FractionWithFactoredDenominatorRing._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 = FractionWithFactoredDenominatorRing._diff_all(
+ H, X, 2 * N, ending=[X[d - 1]], sub_final=P)
+ if verbose:
+ print("Computing derivatives of auxiallary functions...")
+ # For convenience in checking if all the nontrivial derivatives of U
+ # at p are zero a few line below, store the value of U(p) in atP
+ # instead of in Uderivs.
+ Uderivs = {}
+ 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 = FractionWithFactoredDenominatorRing._diff_prod(
+ Hderivs, U, Hcheck, X,
+ range(1, k + 1), end, Uderivs, atP)
+ # Check for a nonzero U derivative.
+ if any(u for u in Uderivs.values()):
+ all_zero = False
+ if all_zero:
+ # Then, using a proposition at the end of [RaWi2012], we can
+ # conclude that all higher derivatives of U are zero.
+ for l in xrange(k + 1, 2 * N + 1):
+ for s in UnorderedTuples(X, l):
+ Uderivs.update({diff(U, s).subs(P): Integer(0)})
+ else:
+ # Have to compute the rest of the derivatives.
+ Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
+ Hderivs, U, Hcheck, X,
+ range(k + 1, 2 * N + 1), end, Uderivs,
+ atP)
+ else:
+ Uderivs = FractionWithFactoredDenominatorRing._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.
+ if verbose:
+ print("Computing derivatives of more auxiliary functions...")
+ AA = function('AA')(t)
+ BB = function('BB')(t)
+ if v.mod(2) == 0:
+ At_derivs = FractionWithFactoredDenominatorRing._diff_all(
+ At, T, 2 * N - 2,
+ sub=hderivs1, sub_final=[Tstar, atP],
+ rekey=AA)
+ Phitu_derivs = FractionWithFactoredDenominatorRing._diff_all(
+ Phitu, T, 2 * N - 2 +v,
+ sub=hderivs1,
+ sub_final=[Tstar, atP],
+ zero_order=v + 1, rekey=BB)
+ else:
+ At_derivs = FractionWithFactoredDenominatorRing._diff_all(
+ At, T, N - 1, sub=hderivs1,
+ sub_final=[Tstar, atP], rekey=AA)
+ Phitu_derivs = FractionWithFactoredDenominatorRing._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)
+ if verbose:
+ print("Computing second order differential operator actions...")
+ DD = FractionWithFactoredDenominatorRing._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(sum([(-1) ** l * gamma((2 * k + v * l + 1) / v) /
+ (factorial(l) * factorial(2 * k + v * l)) *
+ DD[(k, l)] for l in xrange(0, 2 * k + 1)]))
+ chunk = (a ** (-1 / v) / (pi * v) *
+ sum([alpha[d - 1] ** (-(2 * k + 1) / v) *
+ L[k] * asy_var ** (-(2 * k + 1) / v)
+ for k in xrange(N)]))
+ else:
+ zeta = exp(I * pi / (2 * v))
+ for k in xrange(N):
+ L.append(sum([(-1) ** l * gamma((k + v * l + 1) / v) /
+ (factorial(l) * factorial(k + v * l)) *
+ (zeta ** (k + v * l + 1) +
+ (-1) ** (k + v * l) *
+ zeta ** (-(k + v * l + 1))) *
+ DD[(k, l)] for l in xrange(0, k + 1)]))
+ chunk = (abs(a) ** (-1 / v) / (2 * pi * v) *
+ sum([alpha[d - 1] ** (-(k + 1) / v) *
+ L[k] * asy_var ** (-(k + 1) / v)
+ for k in xrange(N)]))
+
+ # Asymptotics for d >= 2 case.
+ # A singular Phit''(Tstar) will cause a crash in this case.
+ else:
+ Phit1 = jacobian(Phit, T).subs(hderivs1)
+ a = jacobian(Phit1, T).subs(hderivs1).subs(Tstar).subs(atP)
+ a_inv = a.inverse()
+ Phitu = (Phit - (1 / QQ(2)) * matrix([T]) *
+ a * matrix([T]).transpose())
+ Phitu = Phitu[0][0]
+ # Compute all partial derivatives of At and Phitu up to
+ # orders 2 * N-2 and 2 * N, respectively.
+ # 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.
+ if verbose:
+ print("Computing derivatives of more auxiliary functions...")
+ AA = function('AA')(*tuple(T))
+ At_derivs = FractionWithFactoredDenominatorRing._diff_all(
+ At, T, 2 * N - 2, sub=hderivs1,
+ sub_final =[Tstar, atP], rekey=AA)
+ BB = function('BB')(*tuple(T))
+ Phitu_derivs = FractionWithFactoredDenominatorRing._diff_all(
+ Phitu, T, 2 * N, sub=hderivs1,
+ 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)
+ if verbose:
+ print("Computing second order differential operator actions...")
+ DD = FractionWithFactoredDenominatorRing._diff_op(
+ AA, BB, AABB_derivs, T, a_inv, 1 , N)
+
+ # Plug above into asymptotic formula.
+ L = []
+ for k in xrange(N):
+ L.append(sum([DD[(0, k, l)] / ((-1) ** k * 2 ** (l + k) *
+ factorial(l) * factorial(l + k))
+ for l in xrange(0, 2 * k + 1)]))
+ chunk = sum([(2 * pi) ** ((1 - d) / Integer(2)) *
+ a.determinant() ** (-Integer(1) / Integer(2)) *
+ alpha[d - 1] ** ((Integer(1) - d) / Integer(2) - k) *
+ L[k] *
+ asy_var ** ((Integer(1) - d) / Integer(2) - k)
+ for k in xrange(N)])
+
+ chunk = chunk.subs(p).simplify()
+ coeffs = chunk.coefficients(asy_var)
+ coeffs.reverse()
+ coeffs = coeffs[:N]
+ if numerical:
+ subexp_part = sum([co[0].subs(p).n(digits=numerical) *
+ asy_var ** co[1] for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)]).n(digits=numerical)
+ else:
+ subexp_part = sum([co[0].subs(p) * asy_var ** co[1]
+ for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)])
+ return (exp_scale ** asy_var * subexp_part, exp_scale, subexp_part)
+
+
+ def asymptotics_multiple(self, p, alpha, N, asy_var, coordinate=None,
+ numerical=0, verbose=False):
+ r"""
+ Same as :meth:`asymptotics`, but only in the case of a
+ convenient multiple point nondegenerate for ``alpha``.
+ Assume also that ``self.dimension >= 2`` and that the
+ ``p.values()`` are not symbolic variables.
+
+ INPUT:
+
+ - ``p`` -- a dictionary with keys that can be coerced to equal
+ ``self.denominator_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, \ldots, 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 is greater than 0, then return a numerical approximation
+ of the Maclaurin ray coefficients of ``self`` with ``numerical``
+ digits of precision; otherwise return exact values
+ - ``verbose`` -- (default: ``False``) print the current state of
+ the algorithm
+
+ OUTPUT:
+
+ The asymptotic expansion.
+
+ The formulas used for computing the asymptotic expansion are
+ Theorem 3.4 and Theorem 3.7 of [RaWi2012]_.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y,z>= PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (4 - 2*x - y - z)*(4 - x -2*y - z)
+ sage: Hfac = H.factor()
+ sage: G = 16/Hfac.unit()
+ sage: F = FFPD(G, Hfac)
+ sage: 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: F.asymptotics_multiple(p, alpha, 2, var('r'), verbose=True) # 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: 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: 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: F.asymptotics_multiple(p, alpha, 2, var('r'), coordinate=1, verbose=True) # long time
+ Creating auxiliary functions...
+ Computing derivatives of auxiliary functions...
+ Computing derivatives of more auxiliary functions...
+ Computing second-order differential operator actions...
+ (1/172872*108^r*(24696*sqrt(7)*sqrt(3)/(sqrt(pi)*sqrt(r))
+ - 1231*sqrt(7)*sqrt(3)/(sqrt(pi)*r^(3/2))),
+ 108,
+ 1/7*sqrt(7)*sqrt(3)/(sqrt(pi)*sqrt(r))
+ - 1231/172872*sqrt(7)*sqrt(3)/(sqrt(pi)*r^(3/2)))
+
+ ::
+
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: 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: F.asymptotics_multiple(p, alpha, 2, var('r')) # long time
+ (3*((1/3)^(-a)*(1/3)^(-b))^r*e^(2/3), (1/3)^(-a)*(1/3)^(-b), 3*e^(2/3))
+ """
+ from itertools import product
+ from sage.calculus.functional import diff
+ from sage.calculus.functions import jacobian
+ from sage.calculus.var import function
+ from sage.calculus.var import var
+ from sage.combinat.combinat import stirling_number1
+ from sage.combinat.tuple import UnorderedTuples
+ from sage.functions.log import exp
+ from sage.functions.log import log
+ from sage.functions.other import factorial
+ from sage.functions.other import sqrt
+ from sage.matrix.constructor import matrix
+ from sage.misc.mrange import xmrange
+ from sage.modules.free_module_element import vector
+ from sage.rings.all import CC
+ from sage.rings.arith import binomial
+ from sage.rings.rational_field import QQ
+ from sage.symbolic.constants import pi
+ from sage.symbolic.relation import solve
+ from sage.symbolic.ring import SR
+
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._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 = {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.
+ if verbose:
+ 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]
+ hsumt = (sum([S[j] * ht[j] for j in xrange(n - 1)]) +
+ (Integer(1) - sum(S)) * ht[n - 1])
+ At = [AA.subs(e).subs({X[d - 1]: hsumt}) for AA in A]
+ Phit = (-log(P[X[d - 1]] * hsumt) +
+ I * sum([alpha[i] / alpha[d - 1] * T[i]
+ for i in xrange(d - 1)]))
+ # atP Stores h and U and all their derivatives evaluated at C.
+ atP = P.copy()
+ atP.update(dict([(hh.subs(P), Integer(1) / P[X[d - 1]]) for hh in h]))
+
+ # Compute the derivatives of h up to order 2 * N and evaluate at P.
+ hderivs1 = {} # First derivatives of h.
+ for (i, j) in xmrange([d - 1, n], tuple):
+ s = solve(diff(H[j].subs({X[d - 1]: Integer(1) / h[j]}), X[i]),
+ diff(h[j], X[i]))[0].rhs().simplify()
+ hderivs1.update({diff(h[j], X[i]): s})
+ atP.update({diff(h[j], X[i]).subs(P): s.subs(P).subs(atP)})
+ hderivs = FractionWithFactoredDenominatorRing._diff_all(
+ h, X[0:d - 1], 2 * N, sub=hderivs1, rekey=h)
+ 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.
+ if verbose:
+ print("Computing derivatives of auxiliary functions...")
+ m = min(n, N)
+ end = [X[d-1] for j in xrange(n)]
+ Hprodderivs = FractionWithFactoredDenominatorRing._diff_all(
+ Hprod, X, 2 * N - 2 + n, ending=end,
+ 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 = FractionWithFactoredDenominatorRing._diff_prod(
+ Hprodderivs, U, Hcheck, X,
+ range(1, k + 1), end, Uderivs, atP)
+ # Check for a nonzero U derivative.
+ if any(u for u in Uderivs.values()):
+ all_zero = False
+ if all_zero:
+ # Then all higher derivatives of U are zero.
+ for l in xrange(k + 1, 2 * N - 2 + m):
+ for s in UnorderedTuples(X, l):
+ Uderivs.update({diff(U, s).subs(P): Integer(0)})
+ else:
+ # Have to compute the rest of the derivatives.
+ Uderivs = FractionWithFactoredDenominatorRing._diff_prod(
+ Hprodderivs, U, Hcheck, X,
+ range(k + 1, 2 * N - 2 + m), end,
+ Uderivs, atP)
+ else:
+ Uderivs = FractionWithFactoredDenominatorRing._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 - (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.
+ if verbose:
+ print("Computing derivatives of more auxiliary functions...")
+ AA = [function('A' + str(j))(*tuple(T + S)) for j in xrange(n)]
+ At_derivs = FractionWithFactoredDenominatorRing._diff_all(
+ At, T + S, 2 * N - 2, sub=hderivs1,
+ sub_final=[thetastar, atP], rekey=AA)
+ BB = function('BB')(*tuple(T + S))
+ Phitu_derivs = FractionWithFactoredDenominatorRing._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)
+
+ if verbose:
+ print("Computing second-order differential operator actions...")
+ DD = FractionWithFactoredDenominatorRing._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)] = sum([DD[(j, k, l)] / ((-1) ** k * 2 ** (k + l) *
+ factorial(l) *
+ factorial(k + l))
+ for l in xrange(2 * k + 1)])
+ det = (a.determinant() ** (-1 / Integer(2)) *
+ (2 * pi) ** ((n - d) / Integer(2)))
+ chunk = det * sum([(alpha[d - 1] * asy_var) ** ((n - d) /
+ Integer(2) - q) *
+ sum([L[(j, k)] * binomial(n - 1, j) *
+ stirling_number1(n - j, n + k - q) *
+ (-1) ** (q - j - k)
+ for (j, k) in product(xrange(min(n - 1, q) + 1),
+ xrange(max(0, q - n),
+ q + 1))
+ if j + k <= q])
+ for q in xrange(N)])
+ chunk = chunk.subs(P).simplify()
+ coeffs = chunk.coefficients(asy_var)
+ coeffs.reverse()
+ coeffs = coeffs[:N]
+ if numerical:
+ subexp_part = sum([co[0].subs(p).n(digits=numerical) *
+ asy_var ** co[1]
+ for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)]).n(digits=numerical)
+ else:
+ subexp_part = sum([co[0].subs(p) * asy_var ** co[1]
+ for co in coeffs])
+ exp_scale = prod([(P[X[i]] ** (-alpha[i])).subs(p)
+ for i in xrange(d)])
+ return (exp_scale ** asy_var * subexp_part, exp_scale, subexp_part)
+
+
+ def _crit_cone_combo(self, p, alpha, coordinate=None):
+ r"""
+ Return an auxiliary point associated to the multiple
+ point ``p`` of the factors ``self``.
+
+ INPUT:
+
+ - ``p`` -- a dictionary with keys that can be coerced to equal
+ ``self.denominator_ring.gens()``
+ - ``alpha`` -- a list of rationals
+
+ OUTPUT:
+
+ A solution of the matrix equation `y \Gamma = \alpha^{\prime}` for `y`,
+ where `\Gamma` is the matrix given by
+ ``[FractionWithFactoredDenominatorRing.direction(v) for v in self.log_grads(p)]`` and
+ `\alpha^{\prime}` is ``FractionWithFactoredDenominatorRing.direction(alpha)``
+
+ .. NOTE::
+
+ For internal use by :meth:`asymptotics_multiple()`.
+
+ .. NOTE::
+
+ Use this function only when `\Gamma` is well-defined and
+ there is a unique solution to the matrix equation
+ `y \Gamma = \alpha'`. Fails otherwise.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ 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: f._crit_cone_combo(p, alpha)
+ [1/3*(2*a - b)/b, -2/3*(a - 2*b)/b]
+ """
+ from sage.calculus.var import var
+ from sage.matrix.constructor import matrix
+ from sage.symbolic.relation import solve
+
+ # Assuming here that each log_grads(f) has nonzero final component.
+ # Then 'direction' will not throw a division by zero error.
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._coerce_point_(R, p)
+
+ d = self.dimension()
+ n = len(self.denominator_factored())
+ Gamma = matrix(
+ [FractionWithFactoredDenominatorRing.direction(v, coordinate)
+ for v in self.log_grads(p)])
+ beta = FractionWithFactoredDenominatorRing.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]
+
+
+ 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.denominator_ring``
+
+ OUTPUT:
+
+ A list.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: p = exp(x)
+ sage: df = [(x^3 + 3*y^2, 5), (x*y, 2), (y, 1)]
+ sage: f = FFPD(p, df)
+ sage: f
+ (e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)])
+ sage: R.gens()
+ (x, y)
+ sage: p = None
+ sage: f.grads(p)
+ [(0, 1), (y, x), (3*x^2, 6*y)]
+
+ sage: p = {x: sqrt(2), y: var('a')}
+ sage: f.grads(p)
+ [(0, 1), (a, sqrt(2)), (6, 6*a)]
+ """
+ from sage.calculus.functional import diff
+
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._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.denominator_ring``
+
+ OUTPUT:
+
+ A list.
+
+ 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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: p = exp(x)
+ sage: df = [(x^3 + 3*y^2, 5), (x*y, 2), (y, 1)]
+ sage: f = FFPD(p, df)
+ sage: f
+ (e^x, [(y, 1), (x*y, 2), (x^3 + 3*y^2, 5)])
+ sage: R.gens()
+ (x, y)
+ sage: p = None
+ sage: 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: f.log_grads(p)
+ [(0, a), (sqrt(2)*a, sqrt(2)*a), (6*sqrt(2), 6*a^2)]
+ """
+ from sage.calculus.functional import diff
+
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._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.denominator_ring.gens()`` and values in a field
+ - ``coordinate`` -- (optional; default: ``None``) a natural number
+
+ OUTPUT:
+
+ A list of vectors.
+
+ This list of vectors 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`` we have
+ ``self.log_grads()[i][j] != 0`` and set ``coordinate = j``.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y,z> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: G = 1
+ sage: H = (1 - x*(1 + y)) * (1 - z*x**2*(1 + 2*y))
+ sage: Hfac = H.factor()
+ sage: G = 1/Hfac.unit()
+ sage: F = FFPD(G, Hfac)
+ sage: p = {x: 1/2, y: 1, z: 4/3}
+ sage: F.critical_cone(p)
+ ([(2, 1, 0), (3, 1, 3/2)], 2-d cone in 3-d lattice N)
+ """
+ from sage.geometry.cone import Cone
+
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._coerce_point_(R, p)
+
+ 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 = [FractionWithFactoredDenominatorRing.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"""
+ Tests if ``p`` is a convenient multiple point of ``self``.
+
+ INPUT:
+
+ - ``p`` -- a dictionary with keys that can be coerced to equal
+ ``self.denominator_ring.gens()``
+
+ OUTPUT:
+
+ A pair ``(verdict, comment)``.
+
+ In case ``p`` is a convenient multiple point, ``verdict = True`` and
+ ``comment`` is a string stating which variables it's convenient to use.
+ In case ``p`` is not, ``verdict = False`` and ``comment`` is a string
+ explaining why ``p`` fails to be a convenient multiple point.
+
+ See [RaWi2012]_ for more details.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y,z> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (1 - x*(1 + y)) * (1 - z*x**2*(1 + 2*y))
+ sage: 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: F.is_convenient_multiple_point(p1)
+ (True, 'convenient in variables [x, y]')
+ sage: F.is_convenient_multiple_point(p2)
+ (False, 'not a singular point')
+ """
+ from sage.combinat.subset import Subsets
+ from sage.matrix.constructor import matrix
+
+ R = self.denominator_ring
+
+ # Coerce keys of p into R.
+ p = FractionWithFactoredDenominatorRing._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()
+ convenientX = [X[i] for i in convenient_coordinates]
+ return (True, 'convenient in variables {}'.format(convenientX))
+
+
+ def singular_ideal(self):
+ r"""
+ Return the singular ideal of ``self``.
+
+ OUTPUT:
+
+ An ideal.
+
+ 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.rings.asymptotic.asymptotics_multivariate_generating_functions import *
+ sage: R.<x,y,z> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (1 - x*(1 + y))^3 * (1 - z*x**2*(1 + 2*y))
+ sage: 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
+ """
+ R = self.denominator_ring
+
+ 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"""
+ Return the smooth critical ideal of ``self``.
+
+ OUTPUT:
+
+ An ideal.
+
+ 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`.
+
+ See [RaWi2012]_ for more details.
+
+ INPUT:
+
+ - ``alpha`` -- a tuple of positive integers and/or symbolic entries
+ of length ``self.denominator_ring.ngens()``
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (1 - x - y - x*y)^2
+ sage: 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: 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
+ """
+ from sage.calculus.functional import diff
+ from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+ from sage.sets.set import Set
+ from sage.symbolic.ring import SR
+
+ R = self.denominator_ring
+ 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 = sorted(Set(indets)) # Delete duplicates in indets.
+ if indets:
+ L = PolynomialRing(K, indets).fraction_field()
+ 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"""
+ Return the Maclaurin coefficients of ``self`` with given
+ ``multi_indices``.
+
+ INPUT:
+
+ - ``multi_indices`` -- a list of tuples of positive integers, where
+ each tuple has length ``self.dimension()``
+ - ``numerical`` -- (optional; default: 0) a natural number; if
+ positive, return numerical approximations of coefficients with
+ ``numerical`` digits of accuracy
+
+ OUTPUT:
+
+ A dictionary whose value of the key ``nu`` are the Maclaurin
+ coefficient of index ``n``u of ``self``.
+
+ .. NOTE::
+
+ Uses iterated univariate Maclaurin expansions. Slow.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = 2 - 3*x
+ sage: Hfac = H.factor()
+ sage: G = 1 / Hfac.unit()
+ sage: F = FFPD(G, Hfac)
+ sage: F
+ (-1/3, [(x - 2/3, 1)])
+ sage: F.maclaurin_coefficients([(2*k,) for k in range(6)])
+ {(0,): 1/2,
+ (2,): 9/8,
+ (4,): 81/32,
+ (6,): 729/128,
+ (8,): 6561/512,
+ (10,): 59049/2048}
+
+ ::
+
+ sage: R.<x,y,z> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: H = (4 - 2*x - y - z) * (4 - x - 2*y - z)
+ sage: Hfac = H.factor()
+ sage: G = 16 / Hfac.unit()
+ sage: F = FFPD(G, Hfac)
+ sage: alpha = vector([3, 3, 2])
+ sage: interval = [1, 2, 4]
+ sage: S = [r*alpha for r in interval]
+ sage: F.maclaurin_coefficients(S, numerical=10)
+ {(3, 3, 2): 0.7849731445,
+ (6, 6, 4): 0.7005249476,
+ (12, 12, 8): 0.5847732654}
+ """
+ from sage.symbolic.ring import SR
+
+ R = self.denominator_ring
+ d = self.dimension()
+ coeffs = {}
+
+ # Deal with the simple univariate case first.
+ if d == 1:
+ f = SR(self.quotient())
+ x = SR(R.gens()[0])
+ m = max(multi_indices)[0]
+ f = f.taylor(x, 0, m)
+ F = R(f)
+ tmp = F.coefficients()
+ for nu in multi_indices:
+ val = tmp[nu[0]]
+ if numerical:
+ val = val.n(digits=numerical)
+ coeffs[tuple(nu)] = val
+ return coeffs
+
+ # 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.
+ X = R.gens()
+ for nu in multi_indices:
+ monomial = prod([X[i] ** nu[i] for i in xrange(d)])
+ val = F.monomial_coefficient(monomial)
+ if numerical:
+ val = val.n(digits=numerical)
+ coeffs[tuple(nu)] = val
+ return coeffs
+
+
+ def relative_error(self, approx, alpha, interval, exp_scale=Integer(1),
+ digits=10):
+ r"""
+ Return 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.denominator_ring.ngens()``
+ - ``interval`` -- a list of positive integers
+ - ``exp_scale`` -- (optional; default: 1) a number
+
+ OUTPUT:
+
+ A list of tuples with properties described below.
+
+ This outputs a list whose entries are a tuple
+ ``(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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ 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
+ [((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])]
+ """
+ from sage.modules.free_module_element import vector
+
+ if not isinstance(approx, (list, tuple)):
+ approx = [approx]
+ if approx[0].variables():
+ av = approx[0].variables()[0]
+ else:
+ av = Integer(1)
+
+ #print "Calculating errors table in the form"
+ #print "exponent, scaled Maclaurin coefficient, scaled asymptotic values, relative errors..."
+
+ # Get Maclaurin coefficients of self.
+ alpha = vector(alpha)
+ multi_indices = [r * alpha for r in interval]
+ mac = self.maclaurin_coefficients(multi_indices, numerical=digits)
+ #mac = self.old_maclaurin_coefficients(alpha, max(interval))
+ mac_approx = {}
+ stats = []
+ for r in interval:
+ exp_s_r = exp_scale ** r
+ beta = tuple(r * alpha)
+ mac[beta] = (mac[beta] / exp_s_r).n(digits=digits)
+ mac_approx[beta] = [(f.subs({av: r}) / exp_s_r).n(digits=digits)
+ for f in approx]
+ stats_row = [beta, mac[beta], mac_approx[beta]]
+ if mac[beta] == 0:
+ stats_row.extend([None for a in mac_approx[beta]])
+ else:
+ stats_row.append([(mac[beta] - a) / mac[beta]
+ for a in mac_approx[beta]])
+ stats.append(tuple(stats_row))
+ return stats
+
+
+ def _add_(left, right):
+ r"""
+ Returns the sum of ``left`` with ``right``.
+
+ INPUT:
+
+ - ``left`` -- the left summand (i.e. self).
+
+ - ``right`` -- the right summand.
+
+ OUTPUT:
+
+ The sum as a new element.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: df = (x, 1), (y, 1), (x*y + 1, 1)
+ sage: f = FFPD(2, df)
+ sage: g = FFPD(2*x*y, df)
+ sage: f + g
+ (2, [(y, 1), (x, 1)])
+ """
+ return FractionWithFactoredDenominatorSum([left, right]).sum()
+
+
+ def _mul_(left, right):
+ r"""
+ Returns the product of ``left`` with ``right``.
+
+ INPUT:
+
+ - ``left`` -- the left factor (i.e. self).
+
+ - ``right`` -- the right factor.
+
+ OUTPUT:
+
+ The product as a new element.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(2, [(x, 1), (x*y + 1, 1), (x*y^2 + 1, 1)])
+ sage: g = FFPD(2*x*y, [(y, 1), (x*y + 1, 1), (x^2*y + 1, 1)])
+ sage: f * g
+ (4, [(x*y + 1, 1), (x*y + 1, 1), (x*y^2 + 1, 1), (x^2*y + 1, 1)])
+ """
+ numer = left.numerator() * right.numerator()
+ df = left.denominator_factored() + right.denominator_factored()
+ return left.parent()(numer, df)
+
+
+def is_FractionWithFactoredDenominatorRing(F):
+ r"""
+ Returns if input is an instance of
+ FractionWithFactoredDenominatorRing.
+
+ INPUT:
+
+ - ``F`` -- the object to test.
+
+ OUTPUT:
+
+ ``True`` or ``False``.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, is_FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R); R
+ Multivariate Polynomial Ring in x, y over Rational Field
+ sage: is_FractionWithFactoredDenominatorRing(FFPD)
+ True
+ """
+ return isinstance(F, FractionWithFactoredDenominatorRing)
+
+
+class FractionWithFactoredDenominatorRing(
+ sage.structure.unique_representation.UniqueRepresentation,
+ sage.rings.ring.Ring):
+ r"""
+ This is the ring of fractions with factored denominator.
+
+ INPUT:
+
+ - ``base`` -- the base ring (a polynomial ring).
+
+ - ``category`` -- (default: ``None``) the category.
+
+ See also the corresponding element :class:`FractionWithFactoredDenominator`.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: df = [x, 1], [y, 1], [x*y+1, 1]
+ sage: f = FFPD(x, df) # indirect doctest
+ sage: f
+ (1, [(y, 1), (x*y + 1, 1)])
+
+ AUTHORS:
+
+ - Daniel Krenn (2014-12-01)
+ """
+ def __init__(self, denominator_ring, numerator_ring=None, category=None):
+ r"""
+ See :class:`FractionWithFactoredDenominatorRing` for details.
+
+ INPUT:
+
+ - ``base`` -- the base ring (a polynomial ring).
+
+ - ``category`` -- (default: ``None``) the category.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: P.<X, Y> = ZZ[]
+ sage: FractionWithFactoredDenominatorRing(P)
+ Ring of fractions with factored denominator
+ over Multivariate Polynomial Ring in X, Y over Integer Ring
+ """
+ if numerator_ring is None:
+ from sage.symbolic.ring import SR
+ numerator_ring = denominator_ring
+ if not numerator_ring.has_coerce_map_from(denominator_ring):
+ raise ValueError('Numerator ring {} has no coercion map from the '
+ 'denominator ring {}.'.format(
+ numerator_ring, denominator_ring))
+
+ self.numerator_ring = numerator_ring
+ self.denominator_ring = denominator_ring
+ super(FractionWithFactoredDenominatorRing, self).__init__(
+ denominator_ring, category=category or sage.categories.rings.Rings())
+
+
+ Element = FractionWithFactoredDenominator
+
+
+ def _repr_(self):
+ r"""
+ Returns a representation.
+
+ OUTPUT:
+
+ A string.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: P.<X, Y> = ZZ[]
+ sage: FractionWithFactoredDenominatorRing(P) # indirect doctest
+ Ring of fractions with factored denominator
+ over Multivariate Polynomial Ring in X, Y over Integer Ring
+ """
+ return ("Ring of fractions with factored denominator "
+ "over %s" % repr(self.base()))
+
+
+ def base_ring(self):
+ r"""
+ Returns the base ring.
+
+ OUTPUT:
+
+ A ring.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: P.<X, Y> = ZZ[]
+ sage: F = FractionWithFactoredDenominatorRing(P); F
+ Ring of fractions with factored denominator
+ over Multivariate Polynomial Ring in X, Y over Integer Ring
+ sage: F.base_ring()
+ Integer Ring
+ sage: F.base()
+ Multivariate Polynomial Ring in X, Y over Integer Ring
+ """
+ return self.base().base_ring()
+
+
+ from sage.misc.decorators import rename_keyword
+ @rename_keyword(deprecation=10519, reduce_='reduce')
+ def _element_constructor_(self, *args, **kwargs):
+ r"""
+ Returns an element of this ring.
+
+ See :class:`FractionWithFactoredDenominator` for details.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: df = [x, 1], [y, 1], [x*y+1, 1]
+ sage: f = FFPD(x, df) # indirect doctest
+ sage: f
+ (1, [(y, 1), (x*y + 1, 1)])
+ """
+ from sage.symbolic.ring import SR
+
+ R = self.base()
+ Q = R.fraction_field()
+
+ # process deprecated keyword arguments
+ hasn = kwargs.has_key('numerator')
+ hasdf = kwargs.has_key('denominator_factored')
+ if hasn:
+ from sage.misc.superseded import deprecation
+ deprecation(10519, "Keyword argument 'numerator' "
+ "is deprecated. "
+ "Ignoring non-keyword argumets (if any). "
+ "Specify numerator and factored denominator "
+ "as first and second argument, i.e., use "
+ "something like FFPD(n, df).")
+ if hasdf:
+ from sage.misc.superseded import deprecation
+ deprecation(10519, "Keyword argument 'denominator_factored' "
+ "is deprecated. "
+ "Ignoring non-keyword argumets (if any). "
+ "Specify numerator and factored denominator "
+ "as first and second argument, i.e., use "
+ "something like FFPD(n, df).")
+ if hasn or hasdf:
+ args = [kwargs.pop('numerator') if hasn else R(0),
+ kwargs.pop('denominator_factored') if hasdf else []]
+
+ hasq = kwargs.has_key('quotient')
+ if hasq:
+ from sage.misc.superseded import deprecation
+ deprecation(10519, "Keyword argument 'quotient' "
+ "is deprecated. "
+ "Ignoring non-keyword argumets (if any). "
+ "Specify numerator and factored denominator "
+ "as first and second argument, i.e., use "
+ "something like FFPD(q).")
+ args = [kwargs.pop('quotient')]
+
+ if (hasn or hasdf) and hasq:
+ raise ValueError('Parameters ambiguous.')
+
+ # process keyword arguments
+ reduce = kwargs.pop('reduce', None)
+
+ if kwargs:
+ raise ValueError('Unknown keyword arguments '
+ '%s given' % (kwargs.keys(),))
+
+ # process arguments
+ if len(args) > 2:
+ raise ValueError('Too many arguments given.')
+
+ elif not args:
+ raise ValueError('No argument given. '
+ 'We are in serious troubles...')
+
+ # At this point we have one or two input arguments.
+
+ x = args[0]
+ try:
+ P = x.parent()
+ except AttributeError:
+ P = None
+
+ denominator_factored = None # init
+ reduce_default = True
+
+ if len(args) == 2:
+ numerator, denominator_factored = args
+ if numerator is None:
+ numerator = R(0)
+ if denominator_factored is None:
+ denominator_factored = []
+
+ from sage.rings.semirings.non_negative_integer_semiring import NN
+ try:
+ denominator_factored = sorted(
+ (R(d[0]), NN(d[1])) for d in denominator_factored)
+ except TypeError:
+ raise TypeError('Factored denominator is not well-formed '
+ 'or of wrong type.')
+
+ # From now on we only have one input arguement;
+ # it's called x and has parent P.
+
+ elif is_FractionWithFactoredDenominatorRing(P):
+ numerator = x._numerator
+ denominator_factored = self._denominator_factored
+ reduce_default = False
+
+ elif P == SR:
+ numerator = x.numerator()
+ denominator = x.denominator()
+ reduce_default = False
+
+ elif x in R:
+ numerator = R(x)
+ denominator_factored = []
+
+ elif x in Q:
+ quotient = Q(x)
+ numerator = quotient.numerator()
+ denominator = quotient.denominator()
+ reduce_default = False
+
+ elif hasattr(x, 'numerator') and hasattr(x, 'denominator'):
+ numerator = x.numerator()
+ denominator = x.denominator()
+ reduce_default = False
+
+ else:
+ raise TypeError('Element %s is not contained '
+ 'in %s.' % (x, self))
+
+ if reduce is None:
+ reduce = reduce_default
+
+ if denominator_factored is None:
+ if denominator not in R:
+ raise TypeError('Extracted denominator %s is not '
+ 'in %s.' % (denominator, self))
+ p = numerator
+ q = R(denominator)
+
+ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
+ from sage.rings.polynomial.multi_polynomial_ring_generic import is_MPolynomialRing
+ if is_PolynomialRing(R) or is_MPolynomialRing(R):
+ if not R(q).is_unit():
+ # Factor denominator
+ try:
+ df = q.factor()
+ except NotImplementedError:
+ # Singular's factor() needs 'proof=False'.
+ df = q.factor(proof=False)
+ numerator = p / df.unit()
+ df = sorted(tuple(t) for t in df) # sort for consitency
+ denominator_factored = df
+ else:
+ # At this point, denominator could not be factored.
+ numerator = p / q
+ denominator_factored = []
+
+ return self.element_class(self,
+ numerator=numerator,
+ denominator_factored=denominator_factored,
+ reduce=reduce)
+
+
+ def _coerce_map_from_(self, P):
+ r"""
+ Checks if there is a coercion from the given parent.
+
+ INPUT:
+
+ - ``P`` -- a parent.
+
+ OUTPUT:
+
+ ``True`` if there is a coercion, otherwise ``False`` or
+ ``None``.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: Q = QQ['x,y']
+ sage: FFPD_QQ = FractionWithFactoredDenominatorRing(Q)
+ sage: FFPD_QQ.has_coerce_map_from(Q)
+ True
+ sage: FFPD_QQ.has_coerce_map_from(QQ)
+ True
+ sage: FFPD_QQ.has_coerce_map_from(ZZ)
+ True
+ sage: FFPD_QQ.has_coerce_map_from(Q.fraction_field())
+ True
+ sage: Z = ZZ['x,y']
+ sage: FFPD_ZZ = FractionWithFactoredDenominatorRing(Z)
+ sage: FFPD_ZZ.has_coerce_map_from(FFPD_QQ)
+ False
+ sage: FFPD_QQ.has_coerce_map_from(FFPD_ZZ)
+ True
+ sage: FFPD_ZZ.has_coerce_map_from(QQ)
+ False
+ sage: FFPD_ZZ.has_coerce_map_from(Z.fraction_field())
+ True
+ sage: FFPD_ZZ.has_coerce_map_from(Q.fraction_field())
+ False
+ sage: FFPD_QQ.has_coerce_map_from(Z.fraction_field())
+ True
+ """
+ if is_FractionWithFactoredDenominatorRing(P):
+ if self.base().has_coerce_map_from(P.base()):
+ return True
+
+ from sage.rings.fraction_field import is_FractionField
+ if is_FractionField(P):
+ B = P.base()
+ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
+ from sage.rings.polynomial.multi_polynomial_ring_generic import is_MPolynomialRing
+ if is_PolynomialRing(B) or is_MPolynomialRing(B):
+ if self.base().has_coerce_map_from(B):
+ return True
+
+ if self.base().has_coerce_map_from(P):
+ return True
+
+
+ def _an_element_(self):
+ r"""
+ Returns an element.
+
+ OUTPUT:
+
+ An element.
+
+ TESTS::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: FFPD.an_element() # indirect doctest
+ (42, [(x, 3)])
+ """
+ from sage.rings.semirings.non_negative_integer_semiring import NN
+ return self(NN.an_element(), [(self.base().an_element(), NN(3))])
+
+
+ @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 :meth:`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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ sage: u = ['a', 'b', 'c', 'd', 'e']
+ sage: s = ['b', 'd']
+ sage: FFDR._permutation_sign(s, u)
+ -1
+ sage: s = ['d', 'b']
+ sage: FFDR._permutation_sign(s, u)
+ 1
+ """
+ from sage.combinat.permutation import Permutation
+ from sage.sets.set import Set
+
+ # 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()
+
+
+ @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`` -- (default: ``False``) boolean; set to ``True`` to
+ simplify the result
+
+ 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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ sage: var('x, y, z')
+ (x, y, z)
+ sage: a = {x:1}
+ sage: b = {y:2}
+ sage: c = {z:3}
+ sage: FFDR.subs_all(x + y + z, a)
+ y + z + 1
+ sage: FFDR.subs_all(x + y + z, [c, a])
+ y + 4
+ sage: FFDR.subs_all([x + y + z, y^2], b)
+ [x + z + 2, 4]
+ sage: FFDR.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: FFDR.subs_all(a, b)
+ {'bar': -1, 'foo': 5}
+ """
+ singleton = False
+ if not isinstance(f, (list, tuple)):
+ f = [f]
+ singleton = True
+ if not isinstance(sub, (list, tuple)):
+ 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)]
+ return g[Integer(0)].simplify()
+
+ if singleton and not simplify:
+ return g[Integer(0)]
+
+ if not singleton and simplify:
+ G = []
+ for gg in g:
+ if isinstance(gg, dict):
+ G.append(gg)
+ else:
+ G.append(gg.simplify())
+ return G
+
+ 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`.
+
+ 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, ..., sr:deriv_r}``.
+
+ Here ``s_1, ..., 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` if and only if ``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, ..., 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::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ sage: f = function('f')(x)
+ sage: dd = FFDR._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 = FFDR._diff_all(f, [x], 3, sub=d1)
+ sage: dd[(x, x, x)]
+ 24*x
+
+ sage: dd = FFDR._diff_all(f, [x], 3, sub=d1, rekey=f)
+ sage: dd[diff(f, x, 3)]
+ 24*x
+
+ sage: a = {x:1}
+ sage: dd = FFDR._diff_all(f, [x], 3, sub=d1, rekey=f, sub_final=a)
+ sage: dd[diff(f, x, 3)]
+ 24
+
+ ::
+
+ sage: X = var('x, y, z')
+ sage: f = function('f')(*X)
+ sage: dd = FFDR._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 = FFDR._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 = FFDR._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)
+ """
+ from sage.calculus.functional import diff
+ from sage.combinat.tuple import UnorderedTuples
+
+ 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)] = FractionWithFactoredDenominatorRing.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] = FractionWithFactoredDenominatorRing.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 = FractionWithFactoredDenominatorRing.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 = FractionWithFactoredDenominatorRing.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] = FractionWithFactoredDenominatorRing.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.
+
+ 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.
+
+ The output is
+ a dictionary whose keys are natural number tuples of the form
+ `(j, k, l)`, where `l \leq 2k`, `j \leq r-1`, and `j+k \leq 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])`.
+
+ .. NOTE::
+
+ For internal use by :meth:`asymptotics_smooth()` and
+ :meth:`asymptotics_multiple()`.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ 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 = FFDR._diff_op(A, B, AB_derivs, T, M, 1, 2)
+ sage: sorted(DD.keys())
+ [(0, 0, 0), (0, 1, 0), (0, 1, 1), (0, 1, 2)]
+ sage: len(DD[(0, 1, 2)])
+ 246
+ """
+ from itertools import product
+ from sage.calculus.functional import diff
+ from sage.combinat.tuple import UnorderedTuples
+ from sage.misc.mrange import xmrange
+
+ if not isinstance(A, list):
+ A = [A]
+
+ # First, compute the necessary product derivatives of A and B.
+ product_derivs = {}
+ for j, k in xmrange([r, N], tuple):
+ if j + k < N:
+ for l in xrange(2 * k + 1):
+ for s in UnorderedTuples(V, 2 * (k + l)):
+ DF = diff(A[j] * B ** l, s).subs(AB_derivs)
+ product_derivs[tuple([j, k, l] + s)] = DF
+
+ # Second, compute DD^(k+l)(A[j]*B^l)(p) and store values in dictionary.
+ DD = {}
+ rows = M.nrows()
+ for j, k in xmrange([r, N], tuple):
+ 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 xmrange([rows, rows], tuple) if b <= a]
+ P = product(S, repeat=k + l)
+ diffo = Integer(0)
+ for t in P:
+ idx = (j, k, l) + FractionWithFactoredDenominatorRing._diff_seq(V, t)
+ if product_derivs[idx] != Integer(0):
+ MM = Integer(1)
+ for (a, b) in t:
+ MM *= M[a][b]
+ if a != b:
+ MM *= Integer(2)
+ diffo += MM * product_derivs[idx]
+ 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``.
+
+ 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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ sage: V = list(var('x, t, z'))
+ sage: FFDR._diff_seq(V,([0, 1],[0, 2, 1],[0, 0]))
+ (x, x, x, x, t, t, z)
+
+ .. NOTE::
+
+ This function is for internal use by :meth:`diff_op()`.
+ """
+ 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`` -- a symbolic variable
+ - ``a`` -- a complex number
+ - ``v, N`` -- natural numbers
+
+ OUTPUT:
+
+ A dictionary.
+
+ The output is
+ a dictionary whose keys are natural number pairs of the form `(k, l)`,
+ where `k < N` and `l \leq 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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ sage: A = function('A')(x)
+ sage: B = function('B')(x)
+ sage: AB_derivs = {}
+ sage: sorted(FFDR._diff_op_simple(A, B, AB_derivs, x, 3, 2, 2).items())
+ [((0, 0), A(x)),
+ ((1, 0), 1/2*I*2^(2/3)*D[0](A)(x)),
+ ((1, 1), 1/4*2^(2/3)*(B(x)*D[0, 0, 0, 0](A)(x)
+ + 4*D[0, 0, 0](A)(x)*D[0](B)(x) + 6*D[0, 0](A)(x)*D[0, 0](B)(x)
+ + 4*D[0](A)(x)*D[0, 0, 0](B)(x) + A(x)*D[0, 0, 0, 0](B)(x)))]
+ """
+ from sage.calculus.functional import diff
+ from sage.functions.other import sqrt
+
+ 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`.
+
+ 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`.
+
+ 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`.
+
+ .. NOTE::
+
+ For internal use by :meth:`asymptotics_multiple()`.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ 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 = FFDR._diff_prod(fd, u, g, [x], [1, 2], [], ud, atc)
+ sage: dd[diff(u, x, 2)(x=2)]
+ 22/9
+ """
+ from sage.calculus.functional import diff
+ from sage.calculus.var import var
+ from sage.combinat.tuple import UnorderedTuples
+ from sage.symbolic.relation import solve
+
+ 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(
+ FractionWithFactoredDenominatorRing.subs_all(D, sol[Integer(0)]))
+ return uderivs
+
+
+ @staticmethod
+ def direction(v, coordinate=None):
+ r"""
+ Return ``[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.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing as FFDR
+ sage: FFDR.direction([2, 3, 5])
+ (2/5, 3/5, 1)
+ sage: FFDR.direction([2, 3, 5], 0)
+ (1, 3/2, 5/2)
+ """
+ if coordinate is None:
+ coordinate = len(v) - 1
+ return tuple([vv / v[coordinate] for vv in v])
+
+
+ @staticmethod
+ def _coerce_point_(R, p):
+ r"""
+ Coerce the keys of the dictionary ``p`` into the ring ``R``.
+
+ Assume that it is possible.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD()
+ sage: p = {SR(x): 1, SR(y): 7/8}
+ sage: p
+ {y: 7/8, x: 1}
+ sage: for k in sorted(p.keys()):
+ ....: print k, k.parent()
+ y Symbolic Ring
+ x Symbolic Ring
+ sage: q = FractionWithFactoredDenominatorRing._coerce_point_(R, p)
+ sage: q
+ {y: 7/8, x: 1}
+ sage: for k in sorted(q.keys()):
+ ....: print k, k.parent()
+ y Multivariate Polynomial Ring in x, y over Rational Field
+ x Multivariate Polynomial Ring in x, y over Rational Field
+ """
+ from sage.symbolic.ring import SR
+
+ 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 FractionWithFactoredDenominatorSum(list):
+ r"""
+ A list representing the sum of :class:`FractionWithFactoredDenominator` objects with distinct
+ denominator factorizations.
+
+ AUTHORS:
+
+ - Alexander Raichev (2012-06-25)
+
+ - Daniel Krenn (2014-12-01)
+ """
+ def __repr__(self):
+ r"""
+ Return a string representation of ``self``.
+
+ OUTPUT:
+
+ A string.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x + y, [(y, 1), (x, 1)])
+ sage: g = FFPD(x**2 + y, [(y, 1), (x, 2)])
+ sage: FractionWithFactoredDenominatorSum([f, g])
+ (x + y, [(y, 1), (x, 1)]) + (x^2 + y, [(y, 1), (x, 2)])
+ """
+ return ' + '.join(repr(r) for r in self)
+
+
+ def __eq__(self, other):
+ r"""
+ Return ``True`` if ``self`` is equal to ``other``.
+
+ OUTPUT:
+
+ A boolean.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x + y, [(y, 1), (x, 1)])
+ sage: g = FFPD(x*(x + y), [(y, 1), (x, 2)])
+ sage: s = FractionWithFactoredDenominatorSum([f]); s
+ (x + y, [(y, 1), (x, 1)])
+ sage: t = FractionWithFactoredDenominatorSum([g]); t
+ (x + y, [(y, 1), (x, 1)])
+ sage: s == t
+ True
+ """
+ return sorted(self) == sorted(other)
+
+
+ def __ne__(self, other):
+ r"""
+ Return ``True`` if ``self`` is not equal to ``other``.
+
+ OUTPUT:
+
+ A boolean.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x + y, [(y, 1), (x, 1)])
+ sage: g = FFPD(x + y, [(y, 1), (x, 2)])
+ sage: s = FractionWithFactoredDenominatorSum([f]); s
+ (x + y, [(y, 1), (x, 1)])
+ sage: t = FractionWithFactoredDenominatorSum([g]); t
+ (x + y, [(y, 1), (x, 2)])
+ sage: s != t
+ True
+ """
+ return not self.__eq__(other)
+
+
+ @property
+ def denominator_ring(self):
+ r"""
+ Return the polynomial ring of the denominators of ``self``.
+
+ OUTPUT:
+
+ A ring or ``None`` if the list is empty.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R)
+ sage: f = FFPD(x + y, [(y, 1), (x, 1)])
+ sage: s = FractionWithFactoredDenominatorSum([f])
+ sage: s.denominator_ring
+ Multivariate Polynomial Ring in x, y over Rational Field
+ sage: g = FFPD(x + y, [])
+ sage: t = FractionWithFactoredDenominatorSum([g])
+ sage: t.denominator_ring
+ Multivariate Polynomial Ring in x, y over Rational Field
+ """
+ for r in self:
+ return r.denominator_ring
+ return None
+
+
+ def whole_and_parts(self):
+ r"""
+ Rewrite ``self`` as a sum
+ of a (possibly zero) polynomial
+ followed by reduced rational expressions.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ Only useful for multivariate decompositions.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: f = x**2 + 3*y + 1/x + 1/y
+ sage: f = FFPD(f); f
+ (x^3*y + 3*x*y^2 + x + y, [(y, 1), (x, 1)])
+ sage: FractionWithFactoredDenominatorSum([f]).whole_and_parts()
+ (x^2 + 3*y, []) + (x + y, [(y, 1), (x, 1)])
+
+ sage: f = cos(x)**2 + 3*y + 1/x + 1/y; f
+ cos(x)^2 + 3*y + 1/x + 1/y
+ sage: G = f.numerator()
+ sage: H = R(f.denominator())
+ sage: f = FFPD(G, H.factor()); f
+ (x*y*cos(x)^2 + 3*x*y^2 + x + y, [(y, 1), (x, 1)])
+ sage: FractionWithFactoredDenominatorSum([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.denominator_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(r.parent()(b, r.denominator_factored(), reduce=False))
+ return FractionWithFactoredDenominatorSum(
+ [r.parent()(whole, ())] + parts) # r.parent() is not the nicest here
+
+
+ def _combine_like_terms_(self):
+ r"""
+ Combine terms in ``self`` with the same denominator.
+ Only useful for multivariate decompositions.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominatorSum`.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: f = FFPD(1/(x * y * (x*y + 1)))
+ sage: g = FFPD(x/(x * y * (x*y + 1)))
+ sage: s = FractionWithFactoredDenominatorSum([f, g, f])
+ sage: t = s._combine_like_terms_()
+ sage: 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: t
+ (1, [(y, 1), (x*y + 1, 1)]) + (2, [(y, 1), (x, 1), (x*y + 1, 1)])
+
+ sage: H = x * y * (x*y + 1)
+ sage: f = FFPD(1, H.factor())
+ sage: g = FFPD(exp(x + y), H.factor())
+ sage: s = FractionWithFactoredDenominatorSum([f, g])
+ sage: 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: 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 = f.parent()(num, temp.denominator_factored())
+ else:
+ # Append temp to new_FFPDs and update temp.
+ new_FFPDs.append(temp)
+ temp = f
+ new_FFPDs.append(temp)
+ return FractionWithFactoredDenominatorSum(new_FFPDs)
+
+
+ def sum(self):
+ r"""
+ Return the sum of the elements in ``self``.
+
+ OUTPUT:
+
+ An instance of :class:`FractionWithFactoredDenominator`.
+
+ EXAMPLES::
+
+ sage: from sage.rings.asymptotic.asymptotics_multivariate_generating_functions import FractionWithFactoredDenominatorRing, FractionWithFactoredDenominatorSum
+ sage: R.<x,y> = PolynomialRing(QQ)
+ sage: FFPD = FractionWithFactoredDenominatorRing(R, SR)
+ sage: df = (x, 1), (y, 1), (x*y + 1, 1)
+ sage: f = FFPD(2, df)
+ sage: g = FFPD(2*x*y, df)
+ sage: FractionWithFactoredDenominatorSum([f, g])
+ (2, [(y, 1), (x, 1), (x*y + 1, 1)]) + (2, [(x*y + 1, 1)])
+ sage: FractionWithFactoredDenominatorSum([f, g]).sum()
+ (2, [(y, 1), (x, 1)])
+
+ sage: f = FFPD(cos(x), [(x, 2)])
+ sage: g = FFPD(cos(y), [(x, 1), (y, 2)])
+ sage: FractionWithFactoredDenominatorSum([f, g])
+ (cos(x), [(x, 2)]) + (cos(y), [(y, 2), (x, 1)])
+ sage: FractionWithFactoredDenominatorSum([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.denominator_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 FractionWithFactoredDenominatorRing(numer.parent())(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 FractionWithFactoredDenominatorRing(numer.parent())(numer, df, reduce=False)
+