summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--src/sage/categories/pushout.py2
-rw-r--r--src/sage/rings/padics/all.py4
-rw-r--r--src/sage/rings/padics/factory.py523
-rw-r--r--src/sage/rings/padics/generic_nodes.py10
-rw-r--r--src/sage/rings/padics/lattice_precision.py1362
-rw-r--r--src/sage/rings/padics/padic_base_generic.py11
-rw-r--r--src/sage/rings/padics/padic_base_leaves.py693
-rw-r--r--src/sage/rings/padics/padic_lattice_element.py896
-rw-r--r--src/sage/rings/padics/padic_printing.pyx7
9 files changed, 3461 insertions, 47 deletions
diff --git a/src/sage/categories/pushout.py b/src/sage/categories/pushout.py
index 2d845b9..ecaf8f0 100644
--- a/src/sage/categories/pushout.py
+++ b/src/sage/categories/pushout.py
@@ -2467,7 +2467,7 @@ class CompletionFunctor(ConstructionFunctor):
return not (self == other)
_real_types = ['Interval','Ball','MPFR','RDF','RLF']
- _dvr_types = [None, 'fixed-mod','floating-point','capped-abs','capped-rel','lazy']
+ _dvr_types = [None, 'fixed-mod','floating-point','capped-abs','capped-rel','lazy','lattice-cap']
def merge(self, other):
"""
diff --git a/src/sage/rings/padics/all.py b/src/sage/rings/padics/all.py
index 809b4aa..cf9be74 100644
--- a/src/sage/rings/padics/all.py
+++ b/src/sage/rings/padics/all.py
@@ -1,7 +1,7 @@
from __future__ import absolute_import
from .generic_nodes import is_pAdicField, is_pAdicRing
-from .factory import Zp, Zq, Zp as pAdicRing, ZpCR, ZpCA, ZpFM, ZpFP, ZqCR, ZqCA, ZqFM, ZqFP #, ZpL, ZqL
-from .factory import Qp, Qq, Qp as pAdicField, QpCR, QpFP, QqCR, QqFP #, QpL, QqL
+from .factory import Zp, Zq, Zp as pAdicRing, ZpCR, ZpCA, ZpFM, ZpFP, ZpLC, ZqCR, ZqCA, ZqFM, ZqFP #, ZpL, ZqL
+from .factory import Qp, Qq, Qp as pAdicField, QpCR, QpFP, QpLC, QqCR, QqFP #, QpL, QqL
from .factory import pAdicExtension
from .padic_generic import local_print_mode
from .pow_computer import PowComputer
diff --git a/src/sage/rings/padics/factory.py b/src/sage/rings/padics/factory.py
index 291172c..792bc48 100644
--- a/src/sage/rings/padics/factory.py
+++ b/src/sage/rings/padics/factory.py
@@ -22,6 +22,7 @@ from __future__ import absolute_import, print_function
from sage.structure.factory import UniqueFactory
from sage.rings.integer import Integer
+from sage.rings.infinity import Infinity
from sage.structure.factorization import Factorization
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
@@ -31,8 +32,10 @@ from .padic_base_leaves import (pAdicRingCappedRelative,
pAdicRingCappedAbsolute,
pAdicRingFixedMod,
pAdicRingFloatingPoint,
+ pAdicRingLattice,
pAdicFieldCappedRelative,
- pAdicFieldFloatingPoint)
+ pAdicFieldFloatingPoint,
+ pAdicFieldLattice)
from . import padic_printing
######################################################
@@ -75,7 +78,7 @@ def _default_show_prec(type, print_mode):
INPUT:
- - ``type`` -- a string: ``'capped-rel'``, ``'capped-abs'``, ``'fixed-mod'`` or ``'floating-point'``
+ - ``type`` -- a string: ``'capped-rel'``, ``'capped-abs'``, ``'fixed-mod'``, ``'floating-point'`` or ``'lattice-cap'``
- ``print_mode`` -- a string: ``'series'``, ``'terse'``, ``'val-unit'``, ``'digits'``, ``'bars'``
EXAMPLES::
@@ -95,7 +98,7 @@ def _default_show_prec(type, print_mode):
else:
return False
-def get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, check, valid_non_lazy_types):
+def get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, check, valid_non_lazy_types, label=None, absprec=None, relprec=None):
"""
This implements create_key for Zp and Qp: moving it here prevents code duplication.
@@ -105,19 +108,72 @@ def get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_se
sage: from sage.rings.padics.factory import get_key_base
sage: get_key_base(11, 5, 'capped-rel', None, None, None, None, ':', None, None, False, True, ['capped-rel'])
- (11, 5, 'capped-rel', 'series', '11', True, '|', (), -1, False)
+ (11, 5, 'capped-rel', 'series', '11', True, '|', (), -1, False, None)
sage: get_key_base(12, 5, 'capped-rel', 'digits', None, None, None, None, None, None, True, False, ['capped-rel'])
- (12, 5, 'capped-rel', 'digits', '12', True, '|', ('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B'), -1, True)
+ (12,
+ 5,
+ 'capped-rel',
+ 'digits',
+ '12',
+ True,
+ '|',
+ ('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B'),
+ -1,
+ True,
+ None)
"""
- if prec is None:
- prec = DEFAULT_PREC
if check:
if not isinstance(p, Integer):
p = Integer(p)
- if not isinstance(prec, Integer):
- prec = Integer(prec)
if not p.is_prime():
raise ValueError("p must be prime")
+ if type == 'lattice-cap':
+ relative_cap = absolute_cap = None
+ if prec is not None:
+ # We first try to unpack
+ try:
+ relative_cap, absolute_cap = prec
+ except (ValueError, TypeError):
+ relative_cap = prec
+ if absprec is not None:
+ if absolute_cap is None:
+ absolute_cap = absprec
+ else:
+ raise ValueError("absolute cap specified twice")
+ if relprec is not None:
+ if relative_cap is None:
+ relative_cap = relprec
+ else:
+ raise ValueError("relative cap specified twice")
+ if relative_cap is not None:
+ if relative_cap is not Infinity:
+ try:
+ relative_cap = Integer(relative_cap)
+ except TypeError:
+ raise TypeError("relative cap must be either a positive integer or infinity")
+ if relative_cap <= 0:
+ raise ValueError("relative cap must be positive")
+ if absolute_cap is not None:
+ try:
+ absolute_cap = Integer(absolute_cap)
+ except TypeError:
+ raise TypeError("absolute cap must be an integer")
+ if relative_cap is None and absolute_cap is None:
+ relative_cap = DEFAULT_PREC
+ absolute_cap = 2 * DEFAULT_PREC
+ elif relative_cap is None:
+ relative_cap = Infinity
+ elif absolute_cap is None:
+ absolute_cap = 2 * relative_cap
+ prec = (relative_cap, absolute_cap)
+ else:
+ if prec is not None:
+ prec = Integer(prec)
+ if prec is None:
+ if type == 'lattice-cap':
+ prec = (DEFAULT_PREC, 2*DEFAULT_PREC)
+ else:
+ prec = DEFAULT_PREC
print_ram_name = ram_name
if isinstance(print_mode, dict):
if 'pos' in print_mode:
@@ -186,9 +242,8 @@ def get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_se
if show_prec is None:
show_prec = _default_show_prec(type, print_mode)
if type in valid_non_lazy_types:
- key = (p, prec, type, print_mode, name, print_pos, print_sep, tuple(print_alphabet), print_max_terms, show_prec)
+ key = (p, prec, type, print_mode, name, print_pos, print_sep, tuple(print_alphabet), print_max_terms, show_prec, label)
else:
- print(type)
raise ValueError("type must be %s"%(", ".join(valid_non_lazy_types)))
return key
@@ -537,7 +592,8 @@ class Qp_class(UniqueFactory):
"""
def create_key(self, p, prec = None, type = 'capped-rel', print_mode = None,
names = None, ram_name = None, print_pos = None,
- print_sep = None, print_alphabet = None, print_max_terms = None, show_prec=None, check = True):
+ print_sep = None, print_alphabet = None, print_max_terms = None, show_prec = None, check = True,
+ label = None, relprec = None, absprec = None): # specific to Lattice precision
"""
Creates a key from input parameters for ``Qp``.
@@ -546,7 +602,7 @@ class Qp_class(UniqueFactory):
TESTS::
sage: Qp.create_key(5,40)
- (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1, True)
+ (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1, True, None)
"""
if isinstance(names, (int, Integer)):
# old pickle; names is what used to be halt.
@@ -556,7 +612,7 @@ class Qp_class(UniqueFactory):
print_alphabet = print_max_terms
print_max_terms = check
check = True
- return get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, check, ['capped-rel', 'floating-point'])
+ return get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, check, ['capped-rel', 'floating-point', 'lattice-cap'], label, relprec, absprec)
def create_object(self, version, key):
"""
@@ -575,8 +631,9 @@ class Qp_class(UniqueFactory):
elif version[0] < 8:
p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms = key
show_prec = None
+ label = None
else:
- p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, label = key
if isinstance(type, Integer):
# lazy
raise NotImplementedError("lazy p-adics need more work. Sorry.")
@@ -592,7 +649,7 @@ class Qp_class(UniqueFactory):
return obj
except KeyError:
pass
- p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, label = key
if type == 'capped-rel':
if print_mode == 'terse':
@@ -608,6 +665,13 @@ class Qp_class(UniqueFactory):
else:
return pAdicFieldFloatingPoint(p, prec, {'mode': print_mode, 'pos': print_pos, 'sep': print_sep, 'alphabet': print_alphabet,
'ram_name': name, 'max_ram_terms': print_max_terms, 'show_prec': show_prec}, name)
+ elif type == 'lattice-cap':
+ if print_mode == 'terse':
+ return pAdicFieldLattice(p, prec, {'mode': print_mode, 'pos': print_pos, 'sep': print_sep, 'alphabet': print_alphabet,
+ 'ram_name': name, 'max_terse_terms': print_max_terms, 'show_prec': show_prec}, name, label)
+ else:
+ return pAdicFieldLattice(p, prec, {'mode': print_mode, 'pos': print_pos, 'sep': print_sep, 'alphabet': print_alphabet,
+ 'ram_name': name, 'max_ram_terms': print_max_terms, 'show_prec': show_prec}, name, label)
else:
raise ValueError("unexpected type")
@@ -1181,6 +1245,30 @@ def QpFP(p, prec = None, *args, **kwds):
"""
return Qp(p, prec, 'floating-point', *args, **kwds)
+#def QpL(p, prec = DEFAULT_PREC, print_mode = None, halt = DEFAULT_HALT, names = None, print_pos = None,
+# print_sep = None, print_alphabet = None, print_max_terms = None, check=True):
+# """
+# A shortcut function to create lazy p-adic fields.
+
+# Currently deactivated. See documentation for Qp for a description of the input parameters.
+
+# EXAMPLES::
+
+def QpLC(p, prec = None, *args, **kwds):
+ """
+ A shortcut function to create `p`-adic fields with lattice precision.
+
+ See :func:`ZpLC` for more information about this model of precision.
+
+ EXAMPLES::
+
+ sage: R = QpLC(2)
+ sage: R
+ 2-adic Field with lattice precision
+
+ """
+ return Qp(p, prec, 'lattice-cap', *args, **kwds)
+
def QqCR(q, prec = None, *args, **kwds):
"""
A shortcut function to create capped relative unramified `p`-adic
@@ -1211,6 +1299,22 @@ def QqFP(q, prec = None, *args, **kwds):
"""
return Qq(q, prec, 'floating-point', *args, **kwds)
+def QpLC(p, prec = None, *args, **kwds):
+ """
+ A shortcut function to create `p`-adic fields with lattice precision.
+
+ See :func:`ZpLC` for more information about this model of precision.
+
+ EXAMPLES::
+
+ sage: R = QpLC(2)
+ sage: R
+ 2-adic Field with lattice precision
+
+ """
+ return Qp(p, prec, 'lattice-cap', *args, **kwds)
+
+
#######################################################################################################
#
# p-Adic Rings
@@ -1229,14 +1333,16 @@ class Zp_class(UniqueFactory):
- ``p`` -- integer: the `p` in `\mathbb{Z}_p`
- ``prec`` -- integer (default: ``20``) the precision cap of the
- ring. Except for the fixed modulus case, individual elements
+ ring. In the lattice capped case, ``prec`` can either be a
+ pair (``relative_cap``, ``absolute_cap``).
+ Except for the fixed modulus case, individual elements
keep track of their own precision. See TYPES and PRECISION
below.
- ``type`` -- string (default: ``'capped-rel'``) Valid types are
``'capped-rel'``, ``'capped-abs'``, ``'fixed-mod'``,
- ``'floating-point'`` and ``'lazy'`` (though lazy is not yet
- implemented). See TYPES and PRECISION below
+ ``'floating-point'`` and ``'lattice-cap'``.
+ See TYPES and PRECISION below
- ``print_mode`` -- string (default: ``None``). Valid modes are
``'series'``, ``'val-unit'``, ``'terse'``, ``'digits'``, and
@@ -1271,9 +1377,9 @@ class Zp_class(UniqueFactory):
TYPES AND PRECISION:
- There are two types of precision for a `p`-adic element. The first
- is relative precision, which gives the number of known `p`-adic
- digits::
+ There are three types of precision.
+ The first is relative precision; it gives the number of known
+ `p`-adic digits::
sage: R = Zp(5, 20, 'capped-rel', 'series'); a = R(675); a
2*5^2 + 5^4 + O(5^22)
@@ -1286,11 +1392,29 @@ class Zp_class(UniqueFactory):
sage: a.precision_absolute()
22
- There are five types of `p`-adic rings: capped relative rings
+ The third one is lattice precision.
+ It is not attached to a single `p`-adic number but is a unique
+ object modeling the precision on a set of `p`-adics, which is
+ typically the set of all elements within the same parent.
+
+ sage: R = ZpLC(17)
+ sage: x = R(1,10); y = R(1,5)
+ sage: R.precision()
+ Precision Lattice on 2 objects
+ sage: R.precision().precision_lattice()
+ [2015993900449 0]
+ [ 0 1419857]
+
+ We refer to the documentation of the function :func:`ZpLC` for
+ more information about this precision model.
+
+ ::
+
+ There are six types of `p`-adic rings: capped relative rings
(type= ``'capped-rel'``), capped absolute rings
(type= ``'capped-abs'``), fixed modulus rings (type= ``'fixed-mod'``),
- floating point rings (type=``'floating-point'``),
- and lazy rings (type= ``'lazy'``).
+ floating point rings (type=``'floating-point'``), lattice capped
+ rings (type=``'lattice-cap'``) and lazy rings (type= ``'lazy'``).
In the capped relative case, the relative precision of an element
is restricted to be at most a certain value, specified at the
@@ -1337,10 +1461,16 @@ class Zp_class(UniqueFactory):
in that elements do not trac their own precision. However, relative
precision is truncated with each operation rather than absolute precision.
+ On the contrary, the lattice type tracks precision using lattices
+ and automatic differentiation. It is rather slow but provides sharp
+ (often optimal) results regarding precision.
+ We refer to the documentation of the function :func:`ZpLC` for a
+ small demonstration of the capabilities of this precision model.
+
The lazy case will eventually support elements that can increase
their precision upon request. It is not currently implemented.
- PRINTING
+ PRINTING:
There are many different ways to print `p`-adic elements. The
way elements of a given ring print is controlled by options
@@ -1603,7 +1733,8 @@ class Zp_class(UniqueFactory):
"""
def create_key(self, p, prec = None, type = 'capped-rel', print_mode = None,
names = None, ram_name = None, print_pos = None, print_sep = None, print_alphabet = None,
- print_max_terms = None, show_prec = None, check = True):
+ print_max_terms = None, show_prec = None, check = True,
+ label = None, relprec = None, absprec = None):
"""
Creates a key from input parameters for ``Zp``.
@@ -1612,9 +1743,19 @@ class Zp_class(UniqueFactory):
TESTS::
sage: Zp.create_key(5,40)
- (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1, True)
+ (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1, True, None)
sage: Zp.create_key(5,40,print_mode='digits')
- (5, 40, 'capped-rel', 'digits', '5', True, '|', ('0', '1', '2', '3', '4'), -1, False)
+ (5,
+ 40,
+ 'capped-rel',
+ 'digits',
+ '5',
+ True,
+ '|',
+ ('0', '1', '2', '3', '4'),
+ -1,
+ False,
+ None)
"""
if isinstance(names, (int, Integer)):
# old pickle; names is what used to be halt.
@@ -1625,7 +1766,8 @@ class Zp_class(UniqueFactory):
print_max_terms = check
check = True
return get_key_base(p, prec, type, print_mode, names, ram_name, print_pos, print_sep, print_alphabet,
- print_max_terms, show_prec, check, ['capped-rel', 'fixed-mod', 'capped-abs', 'floating-point'])
+ print_max_terms, show_prec, check, ['capped-rel', 'fixed-mod', 'capped-abs', 'floating-point', 'lattice-cap'],
+ label=label, relprec=relprec, absprec=absprec)
def create_object(self, version, key):
"""
@@ -1645,8 +1787,9 @@ class Zp_class(UniqueFactory):
elif version[0] < 8:
p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms = key
show_prec = None
+ label = None
else:
- p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, label = key
if isinstance(type, Integer):
# lazy
raise NotImplementedError("lazy p-adics need more work. Sorry.")
@@ -1655,14 +1798,14 @@ class Zp_class(UniqueFactory):
# keys changed in order to reduce irrelevant duplications: e.g. two Zps with print_mode 'series'
# that are identical except for different 'print_alphabet' now return the same object.
key = get_key_base(p, prec, type, print_mode, name, None, print_pos, print_sep, print_alphabet,
- print_max_terms, None, False, ['capped-rel', 'fixed-mod', 'capped-abs'])
+ print_max_terms, None, False, ['capped-rel', 'fixed-mod', 'capped-abs', 'lattice-cap'])
try:
obj = self._cache[version, key]()
if obj is not None:
return obj
except KeyError:
pass
- p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, show_prec, label = key
if type == 'capped-rel':
return pAdicRingCappedRelative(p, prec, {'mode': print_mode, 'pos': print_pos, 'sep': print_sep, 'alphabet': print_alphabet,
'ram_name': name, 'max_ram_terms': print_max_terms, 'show_prec': show_prec}, name)
@@ -1675,6 +1818,9 @@ class Zp_class(UniqueFactory):
elif type == 'floating-point':
return pAdicRingFloatingPoint(p, prec, {'mode': print_mode, 'pos': print_pos, 'sep': print_sep, 'alphabet': print_alphabet,
'ram_name': name, 'max_ram_terms': print_max_terms, 'show_prec': show_prec}, name)
+ elif type == 'lattice-cap':
+ return pAdicRingLattice(p, prec, {'mode': print_mode, 'pos': print_pos, 'sep': print_sep, 'alphabet': print_alphabet,
+ 'ram_name': name, 'max_ram_terms': print_max_terms, 'show_prec': show_prec}, name, label)
else:
raise ValueError("unexpected type")
@@ -2343,6 +2489,317 @@ def ZqFP(q, prec = None, *args, **kwds):
"""
return Zq(q, prec, 'floating-point', *args, **kwds)
+def ZpLC(p, prec=None, *args, **kwds):
+ """
+ A shortcut function to create `p`-adic rings with lattice precision
+ with cap.
+
+ DEMONSTRATION:
+
+ Below is a small demo of the features by this model of precision::
+
+ sage: R = ZpLC(3, print_mode='terse')
+ sage: x = R(1,10)
+
+ Of course, when we multiply by 3, we gain one digit of absolute
+ precision::
+
+ sage: 3*x
+ 3 + O(3^11)
+
+ The lattice precision machinery sees this even if we decompose
+ the computation into several steps::
+
+ sage: y = x+x
+ sage: y
+ 2 + O(3^10)
+ sage: x + y
+ 3 + O(3^11)
+
+ The same works for the multiplication::
+
+ sage: z = x^2
+ sage: z
+ 1 + O(3^10)
+ sage: x*z
+ 1 + O(3^11)
+
+ This comes more funny when we are working with elements given
+ at different precisions::
+
+ sage: R = ZpLC(2, print_mode='terse')
+ sage: x = R(1,10)
+ sage: y = R(1,5)
+ sage: z = x+y; z
+ 2 + O(2^5)
+ sage: t = x-y; t
+ 0 + O(2^5)
+ sage: z+t # observe that z+t = 2*x
+ 2 + O(2^11)
+ sage: z-t # observe that z-t = 2*y
+ 2 + O(2^6)
+
+ sage: x = R(28888,15)
+ sage: y = R(204,10)
+ sage: z = x/y; z
+ 242 + O(2^9)
+ sage: z*y # which is x
+ 28888 + O(2^15)
+
+ The SOMOS sequence is the sequence defined by the recurrence::
+
+ ..MATH::
+
+ u_n = \frac {u_{n-1} u_{n-3} + u_{n-2}^2} {u_{n-4}}
+
+ It is known for its numerical instability.
+ On the one hand, one can show that if the initial values are
+ invertible in `\mathbb{Z}_p` and known at precision `O(p^N)`
+ then all the next terms of the SOMOS sequence will be known
+ at the same precision as well.
+ On the other hand, because of the division, when we unroll
+ the recurrence, we loose a lot of precision. Observe::
+
+ sage: R = Zp(2, 30, print_mode='terse')
+ sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 4 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 13 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 55 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 21975 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 6639 + O(2^13)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 7186 + O(2^13)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 569 + O(2^13)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 253 + O(2^13)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 4149 + O(2^13)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 2899 + O(2^12)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 3072 + O(2^12)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 349 + O(2^12)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 619 + O(2^12)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 243 + O(2^12)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 3 + O(2^2)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 2 + O(2^2)
+
+ If instead, we use the lattice precision, everything goes well::
+
+ sage: R = ZpLC(2, 30, print_mode='terse')
+ sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 4 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 13 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 55 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 21975 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 23023 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 31762 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 16953 + O(2^15)
+ sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print(d)
+ 16637 + O(2^15)
+
+ sage: for _ in range(100):
+ ....: a,b,c,d = b,c,d,(b*d+c*c)/a
+ sage: a
+ 15519 + O(2^15)
+ sage: b
+ 32042 + O(2^15)
+ sage: c
+ 17769 + O(2^15)
+ sage: d
+ 20949 + O(2^15)
+
+ BEHIND THE SCENE:
+
+ The precision is global.
+ It is encoded by a lattice in a huge vector space whose dimension
+ is the number of elements having this parent.
+
+ Concretely, this precision datum is an instance of the class
+ :class:`sage.rings.padic.lattice_precision.PrecisionLattice`.
+ It is attached to the parent and is created at the same time
+ as the parent.
+ (It is actually a bit more subtle because two different parents
+ may share the same instance; this happens for instance for a
+ `p`-adic ring and its field of fractions.)
+
+ This precision datum is accessible through the method :meth:`precision`::
+
+ sage: R = ZpLC(5, print_mode='terse')
+ sage: prec = R.precision()
+ sage: prec
+ Precision Lattice on 0 object
+
+ This instance knows about all elements of the parent, it is
+ automatically updated when a new element (of this parent) is
+ created::
+
+ sage: x = R(3513,10)
+ sage: prec
+ Precision Lattice on 1 object
+ sage: y = R(176,5)
+ sage: prec
+ Precision Lattice on 2 objects
+ sage: z = R.random_element()
+ sage: prec
+ Precision Lattice on 3 objects
+
+ The method :meth:`tracked_elements` provides the list of all
+ tracked elements::
+
+ sage: prec.tracked_elements()
+ [3513 + O(5^10), 176 + O(5^5), ...]
+
+ Similarly, when a variable is collected by the garbage collector,
+ the precision lattice is updated. Note however that the update
+ might be delayed. We can force it with the method :meth:`del_elements`::
+
+ sage: z = 0
+ sage: prec
+ Precision Lattice on 3 objects
+ sage: prec.del_elements()
+ sage: prec
+ Precision Lattice on 2 objects
+
+ The method :meth:`precision_lattice` returns (a matrix defining)
+ the lattice that models the precision. Here we have::
+
+ sage: prec.precision_lattice()
+ [9765625 0]
+ [ 0 3125]
+
+ Observe that `5^10 = 9765625` and `5^5 = 3125`.
+ The above matrix then reflects the precision on `x` and `y`.
+
+ Now, observe how the precision lattice changes while performing
+ computations::
+
+ sage: x, y = 3*x+2*y, 2*(x-y)
+ sage: prec.del_elements()
+ sage: prec.precision_lattice()
+ [ 3125 48825000]
+ [ 0 48828125]
+
+ The matrix we get is no longer diagonal, meaning that some digits
+ of precision are diffused among the two new elements `x` and `y`.
+ They nevertheless show up when we compute for instance `x+y`::
+
+ sage: x
+ 1516 + O(5^5)
+ sage: y
+ 424 + O(5^5)
+ sage: x+y
+ 17565 + O(5^11)
+
+ These diffused digits of precision (which are tracked but
+ do not appear on the printing) allow to be always sharp on
+ precision.
+
+ PERFORMANCES:
+
+ Each elementary operation requires significant manipulations
+ on the lattice precision and then is costly. Precisely:
+
+ - The creation of a new element has a cost `O(n)` where `n`
+ is the number of tracked elements.
+
+ - The destruction of one element has a cost `O(m^2)` where
+ `m` is the distance between the destroyed element and
+ the last one. Fortunately, it seems that `m` tends to
+ be small in general (the dynamics of the list of tracked
+ elements is rather close to that of a stack).
+
+ It is nevertheless still possible to manipulate several
+ hundred variables (e.g. square matrices of size 5 or
+ polynomials of degree 20 are accessible).
+
+ The class :class:`PrecisionLattice` provides several
+ features for introspection (especially concerning timings).
+ If enabled, it maintains an history of all actions and stores
+ the wall time of each of them::
+
+ sage: R = ZpLC(3)
+ sage: prec = R.precision()
+ sage: prec.history_enable()
+ sage: M = random_matrix(R, 5)
+ sage: d = M.determinant()
+ sage: print(prec.history()) # somewhat random
+ ---
+ 0.004212s oooooooooooooooooooooooooooooooooooo
+ 0.000003s oooooooooooooooooooooooooooooooooo~~
+ 0.000010s oooooooooooooooooooooooooooooooooo
+ 0.001560s ooooooooooooooooooooooooooooooooooooooooo
+ 0.000004s ooooooooooooooooooooooooooooo~oooo~oooo~o
+ 0.002168s oooooooooooooooooooooooooooooooooooooo
+ 0.001787s ooooooooooooooooooooooooooooooooooooooooo
+ 0.000004s oooooooooooooooooooooooooooooooooooooo~~o
+ 0.000198s ooooooooooooooooooooooooooooooooooooooo
+ 0.001152s ooooooooooooooooooooooooooooooooooooooooo
+ 0.000005s ooooooooooooooooooooooooooooooooo~oooo~~o
+ 0.000853s oooooooooooooooooooooooooooooooooooooo
+ 0.000610s ooooooooooooooooooooooooooooooooooooooo
+ ...
+ 0.003879s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.000006s oooooooooooooooooooooooooooooooooooooooooooooooooooo~~~~~
+ 0.000036s oooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.006737s oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.000005s oooooooooooooooooooooooooooooooooooooooooooooooooooo~~~~~ooooo
+ 0.002637s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.007118s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.000008s oooooooooooooooooooooooooooooooooooooooooooooooooooo~~~~o~~~~oooo
+ 0.003504s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.005371s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.000006s ooooooooooooooooooooooooooooooooooooooooooooooooooooo~~~o~~~ooo
+ 0.001858s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.003584s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.000004s oooooooooooooooooooooooooooooooooooooooooooooooooooooo~~o~~oo
+ 0.000801s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.001916s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
+ 0.000022s ooooooooooooooooooooooooooooo~~~~~~~~~~~~~~~~~~~~~~oooo~o~o
+ 0.014705s ooooooooooooooooooooooooooooooooooo
+ 0.001292s ooooooooooooooooooooooooooooooooooooo
+ 0.000002s ooooooooooooooooooooooooooooooooooo~o
+
+ The symbol `o` symbolized a tracked element.
+ The symbol `~` means that the element is marked for deletion.
+
+ The global timings are also accessible as follows::
+
+ sage: prec.timings() # somewhat random
+ {'add': 0.25049376487731934,
+ 'del': 0.11911273002624512,
+ 'mark': 0.0004909038543701172,
+ 'partial reduce': 0.0917658805847168}
+
+ """
+ return Zp(p, prec, 'lattice-cap', *args, **kwds)
+
+#def ZpL(p, prec = DEFAULT_PREC, print_mode = None, halt = DEFAULT_HALT, names = None, print_pos = None,
+# print_sep = None, print_alphabet = None, print_max_terms = None, check=True):
+# """
+# A shortcut function to create lazy `p`-adic rings.
+
+# Currently deactivated. See documentation for Zp for a description of the input parameters.
+
+
#######################################################################################################
#
# The Extension Factory -- creates extensions of p-adic rings and fields
diff --git a/src/sage/rings/padics/generic_nodes.py b/src/sage/rings/padics/generic_nodes.py
index 656faaa..ba6ad41 100644
--- a/src/sage/rings/padics/generic_nodes.py
+++ b/src/sage/rings/padics/generic_nodes.py
@@ -389,10 +389,12 @@ class pAdicRingBaseGeneric(pAdicBaseGeneric, pAdicRingGeneric):
True
"""
from sage.categories.pushout import CompletionFunctor
- return (CompletionFunctor(self.prime(),
- self.precision_cap(),
- {'print_mode':self._printer.dict(), 'type':self._prec_type(), 'names':self._names}),
- ZZ)
+ extras = {'print_mode':self._printer.dict(), 'type':self._prec_type(), 'names':self._names}
+ try:
+ extras['label'] = self._label
+ except AttributeError:
+ pass
+ return (CompletionFunctor(self.prime(), self.precision_cap(), extras), ZZ)
def random_element(self, algorithm='default'):
r"""
diff --git a/src/sage/rings/padics/lattice_precision.py b/src/sage/rings/padics/lattice_precision.py
new file mode 100644
index 00000000..afba9a2
--- /dev/null
+++ b/src/sage/rings/padics/lattice_precision.py
@@ -0,0 +1,1362 @@
+import _weakref as weakref
+from sage.misc.misc import walltime
+
+from sage.structure.sage_object import SageObject
+from sage.structure.unique_representation import UniqueRepresentation
+
+from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+from sage.rings.infinity import Infinity
+
+
+# Class pRational
+#################
+
+# It is a temporary class implementing rational numbers
+# as approximations of p-adic numbers
+
+class pRational:
+ def __init__(self, p, x, exponent=0, valuation=None):
+ self.p = p
+ if x in ZZ:
+ self.x = ZZ(x)
+ else:
+ self.x = x
+ self.exponent = exponent
+ self._valuation = valuation
+
+ def __repr__(self):
+ if self.exponent == 0:
+ return str(self.x)
+ else:
+ return "%s^%s * %s" % (self.p, self.exponent, self.x)
+
+ def reduce(self, prec):
+ if prec is Infinity:
+ return self
+ x = self.x
+ exp = self.exponent
+ if x.parent() is ZZ:
+ if prec > exp:
+ x = x % (self.p ** (prec-exp))
+ else:
+ x = 0
+ elif x.parent() is QQ:
+ num = x.numerator()
+ denom = x.denominator()
+ valdenom = denom.valuation(self.p)
+ denom /= self.p ** valdenom
+ exp -= valdenom
+ modulo = self.p ** (prec - exp)
+ # probably we should use Newton iteration instead
+ # (but it is actually slower for now - Python implementation)
+ _, inv, _ = denom.xgcd(modulo)
+ x = (num*inv) % modulo
+ if self.x == 0:
+ val = Infinity
+ else:
+ val = self._valuation
+ return self.__class__(self.p, x, exp, valuation=val)
+
+ def normalize(self):
+ if self.x == 0:
+ self.exponent = 0
+ else:
+ val = self.valuation()
+ exp = self.exponent
+ self.x /= self.p ** (val-exp)
+ if self.x in ZZ:
+ self.x = ZZ(self.x)
+ self.exponent = val
+
+ def valuation(self):
+ if self._valuation is None:
+ valx = self.x.valuation(self.p)
+ self._valuation = self.exponent + valx
+ return self._valuation
+
+ def is_p_power(self):
+ self.normalize()
+ return self.x == 1
+
+ def is_zero(self):
+ return self.x == 0
+
+ def __add__(self, other):
+ p = self.p
+ sexp = self.exponent
+ oexp = other.exponent
+ if self._valuation is None or other._valuation is None:
+ val = None
+ elif self._valuation < other._valuation:
+ val = self._valuation
+ elif self._valuation > other._valuation:
+ val = other._valuation
+ else:
+ val = None
+ if sexp < oexp:
+ return self.__class__(p, self.x + other.x * p**(oexp-sexp), sexp, valuation=val)
+ else:
+ return self.__class__(p, self.x * p**(sexp-oexp) + other.x, oexp, valuation=val)
+
+ def __sub__(self, other):
+ p = self.p
+ sexp = self.exponent
+ oexp = other.exponent
+ if self._valuation is None or other._valuation is None:
+ val = None
+ elif self._valuation < other._valuation:
+ val = self._valuation
+ elif self._valuation > other._valuation:
+ val = other._valuation
+ else:
+ val = None
+ if sexp < oexp:
+ return self.__class__(p, self.x - other.x * p**(oexp-sexp), sexp, valuation=val)
+ else:
+ return self.__class__(p, self.x * p**(sexp-oexp) - other.x, oexp, valuation=val)
+
+ def __neg__(self):
+ return self.__class__(self.p, -self.x, self.exponent, valuation=self._valuation)
+
+ def __mul__(self, other):
+ if self._valuation is None or other._valuation is None:
+ val = None
+ else:
+ val = self._valuation + other._valuation
+ return self.__class__(self.p, self.x * other.x, self.exponent + other.exponent, valuation=val)
+
+ def __div__(self, other):
+ if self._valuation is None or other._valuation is None:
+ val = None
+ else:
+ val = self._valuation - other._valuation
+ return self.__class__(self.p, self.x / other.x, self.exponent - other.exponent, valuation=val)
+
+ def __lshift__(self, n):
+ if self._valuation is None:
+ val = None
+ else:
+ val = self._valuation + n
+ return self.__class__(self.p, self.x, self.exponent + n, valuation=val)
+
+ def __rshift__(self, n):
+ return self << (-n)
+
+ def unit_part(self):
+ if self.is_zero():
+ raise ValueError("the unit part of zero is not defined")
+ p = self.p
+ val = self.valuation()
+ x = self.x / (p ** (val-self.exponent))
+ return self.__class__(p, x, 0, valuation=0)
+
+ def xgcd(self,other):
+ p = self.p
+ sexp = self.exponent
+ oexp = other.exponent
+ if sexp < oexp:
+ a = ZZ(self.x)
+ b = ZZ(other.x * (p ** (oexp-sexp)))
+ exp = sexp
+ else:
+ a = ZZ(self.x * (p ** (sexp-oexp)))
+ b = ZZ(other.x)
+ exp = oexp
+ d, u, v = a.xgcd(b)
+ if self._valuation is None or other._valuation is None:
+ val = None
+ else:
+ val = min(self._valuation, other._valuation)
+ d = self.__class__(p, d, exp, valuation=val)
+ u = self.__class__(p, u)
+ v = self.__class__(p, v)
+ return d, u, v
+
+ def value(self):
+ return (self.p ** self.exponent) * self.x
+
+ def list(self, prec):
+ if self.x not in ZZ:
+ raise NotImplementedError
+ val = self.valuation()
+ p = self.p
+ x = ZZ(self.x * p**(self.exponent - val))
+ l = [ ]; i = val
+ while i < prec:
+ x, digit = x.quo_rem(p)
+ l.append(digit)
+ i += 1
+ return l
+
+
+# Class PrecisionLattice
+########################
+
+def list_of_padics(elements):
+ """
+ Convert a list of p-adic composed elements (as polynomials, matrices)
+ to a list of weak refererences of their p-adic coefficients.
+
+ This is an helper function for the methods :meth:`precision_lattice`
+
+ TESTS::
+
+ sage: from sage.rings.padics.lattice_precision import list_of_padics
+ sage: R = ZpLC(2)
+ sage: M = random_matrix(R,2,2)
+ sage: list_of_padics(M)
+ [<weakref at 0x...; to 'pAdicLatticeElement' at 0x...>,
+ <weakref at 0x...; to 'pAdicLatticeElement' at 0x...>,
+ <weakref at 0x...; to 'pAdicLatticeElement' at 0x...>,
+ <weakref at 0x...; to 'pAdicLatticeElement' at 0x...>]
+ """
+ from sage.rings.padics.padic_lattice_element import pAdicLatticeElement
+ if isinstance(elements, pAdicLatticeElement):
+ return [ weakref.ref(elements) ]
+ if not isinstance(elements, list):
+ elements = list(elements)
+ ans = [ ]
+ for x in elements:
+ ans += list_of_padics(x)
+ return ans
+
+def format_history(tme, status, timings):
+ """
+ Return a formated output for the history.
+
+ This is an helper function for the methods :meth:`history`.
+
+ TESTS::
+
+ sage: from sage.rings.padics.lattice_precision import format_history
+ sage: format_history(1.23456789, ['o','o','o','o','o','o','~','o','o'], true)
+ '1.234568s oooooo~oo'
+ sage: format_history(1.23456789, ['o','o','o','o','o','o','~','o','o'], false)
+ 'oooooo~oo'
+
+ sage: format_history(12.3456789, ['o','o','o','o','o','o','~','o','o'], true)
+ ' >= 10s oooooo~oo'
+ sage: format_history(10^(-10), ['o','o','o','o','o','o','~','o','o'], true)
+ ' --- oooooo~oo'
+ sage: format_history(-1, ['o','o','o','o','o','o','~','o','o'], true)
+ ' Timings oooooo~oo'
+ """
+ status = ''.join(status)
+ if timings:
+ if tme < 0:
+ s = " Timings "
+ elif tme < 0.000001:
+ s = " --- "
+ elif tme >= 10:
+ s = " >= 10s "
+ else:
+ s = "%.6fs" % tme
+ return s + " " + status
+ else:
+ return status
+
+class PrecisionLattice(UniqueRepresentation, SageObject):
+ """
+ A class for handling precision lattices which are used to
+ track precision in the ZpLC model.
+
+ The precision lattice is stored as a triangular matrix whose
+ rows are generators of the lattice.
+ """
+ # Internal variables:
+ # . self._cap
+ # a cap for the (working) precision
+ # meaning that the precision lattice always contains p^(self._cap)
+ # . self._elements
+ # list of weak references of elements in this parent
+ # . self._lattice
+ # an upper triangular matrix over ZZ representing the
+ # lattice of precision
+ # (its columns are indexed by self._elements)
+ # . self._absolute_precisions
+ def __init__(self, p, label):
+ self._p = p
+ self._label = label
+ self._elements = [ ]
+ self._capped = { }
+ self._lattice = { }
+ self._absolute_precisions = { }
+ self._marked_for_deletion = [ ]
+ self._approx_zero = pRational(p, ZZ(0))
+ # History
+ self._history_init = None
+ self._history = None
+
+ def label(self):
+ """
+ Return the label of the parent to which this precision lattice
+ corresponds.
+
+ EXAMPLE::
+
+ sage: R = ZpLC(2, label="mylabel")
+ sage: R.precision().label()
+ 'mylabel'
+ """
+ return self._label
+
+ def _repr_(self):
+ """
+ Return a string representation of this precision lattice
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2)
+ sage: R.precision()
+ Precision Lattice on ... objects
+
+ If a label has been specified, it is included in the representation
+
+ sage: R = ZpLC(2, label="mylabel")
+ sage: R.precision()
+ Precision Lattice on 0 object (label: mylabel)
+ """
+ n = len(self._elements)
+ if self._label is None:
+ if n > 1:
+ return "Precision Lattice on %s objects" % len(self._elements)
+ else:
+ return "Precision Lattice on %s object" % len(self._elements)
+ else:
+ if n > 1:
+ return "Precision Lattice on %s objects (label: %s)" % (len(self._elements), self._label)
+ else:
+ return "Precision Lattice on %s object (label: %s)" % (len(self._elements), self._label)
+
+ def prime(self):
+ """
+ Return the underlying prime number attached to this precision lattice.
+
+ EXAMPLE::
+
+ sage: R = ZpLC(2, label="mylabel")
+ sage: R.precision().prime()
+ 2
+ """
+ return self._p
+
+ def reduce(self, index=0, partial=False):
+ """
+ Reduce the size of the entries above the diagonal of the precision matrix
+
+ INPUT:
+
+ - ``index`` -- an integer, the starting row for which the reduction
+ is performed
+
+ - ``partial`` -- a boolean (default: False) specifying whether a
+ partial or a full Hermite reduction should be performed
+
+ NOTE:
+
+ The partial reduction has cost `O(m^2)` where `m` is the number of
+ rows that need to be reduced (that is the difference between the
+ total number of rows and ``index``).
+
+ The full Hermite reduction has cost `O(m^3)`.
+
+ NOTE:
+
+ The software ensures that the precision lattice is always
+ partially reduced.
+ Calling the function manually with the argument ``partial=True``
+ should then just do nothing.
+
+ TESTS::
+
+ sage: R = ZpLC(2)
+ sage: x = R.random_element()
+ sage: del x
+ sage: R.precision().del_elements() # indirect doctest
+ """
+ n = len(self._elements)
+ if index >= n-1:
+ return
+ if partial:
+ # Partial reduction
+ # Cost: O(m^2) with m = n-index
+ tme = walltime()
+ diffval = (n-index) * [0]
+ for j in range(n-1, index, -1):
+ col = self._lattice[self._elements[j]]
+ prec = col[j].valuation() - diffval[j-index]
+ for i in range(index,j):
+ col[i] = col[i].reduce(prec)
+ col[i].normalize() # seems to be faster then
+ dval = col[i].valuation() - prec
+ if dval < diffval[i-index]:
+ diffval[i-index] = dval
+ # We update history
+ if self._history is not None:
+ self._history.append(('partial reduce', index, walltime(tme)))
+ else:
+ # Full Hermite reduction
+ # Cost: O(m^3) with m = n-index
+ tme = walltime()
+ for j in range(index+1, n):
+ # In what follows, we assume that col[j] is a power of p
+ col = self._lattice[self._elements[j]]
+ valpivot = col[j].valuation()
+ for i in range(index, j):
+ reduced = col[i].reduce(valpivot)
+ scalar = (col[i] - reduced) >> valpivot
+ if scalar.is_zero(): continue
+ col[i] = reduced
+ col[i].normalize()
+ for j2 in range(j+1, n):
+ col2 = self._lattice[self._elements[j2]]
+ col2[i] -= scalar*col2[i]
+ col2[i].normalize()
+ # We update history
+ if self._history is not None:
+ self._history.append(('full reduce', index, walltime(tme)))
+
+ def new_element(self, x, dx, bigoh, dx_mode='linear_combinaison', capped=False):
+ """
+ Update the lattice when a new element is created.
+
+ This function is not meant to be called manually.
+ It is automatically called by the parent when a new
+ element is created.
+
+ INPUT:
+
+ - ``x`` -- the newly created element
+
+ - ``dx`` -- a dictionary representing the differential of ``x``
+
+ - ``dx_mode`` -- a string, either ``linear_combinaison`` (the default)
+ or `values`
+
+ - ``capped`` -- a boolean, whether this element has been capped
+ according to the parent's cap
+
+ If ``dx_mode`` is ``linear_combinaison``, the dictionary ``dx``
+ encodes the expression of the differential of ``x``.
+ For example, if ``x`` was defined as ``x = y*z`` then:
+
+ .. MATH::
+
+ dx = y dz + z dy
+
+ and the corresponding dictionary is ``{z: y, y: z}`` (except
+ that the keys are not the elements themselves but weak references
+ to them).
+
+ If ``dx_mode`` is ``values``, the dictionary ``dx`` directly
+ specifies the entries that have to stored in the precision lattice.
+ This mode is only used for multiple conversion between different
+ parents (see :meth:`multiple_conversion`).
+
+ TESTS::
+
+ sage: R = ZpLC(2)
+ sage: x = R.random_element()
+ sage: y = R.random_element()
+ sage: z = x*y # indirect doctest
+ """
+ # First we delete some elements marked for deletion
+ if self._marked_for_deletion:
+ self.del_elements(thresold=50)
+
+ # Then we add the new element
+ tme = walltime()
+ p = self._p
+ n = len(self._elements)
+ x_ref = weakref.ref(x, self.mark_for_deletion)
+ self._elements.append(x_ref)
+ col = n * [self._approx_zero]
+ if dx_mode == 'linear_combinaison':
+ for elt,scalar in dx:
+ ref = weakref.ref(elt)
+ if not isinstance(scalar, pRational):
+ scalar = pRational(p, scalar)
+ c = self._lattice[ref]
+ for i in range(len(c)):
+ col[i] += scalar * c[i]
+ elif dx_mode == 'values':
+ for elt,scalar in dx:
+ ref = weakref.ref(elt)
+ if not isinstance(scalar, pRational):
+ scalar = pRational(p, scalar)
+ i = len(self._lattice[ref]) - 1
+ col[i] = scalar
+ else:
+ raise ValueError("dx_mode must be either 'linear_combinaison' or 'values'")
+ for i in range(n):
+ col[i] = col[i].reduce(bigoh)
+ col.append(pRational(p, ZZ(1), bigoh))
+ self._lattice[x_ref] = col
+ self._capped[x_ref] = capped
+
+ # We update history
+ if self._history is not None:
+ self._history.append(('add', None, walltime(tme)))
+
+ #def _set_precision(self, x, column={}):
+ # p = self._p
+ # x_ref = weakref.ref(x)
+ # index = len(self._lattice[x_ref]) - 1
+ # n = len(self._elements)
+ # col = n * [self._approx_zero]
+ # self._lattice[x_ref] = col
+ # self._absolute_precisions[x_ref] = min([ c.valuation() for c in col ])
+
+ def mark_for_deletion(self, ref):
+ """
+ Mark an element for deletion.
+
+ This function is not meant to be called manually.
+ It is automatically called by the garbage collection when
+ an element is collected.
+
+ INPUT:
+
+ - ``ref`` -- a weak reference to the destroyed element
+
+ NOTE::
+
+ This method does not update the precision lattice.
+ The actual update is performed when the method :meth:`del_elements`
+ is called. This is automatically done at the creation of a new
+ element but can be done manually as well.
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2, label='markdel')
+ sage: prec = R.precision()
+ sage: x = R(1,10)
+ sage: prec
+ Precision Lattice on 1 object (label: markdel)
+ sage: del x # indirect doctest: x is here marked for deletion
+ sage: prec
+ Precision Lattice on 1 object (label: markdel)
+ sage: prec.del_elements() # x is indeed deleted
+ sage: prec
+ Precision Lattice on 0 object (label: markdel)
+ """
+ tme = walltime()
+ try:
+ index = len(self._lattice[ref]) - 1
+ except (IndexError, KeyError):
+ return
+ self._marked_for_deletion.append(index)
+ if self._history is not None:
+ self._history.append(('mark', index, walltime(tme)))
+
+ def del_elements(self, thresold=None):
+ """
+ Erase columns of the lattice precision matrix corresponding to
+ elements which are marked for deletion and reduce the matrix
+ in order to keep it upper triangular.
+
+ INPUT:
+
+ - ``thresold`` -- an integer or ``None`` (default: ``None``):
+ a column whose distance to the right at greater than the
+ thresold is not erased
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2, label='delelts')
+ sage: prec = R.precision()
+
+ sage: x = R(1,10)
+ sage: prec
+ Precision Lattice on 1 object (label: delelts)
+ sage: prec.precision_lattice()
+ [1024]
+
+ sage: del x
+ sage: prec
+ Precision Lattice on 1 object (label: delelts)
+ sage: prec.precision_lattice()
+ [1024]
+
+ sage: prec.del_elements()
+ sage: prec
+ Precision Lattice on 0 object (label: delelts)
+ sage: prec.precision_lattice()
+ []
+ """
+ p = self._p
+ n = len(self._elements)
+
+ self._marked_for_deletion.sort(reverse=True)
+ count = 0
+ for index in self._marked_for_deletion:
+ if thresold is not None and index < n - thresold: break
+ n -= 1; count += 1
+
+ tme = walltime()
+ del self._lattice[self._elements[index]]
+ del self._elements[index]
+
+ # Now, we echelonize
+ for i in range(index,n):
+ col = self._lattice[self._elements[i]]
+ vali = col[i].valuation()
+ valj = col[i+1].valuation()
+ d, u, v = col[i].xgcd(col[i+1])
+ up, vp = col[i+1]/d, col[i]/d
+ col[i] = d
+ del col[i+1]
+ for j in range(i+1,n):
+ col = self._lattice[self._elements[j]]
+ col[i], col[i+1] = u*col[i] + v*col[i+1], up*col[i] - vp*col[i+1]
+
+ # We update history
+ if self._history is not None:
+ self._history.append(('del', index, walltime(tme)))
+
+ # And we reduce a bit
+ # (we do not perform a complete reduction because it is costly)
+ self.reduce(index, partial=True)
+
+ del self._marked_for_deletion[:count]
+
+ def lift_to_precision(self, x, prec):
+ """
+ Lift the specified element to the specified precision
+
+ INPUT:
+
+ - ``x`` -- the element whose precision has to be lifted
+
+ - ``prec`` -- the new precision
+
+ NOTE:
+
+ The new precision lattice is computed as the intersection
+ of the current precision lattice with the subspace
+
+ ..MATH::
+
+ p^{prec} \Z_p dx \oplus \bigoplus_{y \neq x} \Q_p dy
+
+ This function may change at the same time the precision of
+ other elements having the same parent.
+
+ NOTE:
+
+ This function is not meant to be called directly.
+ You should prefer call the method :meth:`lift_to_precision`
+ of ``x`` instead.
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2)
+ sage: x = R(1,10); x
+ 1 + O(2^10)
+ sage: y = R(1,5); y
+ 1 + O(2^5)
+ sage: z = x + y; z
+ 2 + O(2^5)
+
+ sage: prec = R.precision()
+ sage: prec.lift_to_precision(z, 12)
+ sage: z
+ 2 + O(2^12)
+ sage: y
+ 1 + O(2^10)
+ """
+ ref = weakref.ref(x)
+ col = self._lattice[ref]
+ n = len(self._elements)
+
+ rows_by_val = { }
+ for i in range(len(col)):
+ v = col[i].valuation()
+ if v >= prec: continue
+ if rows_by_val.has_key(v):
+ rows_by_val[v].append(i)
+ else:
+ rows_by_val[v] = [i]
+ vals = rows_by_val.keys()
+ vals.sort()
+ vals.append(prec)
+
+ for t in range(len(vals)-1):
+ v, w = vals[t], vals[t+1]
+ rows = rows_by_val[v]
+ piv = max(rows)
+ alpha = col[piv].unit_part()
+ for i in rows:
+ if i == piv: continue
+ # We clear the entry on the i-th line
+ beta = col[i].unit_part()
+ for j in range(piv,n):
+ col_cur = self._lattice[self._elements[j]]
+ col_cur[i] = alpha*col_cur[i] - beta*col_cur[piv]
+ # We rescale the piv-th line
+ for j in range(piv,n):
+ col_cur = self._lattice[self._elements[j]]
+ col_cur[piv] = col_cur[piv] << (w-v)
+ # Now the entry on the piv-th line has valuation w
+ # We update the dictionary accordingly
+ if w < prec:
+ rows_by_val[w].append(piv)
+
+ # We clear the cached absolute precisions
+ self._absolute_precisions = { }
+
+
+ def _compute_precision_absolute(self, ref):
+ """
+ Compute the absolute precision of the given element and cache it
+
+ For internal use.
+ """
+ col = self._lattice[ref]
+ absprec = Infinity
+ capped = False
+ for i in range(len(col)):
+ v = col[i].valuation()
+ if v < absprec:
+ absprec = v
+ capped = self._capped[self._elements[i]]
+ elif v == absprec:
+ capped = capped and self._capped[self._elements[i]]
+ self._absolute_precisions[ref] = [absprec, capped]
+
+ def precision_absolute(self, x):
+ """
+ Return the absolute precision of the given element
+
+ INPUT:
+
+ - ``x`` -- the element whose absolute precision is requested
+
+ NOTE:
+
+ The absolute precision is obtained by projecting the precision
+ lattice onto the line of coordinate ``dx``
+
+ NOTE:
+
+ This function is not meant to be called directly.
+ You should prefer call the method :meth:`precision_absolute`
+ of ``x`` instead.
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2)
+ sage: x = R(1,10); x
+ 1 + O(2^10)
+ sage: y = R(1,5); y
+ 1 + O(2^5)
+ sage: z = x + y; z
+ 2 + O(2^5)
+ sage: z.precision_absolute() # indirect doctest
+ 5
+ """
+ ref = weakref.ref(x)
+ if not self._absolute_precisions.has_key(ref):
+ self._compute_precision_absolute(ref)
+ return self._absolute_precisions[ref][0]
+
+ def is_precision_capped(self, x):
+ """
+ Return whether the absolute precision on the given
+ results from a cap coming from the parent
+
+ INPUT:
+
+ - ``x`` -- the element
+
+ This function is not meant to be called directly.
+ You should prefer call the method :meth:`is_precision_capped`
+ of ``x`` instead.
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2)
+ sage: x = R(1,10); x
+ 1 + O(2^10)
+ sage: x.is_precision_capped() # indirect doctest
+ False
+
+ sage: y = x-x; y
+ O(2^40)
+ sage: y.is_precision_capped() # indirect doctest
+ True
+ """
+ ref = weakref.ref(x)
+ if not self._absolute_precisions.has_key(ref):
+ self._compute_precision_absolute(ref)
+ return self._absolute_precisions[ref][1]
+
+ def precision_lattice(self, elements=None, echelon=True):
+ """
+ Return a matrix representing the precision lattice on a
+ subset of elements.
+
+ INPUT:
+
+ - ``elements`` -- a list of elements or ``None`` (default: ``None``)
+
+ - ``echelon`` -- a boolean (default: ``True``); specify whether
+ the result should be in echelon form
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2, label='preclattice')
+ sage: prec = R.precision()
+ sage: x = R(1,10); y = R(1,5)
+ sage: u = x + y
+ sage: v = x - y
+ sage: prec.precision_lattice()
+ [ 1024 0 1024 1024]
+ [ 0 32 32 1099511627744]
+ [ 0 0 2097152 0]
+ [ 0 0 0 1099511627776]
+ sage: prec.precision_lattice([u,v])
+ [ 32 2016]
+ [ 0 2048]
+
+ Here is another example with matrices::
+
+ sage: M = matrix(R, 2, 2, [R(3,5),R(7,5),R(1,5),R(11,1)])
+ sage: N = M^10
+ sage: prec.precision_lattice()
+ 23 x 23 dense matrix over Integer Ring (use the '.str()' method to see the entries)
+
+ The next syntax provides as easy way to select an interesting
+ subset of variables (the selected subset consists of the four
+ entries of the matrix ``N``)::
+
+ sage: prec.precision_lattice(N)
+ [ 2048 512 28160 230400]
+ [ 0 2048 14336 258048]
+ [ 0 0 65536 65536]
+ [ 0 0 0 262144]
+
+ We can give a list of matrices as well::
+
+ sage: prec.precision_lattice([M,N])
+ [ 32 0 0 0 226115584 96788480 52174848 82804736]
+ [ 0 32 0 0 52174848 121765888 11829248 28516352]
+ [ 0 0 32 0 96788480 42762240 121765888 199614464]
+ [ 0 0 0 2 5175296 12475904 1782272 4045824]
+ [ 0 0 0 0 268435456 0 0 0]
+ [ 0 0 0 0 0 268435456 0 0]
+ [ 0 0 0 0 0 0 268435456 0]
+ [ 0 0 0 0 0 0 0 268435456]
+ """
+ if elements is None:
+ elements = self._elements
+ else:
+ elements = list_of_padics(elements)
+ n = len(self._elements)
+ rows = [ ]; val = 0
+ for ref in elements:
+ col = self._lattice[ref]
+ row = [ x.value() for x in col ]
+ valcol = min([ x.valuation() for x in col ])
+ if valcol < val: val = valcol
+ row += (n-len(row)) * [ZZ(0)]
+ rows.append(row)
+ from sage.matrix.constructor import matrix
+ M = matrix(rows).transpose()
+ if val < 0:
+ M *= self._p ** (-val)
+ if echelon:
+ M = M.change_ring(ZZ)
+ M.echelonize()
+ n = len(elements)
+ M = M.submatrix(0,0,n,n)
+ if val < 0:
+ M *= self._p ** val
+ return M
+
+ def number_of_diffused_digits(self, elements=None):
+ """
+ Return the number of diffused digits of precision within a
+ subset of elements.
+
+ NOTE:
+
+ A diffused digit of precision is a known digit which is not
+ located on a single variable but only atppears on a suitable
+ linear combinaison of variables.
+
+ The number of diffused digits of precision quantifies the
+ quality of the approximation of the lattice precision by a
+ jagged precision (that is a precision which is split over
+ all variables).
+
+ We refer to [] for a detail exposition of the notion of
+ diffused digits.
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2)
+ sage: prec = R.precision()
+ sage: x = R(1,10); y = R(1,5)
+ sage: u = x + y
+ sage: v = x - y
+
+ sage: prec.number_of_diffused_digits([x,y])
+ 0
+ sage: prec.number_of_diffused_digits([u,v])
+ 6
+
+ The elements `u` and `v` are known at absolute precision `O(2^5)`.
+ However, the sum `u + v = 2x` is known at precision `O(2^11)`, that
+ is with `6` more digits.
+ That is where the `6` diffused digits of precision comes from.
+
+ Here is another example with matrices::
+
+ sage: M = matrix(R, 2, 2, [R(3,5),R(7,5),R(1,5),R(11,1)])
+ sage: N = M^10
+
+ The next syntax provides as easy way to select an interesting
+ subset of variables (the selected subset consists of the four
+ entries of the matrix ``N``)::
+
+ sage: prec.number_of_diffused_digits(N)
+ 17
+ """
+ M = self.precision_lattice(elements)
+ n = M.nrows()
+ p = self._p
+ diffused = 0
+ for j in range(n):
+ val = minval = M[j,j].valuation(p)
+ for i in range(j):
+ v = M[i,j].valuation(p)
+ if v < minval: minval = v
+ diffused += val - minval
+ return diffused
+
+ def number_of_tracked_elements(self, dead=True):
+ """
+ Return the number of tracked elements through this precision
+ lattice
+
+ INPUT:
+
+ - ``dead`` -- a boolean (default: ``True``); whether dead
+ elements for which the corresponding column is still not
+ erased should be counted or not
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2, label='count')
+ sage: prec = R.precision()
+ sage: x = R(1,10); y = R(1,5)
+ sage: prec.number_of_tracked_elements()
+ 2
+
+ sage: u = x + y
+ sage: v = x - y
+ sage: prec.number_of_tracked_elements()
+ 4
+
+ sage: del x; del y
+ sage: prec.number_of_tracked_elements()
+ 4
+ sage: prec.number_of_tracked_elements(dead=False)
+ 2
+
+ sage: prec.del_elements()
+ sage: prec.number_of_tracked_elements()
+ 2
+ """
+ if dead:
+ return len(self._elements)
+ else:
+ count = 0
+ for x_ref in self._elements:
+ if x_ref() is not None: count += 1
+ return count
+
+
+ def tracked_elements(self, values=True, dead=True):
+ """
+ Return the list of tracked elements
+
+ INPUT:
+
+ - ``values`` -- a boolean (default: ``True``); if false,
+ the method returns a list of weak references on tracked
+ elements instead
+
+ - ``dead`` -- a boolean (default: ``True``); whether dead
+ elements for which the corresponding column is still not
+ erased should be listed or not
+
+ EXAMPLES::
+
+ sage: R = ZpLC(2, label='tracked')
+ sage: prec = R.precision()
+ sage: x = R(1,10); y = R(1,5)
+ sage: prec.tracked_elements()
+ [1 + O(2^10), 1 + O(2^5)]
+ sage: prec