Skip to content
This repository has been archived by the owner on Jan 30, 2023. It is now read-only.

Commit

Permalink
Merge branch 'u/dkrenn/asy/generators-binomial' of trac.sagemath.org:…
Browse files Browse the repository at this point in the history 
…sage into t/19898/asy/harmonic-number

* 'u/dkrenn/asy/generators-binomial' of trac.sagemath.org:sage:
  Trac #19510: major rewrite of binomial kn over n
  Trac #19510: separate calculation of negative powers in log_Stirling
  Trac #19946: fix typo
  Trac #19946: link from general doc to detailed explaination
  Trac #19946: rewrite and improve explaination of 1b62954
  Trac #19946: additional doctest
  Trac #19946: reviewer commit: ReSt error
  Trac #19961: document rpow
  Trac #19946: add doctests to document behavior
  Trac #19946: fix _pushout_ for cartesian product of growth groups
  Trac #19510: state constant factors/summands in Stirling and binomial
  Trac #19510: modify constant factor: factor sqrt(k-1) is removed, too.
  • Loading branch information
@dkrenn
dkrenn committed Jan 27, 2016
2 parents 27e9e69 + 10cd04b commit 209329d
Show file tree
Hide file tree
Showing 5 changed files with 208 additions and 91 deletions.
201 changes: 114 additions & 87 deletions src/sage/rings/asymptotic/asymptotic_expansion_generators.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ def Stirling(var, precision=None, skip_constant_factor=False):
the default precision of the asymptotic ring is used.
- ``skip_constant_factor`` -- (default: ``False``) a
boolean. If set, then the constant factor is left out.
boolean. If set, then the constant factor `\sqrt{2\pi}` is left out.
As a consequence, the coefficient ring of the output changes
from ``Symbolic Constants Subring`` (if ``False``) to
``Rational Field`` (if ``True``).
Expand Down Expand Up @@ -193,7 +193,7 @@ def log_Stirling(var, precision=None, skip_constant_summand=False):
the default precision of the asymptotic ring is used.
- ``skip_constant_summand`` -- (default: ``False``) a
boolean. If set, then the constant summand is left out.
boolean. If set, then the constant summand `\log(2\pi)/2` is left out.
As a consequence, the coefficient ring of the output changes
from ``Symbolic Constants Subring`` (if ``False``) to
``Rational Field`` (if ``True``).
Expand Down Expand Up @@ -279,10 +279,8 @@ def log_Stirling(var, precision=None, skip_constant_summand=False):
if precision >= 4 and not skip_constant_summand:
result += log(2*coefficient_ring('pi')) / 2

from sage.misc.misc import srange
from sage.arith.all import bernoulli
for k in srange(2, 2*precision - 6, 2):
result += bernoulli(k) / k / (k-1) / n**(k-1)
result += AsymptoticExpansionGenerators._log_StirlingNegativePowers_(
var, precision - 4)

if precision < 1:
result += (n * log(n)).O()
Expand All @@ -292,12 +290,60 @@ def log_Stirling(var, precision=None, skip_constant_summand=False):
result += log(n).O()
elif precision == 3:
result += A(1).O()
else:
result += (1 / n**(2*precision - 7)).O()

return result


@staticmethod
def _log_StirlingNegativePowers_(var, precision):
r"""
Helper function to calculate the logarithm of Stirling's approximation
formula from the negative powers of ``var`` on, i.e., it skips the
summands `n \log n - n + (\log n)/2 + \log(2\pi)/2`.
INPUT:
- ``var`` -- a string for the variable name.
- ``precision`` -- an integer specifying the number of exact summands.
If this is negative, then the result is `0`.
OUTPUT:
An asymptotic expansion.
TESTS::
sage: asymptotic_expansions._log_StirlingNegativePowers_(
....: 'm', precision=-1)
0
sage: asymptotic_expansions._log_StirlingNegativePowers_(
....: 'm', precision=0)
O(m^(-1))
sage: asymptotic_expansions._log_StirlingNegativePowers_(
....: 'm', precision=3)
1/12*m^(-1) - 1/360*m^(-3) + 1/1260*m^(-5) + O(m^(-7))
sage: _.parent()
Asymptotic Ring <m^ZZ> over Rational Field
"""
from asymptotic_ring import AsymptoticRing
from sage.rings.rational_field import QQ

A = AsymptoticRing(growth_group='{n}^ZZ'.format(n=var),
coefficient_ring=QQ)
if precision < 0:
return A.zero()
n = A.gen()

