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

Commit

Permalink
Trac #19969: unique code for all three cases
Browse files Browse the repository at this point in the history
  • Loading branch information
@cheuberg
cheuberg committed Feb 7, 2016
1 parent 9488698 commit 8938c97
Showing 1 changed file with 82 additions and 149 deletions.
231 changes: 82 additions & 149 deletions src/sage/rings/asymptotic/asymptotic_expansion_generators.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -733,9 +733,10 @@ def SingularityAnalysis(var, zeta=1, alpha=0, beta=0, delta=0,
sage: asymptotic_expansions.SingularityAnalysis(
....: 'n', alpha=0)
0
sage: _.parent()
Asymptotic Ring <n^ZZ> over Rational Field
Traceback (most recent call last):
...
NotImplementedOZero: The result is O(0) which means 0
for sufficiently large n
sage: asymptotic_expansions.SingularityAnalysis(
....: 'n', alpha=-1)
Traceback (most recent call last):
Expand Down Expand Up @@ -780,13 +781,38 @@ def SingularityAnalysis(var, zeta=1, alpha=0, beta=0, delta=0,
1/sqrt(pi)*(1/2)^n*n^(-1/2) - 1/8/sqrt(pi)*(1/2)^n*n^(-3/2)
+ 1/128/sqrt(pi)*(1/2)^n*n^(-5/2) + O((1/2)^n*n^(-7/2))
"""
from itertools import islice, count
from asymptotic_ring import AsymptoticRing
from sage.functions.other import gamma
from growth_group import ExponentialGrowthGroup, \
MonomialGrowthGroup
from sage.categories.cartesian_product import cartesian_product
from sage.functions.other import binomial, gamma
from sage.calculus.calculus import limit
from sage.misc.cachefunc import cached_function
from sage.misc.misc import srange
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.symbolic.ring import SR

SCR = SR.subring(no_variables=True)
s = SR('s')

@cached_function
def inverse_gamma_derivative(alpha, r):
"""
Return value of `r`-th derivative of 1/Gamma
at alpha.
"""
if r == 0:
result = 1/gamma(alpha)
else:
result = limit((1/gamma(s)).diff(s, r), s=alpha)

try:
return SCR(result)
except TypeError:
return result

if isinstance(alpha, int):
alpha = ZZ(alpha)
if isinstance(beta, int):
Expand All @@ -797,162 +823,69 @@ def SingularityAnalysis(var, zeta=1, alpha=0, beta=0, delta=0,
if precision is None:
precision = AsymptoticRing.__default_prec__

SCR = SR.subring(no_variables=True)

if delta != 0:
raise NotImplementedError

elif beta != 0:
from growth_group import ExponentialGrowthGroup, \
MonomialGrowthGroup
from sage.categories.cartesian_product import cartesian_product

if zeta == 1:
group = cartesian_product(
[MonomialGrowthGroup(alpha.parent(), var),
MonomialGrowthGroup(beta.parent(), 'log({})'.format(var))])
else:
group = cartesian_product(
[ExponentialGrowthGroup((1/zeta).parent(), var),
MonomialGrowthGroup(alpha.parent(), var),
MonomialGrowthGroup(beta.parent(), 'log({})'.format(var))])

A = AsymptoticRing(growth_group=group, coefficient_ring=SR,
default_prec=precision)
n = A.gen()

if zeta == 1:
exponential_factor = 1
else:
exponential_factor = n.rpow(1/zeta)

if beta in ZZ and beta > 0:
from itertools import chain, islice, count
from sage.functions.other import binomial
from sage.calculus.calculus import limit
L = _sa_coefficients_lambda_(precision, beta=beta)
s = SR('s')
summands = chain.from_iterable(
iter(
iter(binomial(beta, r) * sum(L[(k, ell)] * (-1)**ell *
limit((1/gamma(s)).diff(s, r),
s=alpha-ell)
for ell in srange(k, 2*k+1)
if (k, ell) in L) *
n**(-k) * n.log()**(-r)
for r in srange(beta+1))
for k in count()))
contribution = sum(islice(summands, precision))
error_term = next(summands)
while error_term.is_zero():
error_term = next(summands)

return exponential_factor * n**(alpha-1) * n.log()**beta * \
(contribution + error_term.O())
else:
from itertools import count, islice
from sage.functions.other import binomial
from sage.calculus.calculus import limit
s = SR('s')
summands = iter(binomial(beta, k) *
limit((1/gamma(s)).diff(s, k), s=alpha) *
n.log() ** (-k) for k in count())

return exponential_factor * n**(alpha-1) * n.log()**beta * \
(sum(islice(summands, precision)) +
(n.log()**(-precision)).O())

elif alpha != 0:

iga = QQ(1) / gamma(alpha)
if iga.parent() is SR:
try:
iga = SCR(iga)
except TypeError:
pass

from growth_group import ExponentialGrowthGroup, \
MonomialGrowthGroup

if zeta == 1:
group = MonomialGrowthGroup(alpha.parent(), var)
else:
from sage.categories.cartesian_product \
import cartesian_product
group = cartesian_product(
[ExponentialGrowthGroup((1/zeta).parent(), var),
MonomialGrowthGroup(alpha.parent(), var)])

A = AsymptoticRing(
growth_group=group,
coefficient_ring=iga.parent(), default_prec=precision)
n = A.gen()

if zeta == 1:
exponential_factor = 1
else:
exponential_factor = n.rpow(1/zeta)
groups = []
if zeta != 1:
groups.append(ExponentialGrowthGroup((1/zeta).parent(), var))

e = _sa_coefficients_e_(precision, alpha)
result = sum(n**(alpha-1-k) * iga * e[k]
for k in srange(precision))
groups.append(MonomialGrowthGroup(alpha.parent(), var))
if beta != 0:
groups.append(MonomialGrowthGroup(beta.parent(), 'log({})'.format(var)))

if alpha in ZZ:
a = ZZ(alpha)
if a > 0 and a <= precision:
return exponential_factor * result
elif a < 0:
assert result.is_zero()
raise NotImplementedOZero(
'The result is O(0) which means 0 for sufficiently '
'large {}'.format(var))

return exponential_factor * (result + (n**(alpha-1-precision)).O())
group = cartesian_product(groups)
A = AsymptoticRing(growth_group=group, coefficient_ring=SCR,
default_prec=precision)
n = A.gen()

if zeta == 1:
exponential_factor = 1
else:
return AsymptoticRing(growth_group='%s^ZZ' % (var,),
coefficient_ring=QQ).zero()


def _sa_coefficients_e_(K, alpha):
r"""
Return the coefficient `e_k` used in singularity analysis.
INPUT:
- ``K`` -- an integer specifying the number of coefficients.
- ``alpha`` -- an object.
OUTPUT:
A tuple of objects.
exponential_factor = n.rpow(1/zeta)

.. SEEALSO::
:meth:`~AsymptoticExpansionGenerators.SingularityAnalysis`
TESTS::
sage: from sage.rings.asymptotic.asymptotic_expansion_generators \
....: import _sa_coefficients_e_
sage: a = SR.var('a')
sage: tuple(c.factor() for c in _sa_coefficients_e_(4, a))
(1,
1/2*(a - 1)*a,
1/24*(3*a - 1)*(a - 1)*(a - 2)*a,
1/48*(a - 1)^2*(a - 2)*(a - 3)*a^2)
"""
from sage.arith.all import falling_factorial
from sage.misc.misc import srange
from sage.rings.integer_ring import ZZ
if beta in ZZ and beta >= 0:
it = ((k, r)
for k in count()
for r in srange(beta+1))
k_max = precision
else:
it = ((0, r)
for r in count())
k_max = 0

L = _sa_coefficients_lambda_(K)
return tuple(sum((-1)**ell * L[(k, ell)] *
falling_factorial(alpha - 1, ell)
for ell in srange(k, 2*k+1) if L.has_key((k, ell)))
for k in srange(K))
if beta != 0:
log_n = n.log()
else:
# avoid construction of log(n)
# because it does not exist in growth group.
log_n = 1

it = reversed(list(islice(it, precision+1)))
L = _sa_coefficients_lambda_(max(1, k_max), beta=beta)
(k, r) = next(it)
result = (n**(-k) * log_n ** (-r)).O()

if alpha in ZZ and beta == 0:
if alpha > 0 and alpha <= precision:
result = A(0)
elif alpha <= 0 and precision > 0:
raise NotImplementedOZero(
'The result is O(0) which means 0 for sufficiently '
'large {}'.format(var))

for (k, r) in it:
result += binomial(beta, r) * \
sum(L[(k, ell)] * (-1)**ell *
inverse_gamma_derivative(alpha-ell, r)
for ell in srange(k, 2*k+1)
if (k, ell) in L) * \
n**(-k) * log_n **(-r)

result *= exponential_factor * n**(alpha-1) * log_n**beta

return result

def _sa_coefficients_lambda_(K, beta=0):
r"""
Expand Down

0 comments on commit 8938c97

Please sign in to comment.