Commits (2)
......@@ -285,7 +285,7 @@ cdef class pAdicCappedAbsoluteElement(CAElement):
mpz_clear(ppow_minus_one)
return infinity
def _log_binary_splitting(self, aprec, mina=0):
def _log_binary_splitting(self, aprec):
r"""
Return ``\log(self)`` for ``self`` equal to 1 in the residue field
......@@ -297,9 +297,6 @@ cdef class pAdicCappedAbsoluteElement(CAElement):
- ``aprec`` -- an integer, the precision to which the result is
correct. ``aprec`` must not exceed the precision cap of the ring over
which this element is defined.
- ``mina`` -- an integer (default: 0), the series will check `n` up to
this valuation (and beyond) to see if they can contribute to the
series.
NOTE::
......
......@@ -338,7 +338,7 @@ cdef class pAdicCappedRelativeElement(CRElement):
mpz_mul(selfvalue.value, self.prime_pow.pow_mpz_t_tmp(self.ordp), self.unit)
return Mod(selfvalue, modulus)
def _log_binary_splitting(self, aprec, mina=0):
def _log_binary_splitting(self, aprec):
r"""
Return ``\log(self)`` for ``self`` equal to 1 in the residue field
......@@ -350,9 +350,6 @@ cdef class pAdicCappedRelativeElement(CRElement):
- ``aprec`` -- an integer, the precision to which the result is
correct. ``aprec`` must not exceed the precision cap of the ring over
which this element is defined.
- ``mina`` -- an integer (default: 0), the series will check `n` up to
this valuation (and beyond) to see if they can contribute to the
series.
NOTE::
......
......@@ -359,7 +359,7 @@ cdef class pAdicFixedModElement(FMElement):
mpz_clear(tmp)
return infinity
def _log_binary_splitting(self, aprec, mina=0):
def _log_binary_splitting(self, aprec):
r"""
Return ``\log(self)`` for ``self`` equal to 1 in the residue field
......@@ -371,9 +371,6 @@ cdef class pAdicFixedModElement(FMElement):
- ``aprec`` -- an integer, the precision to which the result is
correct. ``aprec`` must not exceed the precision cap of the ring over
which this element is defined.
- ``mina`` -- an integer (default: 0), the series will check `n` up to
this valuation (and beyond) to see if they can contribute to the
series.
NOTE::
......
......@@ -28,7 +28,7 @@ from sage.rings.finite_rings.integer_mod import Mod
cdef extern from "sage/rings/padics/transcendantal.c":
cdef void padicexp(mpz_t ans, const mpz_t a, unsigned long p, unsigned long prec, const mpz_t modulo)
cdef void padicexp_Newton(mpz_t ans, const mpz_t a, unsigned long p, unsigned long prec, unsigned long precinit, const mpz_t modulo)
cdef void padiclog(mpz_t ans, const mpz_t a, unsigned long p, unsigned long prec, const mpz_t modulo)
cdef class PowComputer_(PowComputer_base):
"""
......@@ -408,6 +408,7 @@ cdef class pAdicFloatingPointElement(FPElement):
sage: x = R(7*w)
sage: x.exp(algorithm="newton") # indirect doctest
1 + w*7 + (4*w + 2)*7^2 + (w + 6)*7^3 + 5*7^4 + O(7^5)
"""
cdef unsigned long p
cdef unsigned long prec = aprec
......@@ -430,3 +431,80 @@ cdef class pAdicFloatingPointElement(FPElement):
return ans
def _log_binary_splitting(self, aprec):
r"""
Return ``\log(self)`` for ``self`` equal to 1 in the residue field
This is a helper method for :meth:`log`.
It uses a fast binary splitting algorithm.
INPUT:
- ``aprec`` -- an integer, the precision to which the result is
correct. ``aprec`` must not exceed the precision cap of the ring over
which this element is defined.
- ``mina`` -- an integer (default: 0), the series will check `n` up to
this valuation (and beyond) to see if they can contribute to the
series.
NOTE::
The function does not check that its argument ``self`` is
1 in the residue field. If this assumption is not fullfiled
the behaviour of the function is not specified.
ALGORITHM:
1. Raise `u` to the power `p^v` for a suitable `v` in order
to make it closer to 1. (`v` is chosen such that `p^v` is
close to the precision.)
2. Write
.. MATH::
u^{p-1} = \prod_{i=1}^\infty (1 - a_i p^{(v+1)*2^i})
with `0 \leq a_i < p^{(v+1)*2^i}` and compute
`\log(1 - a_i p^{(v+1)*2^i})` using the standard Taylor expansion
.. MATH::
\log(1 - x) = -x - 1/2 x^2 - 1/3 x^3 - 1/4 x^4 - 1/5 x^5 - \cdots
together with a binary splitting method.
3. Divide the result by `p^v`
The complexity of this algorithm is quasi-linear.
EXAMPLES::
sage: r = Qp(5,prec=4)(6)
sage: r._log_binary_splitting(2)
5 + O(5^2)
sage: r._log_binary_splitting(4)
5 + 2*5^2 + 4*5^3 + O(5^4)
sage: r._log_binary_splitting(100)
5 + 2*5^2 + 4*5^3 + O(5^4)
sage: r = Zp(5,prec=4,type='fixed-mod')(6)
sage: r._log_binary_splitting(5)
5 + 2*5^2 + 4*5^3 + O(5^4)
"""
cdef unsigned long p
cdef unsigned long prec = aprec
cdef pAdicFloatingPointElement ans
if mpz_fits_slong_p(self.prime_pow.prime.value) == 0:
raise NotImplementedError("The prime %s does not fit in a long" % self.prime_pow.prime)
p = self.prime_pow.prime
ans = self._new_c()
ans.ordp = max((self-1).valuation(), (self+1).valuation())
sig_on()
padiclog(ans.unit, self.unit, p, prec, self.prime_pow.pow_mpz_t_tmp(prec))
mpz_divexact(ans.unit, ans.unit, self.prime_pow.pow_mpz_t_tmp(ans.ordp))
sig_off()
return ans
......@@ -1523,7 +1523,7 @@ cdef class pAdicGenericElement(LocalGenericElement):
r = rational_reconstruction(alpha, m)
return (Rational(p)**self.valuation())*r
def _log_generic(self, aprec, mina=0):
def _log_generic(self, aprec):
r"""
Return ``\log(self)`` for ``self`` equal to 1 in the residue field
......@@ -1535,10 +1535,6 @@ cdef class pAdicGenericElement(LocalGenericElement):
correct. ``aprec`` must not exceed the precision cap of the ring over
which this element is defined.
- ``mina`` -- an integer (default: 0), the series will check `n` up to
this valuation (and beyond) to see if they can contribute to the
series.
ALGORITHM:
The computation uses the series
......@@ -1588,6 +1584,28 @@ cdef class pAdicGenericElement(LocalGenericElement):
e=R.ramification_index()
p=R.prime()
mina = 0
if e != 1:
xval = x.valuation()
lamb = aprec - xval
if lamb > 0 and lamb*(p-1) <= e:
# This is the precision region where the absolute
# precision of the answer might be less than the
# absolute precision of the input
# kink is the number of times we multiply the relative
# precision by p before starting to add e instead.
kink = (e // (lamb * (p-1))).exact_log(p) + 1
# deriv0 is within 1 of the n yielding the minimal
# absolute precision
deriv0 = (e / (aprec * p.log(prec=53))).floor().exact_log(p)
# These are the absolute precisions of x^(p^n) at potential minimum points
L = [(aprec * p**n - n * e, n) for n in [0, kink, deriv0, deriv0+1]]
L.sort()
mina = L[0][1]
# we sum all terms of the power series of log into total
total=R.zero()
......@@ -1656,7 +1674,7 @@ cdef class pAdicGenericElement(LocalGenericElement):
return total.add_bigoh(aprec)
def _log_binary_splitting(self, aprec, mina=0):
def _log_binary_splitting(self, aprec):
r"""
Return ``\log(self)`` for ``self`` equal to 1 in the residue field
......@@ -1719,7 +1737,7 @@ cdef class pAdicGenericElement(LocalGenericElement):
5 + 2*5^2 + 4*5^3 + O(5^4)
"""
raise NotImplementedError
raise NotImplementedError("The binary splitting algorithm is not implemented over the ring: %s" % self.parent())
def log(self, p_branch=None, pi_branch=None, aprec=None, change_frac=False, algorithm=None):
r"""
......@@ -2097,56 +2115,46 @@ cdef class pAdicGenericElement(LocalGenericElement):
pi_branch = (p_branch - R._log_unit_part_p()) / R.e()
total = self.valuation() * pi_branch
y = self.unit_part()
x = 1 - y
if x.valuation()>0:
if (y-1).valuation() > 0:
denom=Integer(1)
else:
y=y**(q-1) # Is it better to multiply it by Teichmuller element?
denom=Integer(q-1)
x = 1 - y
minaprec = y.precision_absolute()
minn = 0
e = R.e()
if e != 1:
xval = x.valuation()
lamb = minaprec - xval
if lamb > 0 and lamb*(p-1) <= e:
# This is the precision region where the absolute
# precision of the answer might be less than the
# absolute precision of the input
# kink is the number of times we multiply the relative
# precision by p before starting to add e instead.
kink = (e // (lamb * (p-1))).exact_log(p) + 1
# deriv0 is within 1 of the n yielding the minimal
# absolute precision
deriv0 = (e / (minaprec * p.log(prec=53))).floor().exact_log(p)
# These are the absolute precisions of x^(p^n) at potential minimum points
L = [(minaprec * p**n - n * e, n) for n in [0, kink, deriv0, deriv0+1]]
L.sort()
minaprec = L[0][0]
minn = L[0][1]
if aprec is None or aprec > minaprec:
aprec=minaprec
if algorithm is None:
try:
# The binary splitting algorithm is supposed to be faster
log_unit = y._log_binary_splitting(aprec, minn)
except NotImplementedError:
log_unit = y._log_generic(aprec, minn)
def _log(y,aprec):
try:
# The binary splitting algorithm is supposed to be faster
return y._log_binary_splitting(aprec)
except NotImplementedError:
return y._log_generic(aprec)
elif algorithm == "generic":
log_unit = y._log_generic(aprec, minn)
_log = self.__class__._log_generic
elif algorithm == "binary_splitting":
log_unit = y._log_binary_splitting(aprec, minn)
_log = self.__class__._log_binary_splitting
else:
raise ValueError("Algorithm must be either 'generic', 'binary_splitting' or None")
if R.is_floating_point() and aprec is None:
val_guessed = (y-1).valuation()
if R.prime() == 2:
val_guessed = max(val_guessed, (y+1).valuation())
if val_guessed is infinity:
log_unit = R(0)
else:
while True: # should we have a maximal precision?
aprec = val_guessed + R.precision_cap()
log_unit = _log(y, aprec)
val = min(aprec, log_unit.valuation())
if val == val_guessed: break
val_guessed = val
else:
minaprec = y.precision_absolute()
if aprec is None or aprec > minaprec:
aprec = minaprec
log_unit = _log(y, aprec)
retval = total + log_unit*R(denom).inverse_of_unit()
if not change_frac:
if retval.valuation() < 0 and not R.is_field():
......