from sage.arith.all import bernoulli
from sage.misc.misc import srange

result = sum((bernoulli(k) / k / (k-1) / n**(k-1)
for k in srange(2, 2*precision + 2, 2)),
A.zero())
return result + (1 / n**(2*precision + 1)).O()


@staticmethod
def HarmonicNumber(var, precision=None, skip_constant_summand=False):
r"""
Expand Down Expand Up @@ -400,8 +446,7 @@ def HarmonicNumber(var, precision=None, skip_constant_summand=False):


@staticmethod
def Binomial_kn_over_n(var, k, precision=None, skip_constant_factor=False,
algorithm=None):
def Binomial_kn_over_n(var, k, precision=None, skip_constant_factor=False):
r"""
Return the asymptotic expansion of the binomial coefficient
`kn` choose `n`.
Expand All @@ -410,126 +455,108 @@ def Binomial_kn_over_n(var, k, precision=None, skip_constant_factor=False,
- ``var`` -- a string for the variable name.
- ``k`` -- a number.
- ``k`` -- a number or symbolic constant.
- ``precision`` -- (default: ``None``) an integer. If ``None``, then
the default precision of the asymptotic ring is used.
- ``skip_constant_factor`` -- (default: ``False``) a
boolean. If set, then the constant summand is left out.
boolean. If set, then the constant factor `\sqrt{k/(2\pi(k-1))}`
is left out.
As a consequence, the coefficient ring of the output changes
from ``Symbolic Constants Subring`` (if ``False``) to
``Rational Field`` (if ``True``).
- ``algorithm`` -- either ``'direct'`` or ``'log'`` (default).
OUTPUT:
An asymptotic expansion.
EXAMPLES::
sage: asymptotic_expansions.Binomial_kn_over_n('n', k=2, precision=5)
1/sqrt(pi)*(e^n)^(2*log(2))*n^(-1/2)
- 1/8/sqrt(pi)*(e^n)^(2*log(2))*n^(-3/2)
+ 1/128/sqrt(pi)*(e^n)^(2*log(2))*n^(-5/2)
+ O((e^n)^(2*log(2))*n^(-7/2))
sage: asymptotic_expansions.Binomial_kn_over_n('n', k=2, precision=3)
1/sqrt(pi)*4^n*n^(-1/2)
- 1/8/sqrt(pi)*4^n*n^(-3/2)
+ 1/128/sqrt(pi)*4^n*n^(-5/2)
+ O(4^n*n^(-7/2))
sage: _.parent()
Asymptotic Ring <(e^(n*log(n)))^(Symbolic Constants Subring)
* (e^n)^(Symbolic Constants Subring)
* n^(Symbolic Constants Subring)
* log(n)^(Symbolic Constants Subring)>
over Symbolic Constants Subring
Asymptotic Ring <QQ^n * n^QQ> over Symbolic Constants Subring
::
sage: asymptotic_expansions.Binomial_kn_over_n('n', k=3, precision=5)
1/2*sqrt(3)/sqrt(pi)*(e^n)^(3*log(3) - 2*log(2))*n^(-1/2)
- 7/144*sqrt(3)/sqrt(pi)*(e^n)^(3*log(3)
- 2*log(2))*n^(-3/2) + 49/20736*sqrt(3)/sqrt(pi)*(e^n)^(3*log(3)
- 2*log(2))*n^(-5/2)
+ O((e^n)^(3*log(3) - 2*log(2))*n^(-7/2))
sage: asymptotic_expansions.Binomial_kn_over_n('n', k=3, precision=3)
1/2*sqrt(3)/sqrt(pi)*(27/4)^n*n^(-1/2)
- 7/144*sqrt(3)/sqrt(pi)*(27/4)^n*n^(-3/2)
+ 49/20736*sqrt(3)/sqrt(pi)*(27/4)^n*n^(-5/2)
+ O((27/4)^n*n^(-7/2))
::
sage: asymptotic_expansions.Binomial_kn_over_n('n', k=7/5, precision=5)
1/2*sqrt(7)/sqrt(pi)*(e^n)^(7/5*log(7/5) - 2/5*log(2/5))*n^(-1/2)
- 13/112*sqrt(7)/sqrt(pi)*(e^n)^(7/5*log(7/5) - 2/5*log(2/5))*n^(-3/2)
+ 169/12544*sqrt(7)/sqrt(pi)*(e^n)^(7/5*log(7/5) - 2/5*log(2/5))*n^(-5/2)
+ O((e^n)^(7/5*log(7/5) - 2/5*log(2/5))*n^(-7/2))
sage: asymptotic_expansions.Binomial_kn_over_n('n', k=7/5, precision=3)
1/2*sqrt(7)/sqrt(pi)*(7/10*7^(2/5)*2^(3/5))^n*n^(-1/2)
- 13/112*sqrt(7)/sqrt(pi)*(7/10*7^(2/5)*2^(3/5))^n*n^(-3/2)
+ 169/12544*sqrt(7)/sqrt(pi)*(7/10*7^(2/5)*2^(3/5))^n*n^(-5/2)
+ O((7/10*7^(2/5)*2^(3/5))^n*n^(-7/2))
sage: _.parent()
Asymptotic Ring <(Symbolic Constants Subring)^n * n^QQ>
over Symbolic Constants Subring
TESTS::
sage: asymptotic_expansions.Binomial_kn_over_n(
....: 'n', k=5, precision=5, skip_constant_factor=True)
1/2*(e^n)^(5*log(5)
- 4*log(4))*n^(-1/2)
- 7/160*(e^n)^(5*log(5) - 4*log(4))*n^(-3/2)
+ 49/25600*(e^n)^(5*log(5) - 4*log(4))*n^(-5/2)
+ O((e^n)^(5*log(5) - 4*log(4))*n^(-7/2))
....: 'n', k=5, precision=3, skip_constant_factor=True)
(3125/256)^n*n^(-1/2)
- 7/80*(3125/256)^n*n^(-3/2)
+ 49/12800*(3125/256)^n*n^(-5/2)
+ O((3125/256)^n*n^(-7/2))
sage: _.parent()
Asymptotic Ring <(e^(n*log(n)))^(Symbolic Constants Subring)
* (e^n)^(Symbolic Constants Subring)
* n^(Symbolic Constants Subring)
* log(n)^(Symbolic Constants Subring)>
over Rational Field
Asymptotic Ring <QQ^n * n^QQ> over Rational Field
sage: asymptotic_expansions.Binomial_kn_over_n(
....: 'n', k=4, precision=1, skip_constant_factor=True)
(256/27)^n*n^(-1/2) + O((256/27)^n*n^(-3/2))
::
sage: S = asymptotic_expansions.Stirling('n', precision=5)
sage: n = S.parent().gen()
sage: all( # long time
....: asymptotic_expansions.Binomial_kn_over_n('n', k=k,
....: precision=5, algorithm='direct').has_same_summands(
....: asymptotic_expansions.Binomial_kn_over_n('n', k=k,
....: precision=5, algorithm='log'))
....: for k in (2, 4))
....: SR(asymptotic_expansions.Binomial_kn_over_n(
....: 'n', k=k, precision=3)).canonicalize_radical() ==
....: SR(S.subs(n=k*n) / (S.subs(n=(k-1)*n) * S)).canonicalize_radical()
....: for k in [2, 3, 4])
True
"""
if algorithm is None:
algorithm = 'log' # 'log' seems to be faster, so we use it
# as a default.

from sage.symbolic.ring import SR
SCR = SR.subring(no_variables=True)
try:
k = SCR.coerce(k)
SCR.coerce(k)
except TypeError as e:
from misc import combine_exceptions
raise combine_exceptions(
TypeError('Cannot use k=%s.' % (k,)), e)

if algorithm == 'direct':
Stirling = AsymptoticExpansionGenerators.Stirling(
var, precision=precision, skip_constant_factor=skip_constant_factor)
n = Stirling.parent().gen()
result = Stirling.subs(n=k*n) / \
(Stirling.subs(n=(k-1)*n) * Stirling)

elif algorithm == 'log':
log_Stirling = AsymptoticExpansionGenerators.log_Stirling(
var, precision=precision, skip_constant_summand=True)
n = log_Stirling.parent().gen()
TypeError('Cannot use k={}.'.format(k)), e)

result = log_Stirling.subs(n=k*n) - \
log_Stirling.subs(n=(k-1)*n) - log_Stirling
S = AsymptoticExpansionGenerators._log_StirlingNegativePowers_(
var, precision=max(precision - 2,0))
n = S.parent().gen()
result = (S.subs(n=k*n) - S.subs(n=(k-1)*n) - S).exp()

P = log_Stirling.parent().change_parameter(
growth_group='(e^(n*log(n)))^QQ * (e^n)^QQ * n^QQ * log(n)^QQ',
coefficient_ring=SCR)
from sage.functions.log import exp
result = exp(P.coerce(result))
from sage.rings.rational_field import QQ

if not skip_constant_factor:
result /= (2*SCR('pi')).sqrt()
P = S.parent().change_parameter(
growth_group='QQ^{n} * {n}^QQ'.format(n=var),
coefficient_ring=QQ)
n = P.gen()

result = result.map_coefficients(
lambda c: c.canonicalize_radical())

else:
raise ValueError('Unknown algorithm %s.' % (algorithm,))
b = k**k / (k-1)**(k-1)
if b.parent() is SR:
b = SCR(b).canonicalize_radical()
result *= n.rpow(b)
result *= n**(-QQ(1)/QQ(2))
if not skip_constant_factor:
result *= (k/((k-1)*2*SCR('pi'))).sqrt()

if skip_constant_factor:
from sage.rings.rational_field import QQ
result = result.parent().change_parameter(
coefficient_ring=QQ)(result / SCR(k.sqrt()))
return result


# Easy access to the asymptotic expansions generators from the command line:
asymptotic_expansions = AsymptoticExpansionGenerators()
r"""
Expand Down
50 changes: 50 additions & 0 deletions src/sage/rings/asymptotic/asymptotic_ring.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,36 @@
sage: (1 + 1/z + O(1/z^5))^(1 + 1/z)
1 + z^(-1) + z^(-2) + 1/2*z^(-3) + 1/3*z^(-4) + O(z^(-5))
.. NOTE::
In the asymptotic ring
::
sage: M.<n> = AsymptoticRing(growth_group='QQ^n * n^QQ', coefficient_ring=ZZ)
the operation
::
sage: (1/2)^n
Traceback (most recent call last):
...
ValueError: 1/2 is not in Exact Term Monoid QQ^n * n^QQ
with coefficients in Integer Ring. ...
fails, since the rational `1/2` is not contained in `M`. You can use
::
sage: n.rpow(1/2)
(1/2)^n
instead. (See also the examples in
:meth:`ExactTerm.rpow() <sage.rings.asymptotic.term_monoid.ExactTerm.rpow>`
for a detailed explanation.)
Another way is to use a larger coefficent ring::
sage: M_QQ.<n> = AsymptoticRing(growth_group='QQ^n * n^QQ', coefficient_ring=QQ)
sage: (1/2)^n
(1/2)^n
Multivariate Arithmetic
^^^^^^^^^^^^^^^^^^^^^^^
Expand Down Expand Up @@ -1416,6 +1446,18 @@ def __pow__(self, exponent, precision=None):
ValueError: Cannot take s + t to the exponent s.
> *previous* ValueError: log(s + t) cannot be constructed since
there are several maximal elements s, t.
Check that :trac:`19946` is fixed::
sage: A.<n> = AsymptoticRing('QQ^n * n^QQ', SR)
sage: e = 2^n; e
2^n
sage: e.parent()
Asymptotic Ring <SR^n * n^SR> over Symbolic Ring
sage: e = A(e); e
2^n
sage: e.parent()
Asymptotic Ring <QQ^n * n^QQ> over Symbolic Ring
"""
if not self.summands:
if exponent == 0:
Expand Down Expand Up @@ -1674,6 +1716,14 @@ def rpow(self, base, precision=None):
ArithmeticError: Cannot construct y^x in Growth Group x^ZZ
> *previous* TypeError: unsupported operand parent(s) for '*':
'Growth Group x^ZZ' and 'Growth Group SR^x'
Check that :trac:`19946` is fixed::
sage: A.<n> = AsymptoticRing('QQ^n * n^QQ', SR)
sage: n.rpow(2)
2^n
sage: _.parent()
Asymptotic Ring <QQ^n * n^SR> over Symbolic Ring
"""
if isinstance(base, AsymptoticExpansion):
return base.__pow__(self, precision=precision)
Expand Down
8 changes: 8 additions & 0 deletions src/sage/rings/asymptotic/growth_group.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -873,6 +873,14 @@ def _rpow_(self, base):
sage: x = G('x')
sage: (x * log(x)).rpow(2) # indirect doctest
2^(x*log(x))
::
sage: n = GrowthGroup('QQ^n * n^QQ')('n')
sage: n.rpow(2)
2^n
sage: _.parent()
Growth Group QQ^n * n^QQ
"""
if base == 0:
raise ValueError('%s is not an allowed base for calculating the '
Expand Down
Loading

0 comments on commit 209329d

Please sign in to comment.