summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorRelease Manager <release@sagemath.org>2016-02-16 18:03:23 +0100
committerVolker Braun <vbraun.name@gmail.com>2016-02-16 18:03:23 +0100
commit28db6f605cf167d3c8a59fdef21e1ad59307348d (patch)
tree5cb940f9f18515f66393e25a158a21f8a35751f4
parentTrac #20057: Add iterator to DisjointSet class (diff)
parentTrac #10519: trivial ReST error (diff)
Trac #10519: analytic combinatorics: new code for computing asymptotics for multivariate generating functions
This code is a collection of functions designed to compute asymptotics of Maclaurin coefficients of certain classes of multivariate generating functions. The main function asymptotics() returns the first `N` terms of the asymptotic expansion of the Maclaurin coefficients `F_{n\alpha}` of the multivariate meromorphic function `F=G/H` as `n\to\infty`. It assumes that `F` is holomorphic in a neighborhood of the origin, that `H` is a polynomial, and that asymptotics in the direction of `\alpha` (a tuple of positive integers) are controlled by smooth or multiple points. URL: http://trac.sagemath.org/10519 Reported by: araichev Ticket author(s): Daniel Krenn, Alex Raichev Reviewer(s): Daniel Krenn, David Loeffler, Travis Scrimshaw
-rw-r--r--src/doc/en/reference/asymptotic/index.rst1
-rw-r--r--src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py4388
2 files changed, 4389 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 00000000..1ab30d0
--- /dev/null
+++ b/src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py
@@ -0,0 +1,4388 @@
+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*. Online Journal of Analytic Combinatorics.
+ Issue 6, (2011). :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 == [{SR(x): 1, SR(y): 1}]
+ True
+ 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')
+ sage: 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/
+#*****************************************************************************
+
+
+from functools import total_ordering
+from itertools import combinations_with_replacement
+from sage.structure.element import RingElement
+from sage.structure.unique_representation import UniqueRepresentation
+from sage.rings.ring import Ring
+from sage.calculus.var import var
+from sage.calculus.functional import diff
+from sage.symbolic.ring import SR
+from sage.misc.misc_c import prod
+from sage.rings.integer import Integer
+from sage.rings.integer_ring import ZZ
+from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+from sage.categories.rings import Rings
+
+
+@total_ordering
+class FractionWithFactoredDenominator(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"""
+ Initialize ``self``.
+
+ 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