summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorDavid Roe <roed.math@gmail.com>2017-08-21 04:49:20 +0000
committerDavid Roe <roed.math@gmail.com>2017-08-21 04:49:20 +0000
commit90a79f74235c9136d6266d221a53f95f91aa7795 (patch)
tree0d6b433322ae87fea36ac4ecb2a33171058329c8
parentUpdated SageMath version to 8.1.beta2 (diff)
parentMerge branch 'u/caruso/lattice_precision' of git://trac.sagemath.org/sage int... (diff)
Merge branch 'u/roed/lattice_precision' of git://trac.sagemath.org/sage into t/23505/lattice_precision
-rw-r--r--src/sage/rings/padics/all.py4
-rw-r--r--src/sage/rings/padics/factory.py368
-rw-r--r--src/sage/rings/padics/lattice_precision.py567
-rw-r--r--src/sage/rings/padics/padic_base_generic.py11
-rw-r--r--src/sage/rings/padics/padic_base_leaves.py206
-rw-r--r--src/sage/rings/padics/padic_lattice_element.py203
-rw-r--r--src/sage/rings/padics/padic_printing.pyx7
7 files changed, 1338 insertions, 28 deletions
diff --git a/src/sage/rings/padics/all.py b/src/sage/rings/padics/all.py
index 809b4aa..0ca7080 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, ZpLP, ZqCR, ZqCA, ZqFM, ZqFP #, ZpL, ZqL
+from .factory import Qp, Qq, Qp as pAdicField, QpCR, QpFP, QpLP, 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 6034161..e07d955 100644
--- a/src/sage/rings/padics/factory.py
+++ b/src/sage/rings/padics/factory.py
@@ -31,8 +31,10 @@ from .padic_base_leaves import (pAdicRingCappedRelative,
pAdicRingCappedAbsolute,
pAdicRingFixedMod,
pAdicRingFloatingPoint,
+ pAdicRingLattice,
pAdicFieldCappedRelative,
- pAdicFieldFloatingPoint)
+ pAdicFieldFloatingPoint,
+ pAdicFieldLattice)
from . import padic_printing
######################################################
@@ -95,7 +97,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):
"""
This implements create_key for Zp and Qp: moving it here prevents code duplication.
@@ -105,9 +107,9 @@ 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
@@ -186,7 +188,7 @@ 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)))
@@ -537,7 +539,7 @@ 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):
"""
Creates a key from input parameters for ``Qp``.
@@ -546,7 +548,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 +558,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'], label)
def create_object(self, version, key):
"""
@@ -566,7 +568,7 @@ class Qp_class(UniqueFactory):
TESTS::
- sage: Qp.create_object((3,4,2),(5, 41, 'capped-rel', 'series', '5', True, '|', (), -1))
+ sage: Qp.create_object((3,4,2),(5, 41, 'capped-rel', 'series', '5', True, '|', (), -1, None))
5-adic Field with capped relative precision 41
"""
if version[0] < 3 or (version[0] == 3 and version[1] < 2) or (version[0] == 3 and version[1] == 2 and version[2] < 3):
@@ -576,7 +578,7 @@ class Qp_class(UniqueFactory):
p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms = key
show_prec = 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 +594,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 +610,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':
+ 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")
@@ -1173,6 +1182,21 @@ def QpFP(p, prec = None, *args, **kwds):
"""
return Qp(p, prec, 'floating-point', *args, **kwds)
+def QpLP(p, prec = None, *args, **kwds):
+ """
+ A shortcut function to create `p`-adic fields with lattice precision.
+
+ See :func:`ZpLP` for more information about this model of precision.
+
+ EXAMPLES::
+
+ sage: R = QpLP(2)
+ sage: R
+ 2-adic Field with lattice precision
+
+ """
+ return Qp(p, prec, 'lattice', *args, **kwds)
+
def QqCR(q, prec = None, *args, **kwds):
"""
A shortcut function to create capped relative unramified `p`-adic
@@ -1595,7 +1619,7 @@ 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):
"""
Creates a key from input parameters for ``Zp``.
@@ -1604,9 +1628,9 @@ 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.
@@ -1617,7 +1641,7 @@ 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'], label=label)
def create_object(self, version, key):
"""
@@ -1627,7 +1651,7 @@ class Zp_class(UniqueFactory):
TESTS::
- sage: Zp.create_object((3,4,2),(5, 41, 'capped-rel', 'series', '5', True, '|', (), -1))
+ sage: Zp.create_object((3,4,2),(5, 41, 'capped-rel', 'series', '5', True, '|', (), -1, None))
5-adic Ring with capped relative precision 41
"""
if (version[0] < 3 or (len(version) > 1 and version[0] == 3 and version[1] < 2) or
@@ -1637,8 +1661,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.")
@@ -1647,14 +1672,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'])
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)
@@ -1667,6 +1692,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':
+ 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")
@@ -2281,6 +2309,308 @@ def ZpFP(p, prec = None, *args, **kwds):
"""
return Zp(p, prec, 'floating-point', *args, **kwds)
+def ZpLP(p, prec = None, *args, **kwds):
+ """
+ A shortcut function to create `p`-adic rings with lattice precision.
+
+ Below is a small demo of the features by this model of precision::
+
+ sage: R = ZpLP(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 = ZpLP(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 = ZpLP(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 = ZpLP(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)
+
+ It is these diffused digits of precision (which are
+ tracked but do not appear on the printing) that 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 = ZpLP(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', *args, **kwds)
+
def ZqCR(q, prec = None, *args, **kwds):
"""
A shortcut function to create capped relative unramified `p`-adic rings.
diff --git a/src/sage/rings/padics/lattice_precision.py b/src/sage/rings/padics/lattice_precision.py
new file mode 100644
index 00000000..1808afe
--- /dev/null
+++ b/src/sage/rings/padics/lattice_precision.py
@@ -0,0 +1,567 @@
+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):
+ 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):
+ status = ''.join(status)
+ if timings:
+ if tme < 0.000001:
+ s = " --- "
+ elif tme >= 10:
+ s = " >= 10s "
+ else:
+ s = "%.6fs" % tme
+ return s + " " + status
+ else:
+ return status
+
+class PrecisionLattice(UniqueRepresentation, SageObject):
+ # 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._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 self._label
+
+ def _repr_(self):
+ 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 self._p
+
+ def reduce(self, index=0, partial=False):
+ 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'):
+ # 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
+
+ # We compute the absolute precision of the new element and cache it
+ self._absolute_precisions[x_ref] = min([ c.valuation() for c in col ])
+
+ # 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):
+ tme = walltime()
+ try:
+ index = len(self._lattice[ref]) - 1
+ except IndexError:
+ 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):
+ 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):
+ 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)
+
+ def precision_absolute(self, x):
+ ref = weakref.ref(x)
+ return self._absolute_precisions[ref]
+
+ def precision_lattice(self, elements=None, echelon=True):
+ 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):
+ 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):
+ return len(self._elements)
+
+ def tracked_elements(self, values=True):
+ if values:
+ return [ ref() for ref in self._elements ]
+ else:
+ return list(self._elements)
+
+
+ # History
+
+ def history_enable(self):
+ if self._history is None:
+ self._history_init = ( len(self._elements), list(self._marked_for_deletion) )
+ self._history = [ ]
+
+ def history_disable(self):
+ self._history = self._history_init = None
+
+ def history_clear(self):
+ if self._history is None:
+ raise ValueError("History is not tracked")
+ self._history_init = ( len(self._elements), list(self._marked_for_deletion) )
+ self._history = [ ]
+
+ def history(self, compact=True, include_reduce=False, timings=True, output_type='asciiart'):
+ if self._history is None:
+ raise ValueError("History is not tracked")
+ total_time = 0
+ if output_type == 'asciiart':
+ # Legend:
+ # o : tracked element
+ # ~ : element marked for deletion
+ # r : partial reduction
+ # R : full Hermite reduction
+ (n, mark) = self._history_init
+ status = n*['o']
+ for index in mark:
+ status[index] = '~'
+ hist = [ ]; oldevent = ''
+ for (event, index, tme) in self._history:
+ if not include_reduce and (event == 'partial reduce' or event == 'full reduce'):
+ continue
+ if not compact or oldevent != event:
+ hist.append(format_history(total_time, status, timings))
+ total_time = 0
+ oldevent = event
+ total_time += tme
+ if event == 'add':
+ status.append('o')
+ elif event == 'mark':
+ status[index] = '~'
+ elif event == 'del':
+ del status[index]
+ elif event == 'partial reduce':
+ status_red = status[:index] + (len(status) - index)*['r']
+ hist.append(format_history(0, status_red, timings))
+ elif event == 'full reduce':
+ status_red = status[:index] + (len(status) - index)*['R']
+ hist.append(format_history(0, status_red, timings))
+ else:
+ raise RuntimeError("event not known in the history")
+ hist.append(format_history(total_time, status, timings))
+ return '\n'.join(hist)
+ else:
+ raise NotImplementedError
+
+ def timings(self, action=None):
+ if self._history is None:
+ raise ValueError("History is not tracked")
+ tme_by_event = { }
+ for (event, _, tme) in self._history:
+ if tme_by_event.has_key(event):
+ tme_by_event[event] += tme
+ else:
+ tme_by_event[event] = tme
+ if action is None:
+ return tme_by_event
+ if tme_by_event.has_key(action):
+ return tme_by_event[action]
+ else:
+ return 0
diff --git a/src/sage/rings/padics/padic_base_generic.py b/src/sage/rings/padics/padic_base_generic.py
index eebbc46..e5461cd 100644
--- a/src/sage/rings/padics/padic_base_generic.py
+++ b/src/sage/rings/padics/padic_base_generic.py
@@ -20,6 +20,9 @@ from __future__ import absolute_import
# http://www.gnu.org/licenses/
#*****************************************************************************
+from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+
from .padic_generic import pAdicGeneric
from .misc import precprint
from sage.rings.padics.pow_computer import PowComputer
@@ -44,9 +47,12 @@ class pAdicBaseGeneric(pAdicGeneric):
if self.is_capped_relative():
coerce_list = [pAdicCoercion_ZZ_CR(self), pAdicCoercion_QQ_CR(self)]
convert_list = []
- else:
+ elif self.is_floating_point():
coerce_list = [pAdicCoercion_ZZ_FP(self), pAdicCoercion_QQ_FP(self)]
convert_list = []
+ else:
+ coerce_list = [QQ]
+ convert_list = []
elif self.is_capped_relative():
coerce_list = [pAdicCoercion_ZZ_CR(self)]
convert_list = [pAdicConvert_QQ_CR(self)]
@@ -60,7 +66,8 @@ class pAdicBaseGeneric(pAdicGeneric):
coerce_list = [pAdicCoercion_ZZ_FP(self)]
convert_list = [pAdicConvert_QQ_FP(self)]
else:
- raise RuntimeError
+ coerce_list = [ZZ]
+ convert_list = []
self._populate_coercion_lists_(coerce_list=coerce_list, convert_list=convert_list, element_constructor=element_class)
def _repr_(self, do_latex=False):
diff --git a/src/sage/rings/padics/padic_base_leaves.py b/src/sage/rings/padics/padic_base_leaves.py
index 702991a..30317bd 100644
--- a/src/sage/rings/padics/padic_base_leaves.py
+++ b/src/sage/rings/padics/padic_base_leaves.py
@@ -197,12 +197,19 @@ from .generic_nodes import pAdicFieldBaseGeneric, \
pAdicFixedModRingGeneric, \
pAdicCappedAbsoluteRingGeneric, \
pAdicFloatingPointRingGeneric, \
- pAdicFloatingPointFieldGeneric
+ pAdicFloatingPointFieldGeneric, \
+ pAdicGeneric
from .padic_capped_relative_element import pAdicCappedRelativeElement
from .padic_capped_absolute_element import pAdicCappedAbsoluteElement
from .padic_fixed_mod_element import pAdicFixedModElement
from .padic_floating_point_element import pAdicFloatingPointElement
+from .padic_lattice_element import pAdicLatticeElement
+from .lattice_precision import PrecisionLattice
+
from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+from sage.rings.infinity import Infinity
+
class pAdicRingCappedRelative(pAdicRingBaseGeneric, pAdicCappedRelativeRingGeneric):
r"""
@@ -719,3 +726,200 @@ class pAdicFieldFloatingPoint(pAdicFieldBaseGeneric, pAdicFloatingPointFieldGene
return True
elif R.precision_cap() == self.precision_cap() and self._printer.richcmp_modes(R._printer, op_LE):
return True
+
+ def _repr_(self, do_latex=False):
+ r"""
+ Print representation.
+
+ EXAMPLES::
+
+ sage: K = QpFP(17); K #indirect doctest
+ 17-adic Field with floating precision 20
+ sage: latex(K)
+ \QQ_{17}
+ """
+ if do_latex:
+ return "\\QQ_{%s}" % self.prime()
+ return "%s-adic Field with floating precision %s"%(self.prime(), self.precision_cap())
+
+
+
+# Lattice precision
+###################
+
+# Maybe the next class should go to sage.rings.padics.generic_nodes but I
+# don't understand quite well the structure of all classes in this directory
+class pAdicLatticeGeneric(pAdicGeneric):
+ def __init__(self, p, prec, label=None, proof=False):
+ if proof:
+ raise NotImplementedError("p-adic with *proved* lattice precision not implemented yet")
+ if label is None:
+ self._label = None
+ else:
+ self._label = str(label)
+ self._prec_cap = prec
+ self._precision = PrecisionLattice(p, label)
+ # We do not use the standard attribute element_class
+ # because we need to be careful with precision
+ # Instead we implement _element_constructor_ (cf below)
+ self._element_class = pAdicLatticeElement
+
+ def _prec_type(self):
+ return 'lattice'
+
+ def precision(self):
+ return self._precision
+
+ def label(self):
+ return self._label
+
+ def _element_constructor_(self, x, prec=None):
+ # We first try the __copy__ method which is sharp on precision
+ try:
+ if prec is None:
+ return x.__copy__(parent=self)
+ else:
+ return x.__copy__(parent=self).add_bigoh(prec)
+ except (TypeError, ValueError, AttributeError):
+ pass
+ return self._element_class(self, x, prec)
+
+ def convert_multiple(self, *elts):
+ p = self.prime()
+
+ # We sort elements by precision lattice
+ elt_by_prec = { }
+ elt_other = [ ]
+ indices = { }
+ for i in range(len(elts)):
+ x = elts[i]; idx = id(x)
+ if indices.has_key(idx):
+ indices[idx].append(i)
+ else:
+ indices[idx] = [i]
+ if isinstance(x, pAdicLatticeElement):
+ prec = x.parent().precision()
+ if prec.prime() != p:
+ raise TypeError("conversion between different p-adic rings not supported")
+ if elt_by_prec.has_key(prec):
+ elt_by_prec[prec].append(x)
+ else:
+ elt_by_prec[prec] = [x]
+ else:
+ elt_other.append(x)
+
+ # We create the elements
+ ans = len(elts)*[None]
+ selfprec = self._precision
+ # First the elements with precision lattice
+ for (prec, L) in elt_by_prec.iteritems():
+ if prec is selfprec:
+ # Here, we use the __copy__ method in order
+ # to be sharp on precision
+ for x in L:
+ y = x.__copy__(parent=self)
+ for i in indices[id(x)]:
+ ans[i] = y
+ else:
+ lattice = prec.precision_lattice(L)
+ for j in range(len(L)):
+ x = L[j]; dx = [ ]
+ for i in range(j):
+ dx.append([L[i], lattice[i,j]])
+ prec = lattice[j,j].valuation(p)
+ y = self._element_class(self, x.value(), prec, dx=dx, dx_mode='values', check=False, reduce=False)
+ for i in indices[id(x)]:
+ ans[i] = y
+ L[j] = y
+ # Now the other elements
+ for x in elt_other:
+ y = self._element_class(self, x)
+ for i in indices[id(x)]:
+ ans[i] = y
+
+ # We return the created elements
+ return ans
+
+
+
+class pAdicRingLattice(pAdicLatticeGeneric, pAdicRingBaseGeneric):
+ def __init__(self, p, prec, print_mode, name, label=None, proof=False):
+ pAdicLatticeGeneric.__init__(self, p, prec, label, proof)
+ pAdicRingBaseGeneric.__init__(self, p, prec, print_mode, str(p), None)
+
+ def _repr_(self, do_latex=False):
+ if do_latex:
+ if self._label is not None:
+ return "\\verb'%s' (\simeq \\mathbb Z_{%s})" % (self._label, self.prime())
+ else:
+ return "\\mathbb Z_{%s}" % self.prime()
+ else:
+ if self._label is not None:
+ return "%s-adic Ring with lattice precision (label: %s)" % (self.prime(), self._label)
+ else:
+ return "%s-adic Ring with lattice precision" % self.prime()
+
+ def _coerce_map_from_(self, R):
+ if isinstance(R, pAdicRingLattice) and R.precision() is self.precision():
+ return True
+
+ def random_element(self, prec=None):
+ if prec is None:
+ prec = self._prec_cap
+ p = self.prime()
+ x = ZZ.random_element(p**prec)
+ return self._element_class(self, x, prec=prec)
+
+ def integer_ring(self, print_mode=None):
+ if print_mode is None:
+ return self
+ from sage.rings.padics.factory import Zp
+ return Zp(self.prime(), self.precision_cap(), 'lattice', print_mode=self._modified_print_mode(print_mode),
+ names=self._uniformizer_print(), label=self._label)
+
+ def fraction_field(self, print_mode=None):
+ from sage.rings.padics.factory import Qp
+ return Qp(self.prime(), self.precision_cap(), 'lattice', print_mode=self._modified_print_mode(print_mode),
+ names=self._uniformizer_print(), label=self._label)
+
+
+
+class pAdicFieldLattice(pAdicLatticeGeneric, pAdicFieldBaseGeneric):
+ def __init__(self, p, prec, print_mode, name, label=None, proof=False):
+ pAdicLatticeGeneric.__init__(self, p, prec, label, proof)
+ pAdicFieldBaseGeneric.__init__(self, p, prec, print_mode, str(p), None)
+
+ def _repr_(self, do_latex=False):
+ if do_latex:
+ if self._label is not None:
+ return "\\verb'%s' (\simeq \\mathbb Q_{%s})" % (self._label, self.prime())
+ else:
+ return "\\mathbb Q_{%s}" % self.prime()
+ else:
+ if self._label is not None:
+ return "%s-adic Field with lattice precision (label: %s)" % (self.prime(), self._label)
+ else:
+ return "%s-adic Field with lattice precision" % self.prime()
+
+ def _coerce_map_from_(self, R):
+ if isinstance(R, (pAdicRingLattice, pAdicFieldLattice)) and R.precision() is self.precision():
+ return True
+
+ def random_element(self, prec=None):
+ if prec is None:
+ prec = self._prec_cap
+ p = self.prime()
+ x = ZZ.random_element(p**prec)
+ return self._element_class(self, x, prec)
+
+ def integer_ring(self, print_mode=None):
+ from sage.rings.padics.factory import Zp
+ return Zp(self.prime(), self.precision_cap(), 'lattice', print_mode=self._modified_print_mode(print_mode),
+ names=self._uniformizer_print(), label=self._label)
+
+ def fraction_field(self, print_mode=None):
+ if print_mode is None:
+ return self
+ from sage.rings.padics.factory import Qp
+ return Qp(self.prime(), self.precision_cap(), 'lattice', print_mode=self._modified_print_mode(print_mode),
+ names=self._uniformizer_print(), label=self._label)
diff --git a/src/sage/rings/padics/padic_lattice_element.py b/src/sage/rings/padics/padic_lattice_element.py
new file mode 100644
index 00000000..289fe05
--- /dev/null
+++ b/src/sage/rings/padics/padic_lattice_element.py
@@ -0,0 +1,203 @@
+import _weakref as weakref
+
+from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+from sage.rings.infinity import Infinity
+
+from sage.rings.padics.generic_nodes import pAdicRingBaseGeneric
+from sage.rings.padics.padic_generic_element import pAdicGenericElement
+from sage.rings.padics.lattice_precision import pRational
+
+
+class pAdicLatticeElement(pAdicGenericElement):
+ def __init__(self, parent, x, prec=None, dx=[], dx_mode='linear_combinaison', valuation=None, check=True, reduce=True):
+ self._parent = parent
+ pAdicGenericElement.__init__(self, parent)
+ self._precision = parent.precision()
+ if check:
+ if isinstance(x, pAdicGenericElement):
+ if parent.prime() != x.parent().prime():
+ raise TypeError("conversion between different p-adic rings/fields not supported")
+ if prec is None:
+ prec = x.precision_absolute()
+ else:
+ prec = min(prec, x.precision_absolute())
+ if isinstance(parent, pAdicRingBaseGeneric):
+ x = ZZ(x)
+ else:
+ x = QQ(x)
+ cap = parent.precision_cap()
+ if prec is None or prec > cap:
+ prec = cap
+ self._precision.new_element(self, dx, bigoh=prec, dx_mode=dx_mode)
+ if isinstance(x, pRational):
+ self._value = x
+ else:
+ self._value = pRational(parent.prime(), QQ(x))
+ if reduce:
+ self._value = self._value.reduce(prec)
+ self._approx_one = pRational(parent.prime(), 1)
+ self._approx_minusone = pRational(parent.prime(), -1)
+
+ def __hash__(self):
+ return id(self)
+
+ def _is_base_elt(self, p):
+ return p == self._parent.prime()
+
+ def approximation(self):
+ prec = self.precision_absolute()
+ app = self._value.reduce(prec)
+ return app.value()
+
+ def value(self):
+ return self._value.value()
+
+ def precision_lattice(self):
+ return self._precision
+
+ def precision_absolute(self):
+ return self._precision.precision_absolute(self)
+
+ def valuation(self, secure=False):
+ p = self._parent.prime()
+ val = self._value.valuation()
+ prec = self.precision_absolute()
+ if val < prec:
+ return val
+ elif secure:
+ raise PrecisionError("Not enough precision")
+ else:
+ return prec
+
+ def precision_relative(self, secure=False):
+ return self.precision_absolute() - self.valuation(secure=secure)
+
+ def _cmp_(self, other):
+ if (self-other).is_zero():
+ return 0
+ else:
+ return cmp(self.lift(), other.lift())
+
+ def _add_(self, other):
+ x = self._value + other._value
+ dx = [ [self, self._approx_one],
+ [other, self._approx_one] ]
+ return self.__class__(self._parent, x, self._parent.precision_cap(), dx=dx, check=False)
+
+ def _sub_(self, other):
+ x = self._value - other._value
+ dx = [ [self, self._approx_one],
+ [other, self._approx_minusone] ]
+ return self.__class__(self._parent, x, self._parent.precision_cap(), dx=dx, check=False)
+
+ def _mul_(self, other):
+ x_self = self._value
+ x_other = other._value
+ x = x_self * x_other
+ dx = [ [self, x_other],
+ [other, x_self ] ]
+ return self.__class__(self._parent, x, self._parent.precision_cap(), dx=dx, check=False)
+
+ def _div_(self, other):
+ if other.is_zero():
+ raise PrecisionError("cannot divide by something indistinguishable from zero")
+ p = self._parent.prime()
+ cap = self._parent.precision_cap()
+ x_self = self._value
+ x_other = other._value
+ x = x_self / x_other
+ # dx = (1/other)*dself - (self/other^2)*dother
+ dx = [ [self, self._approx_one/x_other],
+ [other, -x_self/(x_other*x_other)] ]
+ return self.__class__(self._parent.fraction_field(), x, self._parent.precision_cap(), dx=dx, check=False)
+
+ def __invert__(self):
+ if self.is_zero():
+ raise PrecisionError("cannot invert something indistinguishable from zero")
+ p = self._parent.prime()
+ cap = self._parent.precision_cap()
+ x_self = self._value
+ x = self._approx_one / x_self
+ # dx = -(1/self^2)*dself
+ dx = [ [self, self._approx_minusone/(x_self*x_self)] ]
+ return self.__class__(self._parent.fraction_field(), x, self._parent.precision_cap(), dx=dx, check=False)
+
+ def add_bigoh(self, prec):
+ x = self._value
+ dx = [ [self, self._approx_one ] ]
+ return self.__class__(self._parent, x, prec, dx=dx, check=False)
+
+ def lift_to_precision(self, prec=None, infer_precision=False):
+ from warnings import warn
+ warn("use lift_to_precision with extreme caution in the framework of lattice precision")
+ parent = self._parent
+ if prec is None:
+ prec = parent.precision_cap()
+ if infer_precision:
+ lift = self.__copy__()
+ parent.precision().lift_to_precision(lift, prec)
+ else:
+ if prec < self.precision_absolute():
+ prec = self.precision_absolute()
+ lift = self.__class__(parent, self._value, prec, check=False)
+ return lift
+
+ def _is_exact_zero(self):
+ return False
+
+ def _is_inexact_zero(self):
+ absprec = self.precision_absolute()
+ return self._value.valuation() >= absprec
+
+ def is_zero(self, prec=None):
+ absprec = self.precision_absolute()
+ if prec is None:
+ prec = absprec
+ else:
+ prec = min(absprec, prec)
+ return self._value.valuation() >= prec
+
+ def value(self):
+ return self._value
+
+ def lift(self):
+ return self._value.value()
+
+ def rational_reconstruction(self):
+ return QQ(self._value.value())
+
+ def __rshift__(self, n):
+ return self << (-n)
+
+ def __lshift__(self, n):
+ p = self._parent.prime()
+ powp = pRational(p, ZZ(1), n)
+ x = self._value * powp
+ dx = [ [self, powp] ]
+ return self.__class__(self._parent, x, self._parent.precision_cap(), dx=dx, check=False)
+
+ def unit_part(self):
+ v = self.valuation(secure=True)
+ return self >> v
+
+ def val_unit(self):
+ v = self.valuation(secure=True)
+ return v, self >> v
+
+ def __copy__(self, parent=None):
+ if parent is None:
+ parent = self._parent
+ else:
+ if parent.precision() is not self._parent.precision():
+ raise TypeError("parent must share the same precision object")
+ if isinstance(parent, pAdicRingBaseGeneric) and self.valuation() < 0:
+ raise ValueError("element of negative valuation cannot be convert to the integer ring")
+ dx = [ [ self, self._approx_one ] ]
+ return self.__class__(parent, self._value, self._parent.precision_cap(), dx=dx, check=False)
+
+ def list(self, lift_mode='simple', start_val=None):
+ # TODO: implement other lift modes
+ p = self._parent.prime()
+ prec = self.precision_absolute()
+ return self._value.list(prec)
diff --git a/src/sage/rings/padics/padic_printing.pyx b/src/sage/rings/padics/padic_printing.pyx
index cae22ca..9a38203 100644
--- a/src/sage/rings/padics/padic_printing.pyx
+++ b/src/sage/rings/padics/padic_printing.pyx
@@ -1033,13 +1033,12 @@ cdef class pAdicPrinter_class(SageObject):
v = elt.valuation()
if v >= 0:
lift_z = <Integer> elt.lift()
+ pprec = self.prime_pow.pow_Integer(mpz_get_ui((<Integer>elt.precision_absolute()).value))
else:
lift_z = <Integer> elt.unit_part().lift()
+ pprec = self.prime_pow.pow_Integer(mpz_get_ui((<Integer>elt.precision_relative()).value))
+ mpz_mod(lift_z.value, lift_z.value, pprec.value)
if not pos:
- if v >= 0:
- pprec = self.prime_pow.pow_Integer(mpz_get_ui((<Integer>elt.precision_absolute()).value))
- else:
- pprec = self.prime_pow.pow_Integer(mpz_get_ui((<Integer>elt.precision_relative()).value))
if lift_z > pprec / 2:
mpz_sub(lift_z.value, lift_z.value, pprec.value)
else: