summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorXavier Caruso <xavier.caruso@univ-rennes1.fr>2017-08-02 01:43:20 +0200
committerXavier Caruso <xavier.caruso@univ-rennes1.fr>2017-08-02 01:43:20 +0200
commit858492b4b0cd87dfb63fba0f3ad239801b98cb62 (patch)
tree23edca6d6222719b80a422f9cdd497dd368a78ee
parentFirst more-or-less working implementation (diff)
Second rough implementation of lattice precision
-rw-r--r--src/sage/rings/padics/all.py6
-rw-r--r--src/sage/rings/padics/factory.py393
-rw-r--r--src/sage/rings/padics/lattice_precision.py832
-rw-r--r--src/sage/rings/padics/padic_base_generic.py6
-rw-r--r--src/sage/rings/padics/padic_base_leaves.py191
-rw-r--r--src/sage/rings/padics/padic_lattice_element.py201
-rw-r--r--src/sage/rings/padics/padic_printing.pyx7
7 files changed, 1312 insertions, 324 deletions
diff --git a/src/sage/rings/padics/all.py b/src/sage/rings/padics/all.py
index 9665553..0ca7080 100644
--- a/src/sage/rings/padics/all.py
+++ b/src/sage/rings/padics/all.py
@@ -1,11 +1,9 @@
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
from .pow_computer_ext import PowComputer_ext_maker
from .discrete_value_group import DiscreteValueGroup
-
-from .lattice_precision import ZpLP
diff --git a/src/sage/rings/padics/factory.py b/src/sage/rings/padics/factory.py
index 41fc26b..332bf60 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
######################################################
@@ -70,7 +72,7 @@ ext_table['u', pAdicFieldFloatingPoint] = UnramifiedExtensionFieldFloatingPoint
#ext_table['u', pAdicRingLazy] = UnramifiedExtensionRingLazy
-def get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, check, valid_non_lazy_types):
+def get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, check, valid_non_lazy_types, label=None):
"""
This implements create_key for Zp and Qp: moving it here prevents code duplication.
@@ -80,9 +82,18 @@ def get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, pr
sage: from sage.rings.padics.factory import get_key_base
sage: get_key_base(11, 5, 'capped-rel', None, 0, None, None, None, ':', None, None, True, ['capped-rel'])
- (11, 5, 'capped-rel', 'series', '11', True, '|', (), -1)
+ (11, 5, 'capped-rel', 'series', '11', True, '|', (), -1, None)
sage: get_key_base(12, 5, 'capped-rel', 'digits', 0, None, None, None, None, None, None, False, ['capped-rel'])
- (12, 5, 'capped-rel', 'digits', '12', True, '|', ('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B'), -1)
+ (12,
+ 5,
+ 'capped-rel',
+ 'digits',
+ '12',
+ True,
+ '|',
+ ('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B'),
+ -1,
+ None)
"""
if check:
if not isinstance(p, Integer):
@@ -157,9 +168,9 @@ def get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, pr
else:
name = str(names)
if type in valid_non_lazy_types:
- key = (p, prec, type, print_mode, name, print_pos, print_sep, tuple(print_alphabet), print_max_terms)
+ key = (p, prec, type, print_mode, name, print_pos, print_sep, tuple(print_alphabet), print_max_terms, label)
elif type == 'lazy':
- key = (p, prec, halt, print_mode, name, print_pos, print_sep, tuple(print_alphabet), print_max_terms)
+ key = (p, prec, halt, print_mode, name, print_pos, print_sep, tuple(print_alphabet), print_max_terms, label)
else:
print(type)
raise ValueError("type must be %s or lazy"%(", ".join(valid_non_lazy_types)))
@@ -471,7 +482,7 @@ class Qp_class(UniqueFactory):
"""
def create_key(self, p, prec = DEFAULT_PREC, type = 'capped-rel', print_mode = None,
halt = DEFAULT_HALT, names = None, ram_name = None, print_pos = None,
- print_sep = None, print_alphabet = None, print_max_terms = None, check = True):
+ print_sep = None, print_alphabet = None, print_max_terms = None, check = True, label = None):
"""
Creates a key from input parameters for ``Qp``.
@@ -480,9 +491,9 @@ class Qp_class(UniqueFactory):
TESTS::
sage: Qp.create_key(5,40)
- (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1)
+ (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1, None)
"""
- return get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, check, ['capped-rel', 'floating-point'])
+ return get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, print_sep, print_alphabet, print_max_terms, check, ['capped-rel', 'floating-point', 'lattice'], label)
def create_object(self, version, key):
"""
@@ -492,14 +503,14 @@ 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):
p, prec, type, print_mode, name = key
print_pos, print_sep, print_alphabet, print_max_terms = None, None, None, None
else:
- p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, label = key
if isinstance(type, Integer):
# lazy
raise NotImplementedError("lazy p-adics need more work. Sorry.")
@@ -515,7 +526,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 = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, label = key
if type == 'capped-rel':
if print_mode == 'terse':
@@ -531,6 +542,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': False}, 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': True}, 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': True}, name, label)
else:
raise ValueError("unexpected type")
@@ -1081,6 +1099,24 @@ def QpFP(p, prec = DEFAULT_PREC, print_mode = None, halt = DEFAULT_HALT, names =
print_pos=print_pos, print_sep=print_sep, print_alphabet=print_alphabet, print_max_terms=print_max_terms,
type = 'floating-point')
+def QpLP(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, label=None):
+ """
+ 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=p, prec=prec, print_mode=print_mode, halt=halt, check=check, names=names,
+ print_pos=print_pos, print_sep=print_sep, print_alphabet=print_alphabet, print_max_terms=print_max_terms,
+ type = 'lattice', label=label)
+
#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):
# """
@@ -1521,7 +1557,7 @@ class Zp_class(UniqueFactory):
"""
def create_key(self, p, prec = DEFAULT_PREC, type = 'capped-rel', print_mode = None, halt = DEFAULT_HALT,
names = None, ram_name = None, print_pos = None, print_sep = None, print_alphabet = None,
- print_max_terms = None, check = True):
+ print_max_terms = None, check = True, label = None):
"""
Creates a key from input parameters for ``Zp``.
@@ -1530,12 +1566,21 @@ class Zp_class(UniqueFactory):
TESTS::
sage: Zp.create_key(5,40)
- (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1)
+ (5, 40, 'capped-rel', 'series', '5', True, '|', (), -1, None)
sage: Zp.create_key(5,40,print_mode='digits')
- (5, 40, 'capped-rel', 'digits', '5', True, '|', ('0', '1', '2', '3', '4'), -1)
+ (5,
+ 40,
+ 'capped-rel',
+ 'digits',
+ '5',
+ True,
+ '|',
+ ('0', '1', '2', '3', '4'),
+ -1,
+ None)
"""
return get_key_base(p, prec, type, print_mode, halt, names, ram_name, print_pos, print_sep, print_alphabet,
- print_max_terms, check, ['capped-rel', 'fixed-mod', 'capped-abs', 'floating-point'])
+ print_max_terms, check, ['capped-rel', 'fixed-mod', 'capped-abs', 'floating-point', 'lattice'], label=label)
def create_object(self, version, key):
"""
@@ -1545,7 +1590,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
@@ -1553,7 +1598,7 @@ class Zp_class(UniqueFactory):
p, prec, type, print_mode, name = key
print_pos, print_sep, print_alphabet, print_max_terms = None, None, None, None
else:
- p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, label = key
if isinstance(type, Integer):
# lazy
raise NotImplementedError("lazy p-adics need more work. Sorry.")
@@ -1562,14 +1607,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, 0, name, None, print_pos, print_sep, print_alphabet,
- print_max_terms, False, ['capped-rel', 'fixed-mod', 'capped-abs'])
+ print_max_terms, 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 = key
+ p, prec, type, print_mode, name, print_pos, print_sep, print_alphabet, print_max_terms, 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}, name)
@@ -1582,6 +1627,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': False}, 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': True}, name, label)
else:
raise ValueError("unexpected type")
@@ -2188,6 +2236,311 @@ def ZpFP(p, prec = DEFAULT_PREC, print_mode = None, halt = DEFAULT_HALT, names =
print_pos=print_pos, print_sep=print_sep, print_alphabet=print_alphabet, print_max_terms=print_max_terms,
type = 'floating-point')
+def ZpLP(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, label=None):
+ """
+ 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, obseve how the precision lattices 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)` when `n`
+ is the number of tracked elements.
+
+ - The destruction of one element has a cost `O(m^2)` when
+ `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. squares matrices of size 5 or
+ polynomials of degree 20 are accessible).
+
+ The class :class:`PrecisionLattice` provides several
+ features for introspection (especially concerning timings).
+ If enables, 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=p, prec=prec, print_mode=print_mode, halt=halt, check=check, names=names,
+ print_pos=print_pos, print_sep=print_sep, print_alphabet=print_alphabet, print_max_terms=print_max_terms,
+ type = 'lattice', label=label)
+
#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):
# """
diff --git a/src/sage/rings/padics/lattice_precision.py b/src/sage/rings/padics/lattice_precision.py
index 640ca37..1808afe 100644
--- a/src/sage/rings/padics/lattice_precision.py
+++ b/src/sage/rings/padics/lattice_precision.py
@@ -1,323 +1,567 @@
-"""
-Toy implementation of lattice precision for p-adic numbers.
-
-Here is a small demo::
-
- sage: R = ZpLP(3)
- 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)
- 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, 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)
- 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)
-"""
-
import _weakref as weakref
-from sage.misc.prandom import randint
+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.padics.generic_nodes import pAdicRingBaseGeneric
-from sage.rings.padics.padic_generic_element import pAdicGenericElement
+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)
-class pAdicRingLattice(pAdicRingBaseGeneric):
+ 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._prec_cap
+ # . self._cap
# a cap for the (working) precision
- # meaning that the precision lattice always contains p^(self._prec_cap)
+ # meaning that the precision lattice always contains p^(self._cap)
# . self._elements
# list of weak references of elements in this parent
- # . self._precision
+ # . self._lattice
# an upper triangular matrix over ZZ representing the
# lattice of precision
# (its columns are indexed by self._elements)
-
- def __init__(self, p, prec=50, print_mode={'mode':'terse'}, label=None):
- if label is None:
- self._label = None
- else:
- self._label = str(label)
+ # . self._absolute_precisions
+ def __init__(self, p, label):
+ self._p = p
+ self._label = label
self._elements = [ ]
- self._precision = { }
- self._prec_cap = prec
- self._modulus = p**prec
- pAdicRingBaseGeneric.__init__(self, p, prec, print_mode, '', pAdicLatticeElement)
-
- def _prec_type(self):
- return 'lattice'
+ 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 precision_cap(self):
- return self._prec_cap
-
def _repr_(self):
- s = "%s-adic Field with lattice precision" % self.prime()
- if self._label is not None:
- s += ", called '%s'" % self._label
- return s
+ 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 _new_element(self, x, prec, dx):
+ 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._del_element)
+ x_ref = weakref.ref(x, self.mark_for_deletion)
self._elements.append(x_ref)
- if prec is None or prec > self._prec_cap:
- prec = self._prec_cap
- modulus = self.prime()**prec
- col = n * [ZZ(0)]
- for ref,scalar in dx:
- c = self._precision[ref]
- for i in range(len(c)):
- col[i] += scalar * c[i]
+ 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] = ZZ(col[i]) % modulus
- # We should handle the case where col[i] is not an integer but I'm lazy...
- col.append(modulus)
- self._precision[x_ref] = col
-
- def _del_element(self, ref):
+ 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._precision[ref]) - 1
+ index = len(self._lattice[ref]) - 1
except IndexError:
return
- del self._elements[index]
- del self._precision[ref]
+ self._marked_for_deletion.append(index)
+ if self._history is not None:
+ self._history.append(('mark', index, walltime(tme)))
- # Now, we echelonize...
- p = self.prime()
+ def del_elements(self, thresold=None):
+ p = self._p
n = len(self._elements)
- modulus = self._modulus
- for i in range(index,n):
- col = self._precision[self._elements[i]]
- 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._precision[self._elements[j]]
- col[i], col[i+1] = (u*col[i] + v*col[i+1]) % modulus, (up*col[i] + vp*col[i+1]) % modulus
- # Actually, the matrix we obtained this way is not in Hermite normal form
- # because its entries above the diagonal are not necessarily completely reduced
- # (though they are at least modulo p^precision_cap)
- # Should we reduce them further (this affects seriously the complexity)?
- def precision_absolute(self, x):
+ 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)
- p = self.prime()
- return min([ c.valuation(p) for c in self._precision[ref] ])
-
- def working_precision(self, x):
- # Do something better here
- return self._prec_cap
-
-ZpLP = pAdicRingLattice
-
-
-# The elements
-##############
-
-class pAdicLatticeElement(pAdicGenericElement):
- # Internal variable:
- # . self._approximation
- # an approximation of this p-adic number
- # it is defined modulo p^(working_precision)
-
- def __init__(self, parent, x, prec=None, dx={}):
- self._parent = parent
- self._hash = hash((x, randint(1, 100000)))
- pAdicGenericElement.__init__(self, parent)
- if prec is None:
- prec = parent.precision_cap()
- self._parent._new_element(self, prec, dx)
- self._approximation = x
-
- def __hash__(self):
- return self._hash
-
- def working_precision(self):
- return self.parent().working_precision(self)
-
- def approximation(self, reduce=True):
- if reduce:
- workprec = self.working_precision()
- p = self._parent.prime()
- self._approximation %= p**workprec
- return self._approximation
-
- def precision_absolute(self):
- return self.parent().precision_absolute(self)
-
- def valuation(self, secure=False):
- p = self._parent.prime()
- val = self._approximation.valuation(p)
- prec = self.precision_absolute()
- if val < prec:
- return val
- elif secure:
- raise PrecisionError("Not enough precision")
- else:
- return prec
+ col = self._lattice[ref]
+ n = len(self._elements)
- def precision_relative(self, secure=False):
- return self.precision_absolute() - self.valuation(secure=secure)
+ 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 _repr_(self):
- prec = self.precision_absolute()
- p = self._parent.prime()
- app = self._approximation % (p**prec)
- return "%s + O(%s^%s)" % (app, p, prec)
-
- def _add_(self, other):
- x = self.approximation() + other.approximation()
- dx = [ [weakref.ref(self), 1],
- [weakref.ref(other), 1] ]
- return self.__class__(self._parent, x, dx=dx)
-
- def _sub_(self, other):
- x = self.approximation() - other.approximation()
- dx = [ [weakref.ref(self), 1 ],
- [weakref.ref(other), -1] ]
- return self.__class__(self._parent, x, dx=dx)
-
- def _mul_(self, other):
- x_self = self.approximation()
- x_other = other.approximation()
- x = x_self * x_other
- dx = [ [weakref.ref(self), x_other],
- [weakref.ref(other), x_self ] ]
- return self.__class__(self._parent, x, dx=dx)
-
- def _div_(self, other):
- p = self._parent.prime()
- x_self = self.approximation()
- x_other = other.approximation()
- wp_other = other.working_precision()
- d, inv, _ = x_other.xgcd(p**wp_other)
- q, r = x_self.quo_rem(d)
- if r != 0:
- raise ValueError("division is not exact")
- x = q * inv
- # dx = (1/other)*dself - (self/other^2)*dother
- dx = [ [weakref.ref(self), inv/d ],
- [weakref.ref(other), -x*inv/d ] ]
- return self.__class__(self._parent, x, dx=dx)
+ 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 e929cbb..568ce4d 100644
--- a/src/sage/rings/padics/padic_base_generic.py
+++ b/src/sage/rings/padics/padic_base_generic.py
@@ -21,6 +21,7 @@ from __future__ import absolute_import
#*****************************************************************************
from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
from .padic_generic import pAdicGeneric
from sage.rings.padics.pow_computer import PowComputer
@@ -45,9 +46,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)]
diff --git a/src/sage/rings/padics/padic_base_leaves.py b/src/sage/rings/padics/padic_base_leaves.py
index 0bb973d..3d95fba 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"""
@@ -810,3 +817,185 @@ class pAdicFieldFloatingPoint(pAdicFieldBaseGeneric, pAdicFloatingPointFieldGene
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)
+ 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..007d2fe
--- /dev/null
+++ b/src/sage/rings/padics/padic_lattice_element.py
@@ -0,0 +1,201 @@
+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):
+ 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.reduce(prec)
+ else:
+ self._value = pRational(parent.prime(), QQ(x)).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):