Lattice precision for p-adics

xcaruso commented 7 years ago

In several recent papers, David Roe, Tristan Vaccon and I explain that lattices allow a sharp track of precision: if f is a function we want to evaluate and x is an input given with some uncertainty modeled by a lattice H, then the uncertainty on the output f(x) is exactly df_x(H).

For much more details, I refer to my lecture notes http://xavier.toonywood.org/papers/publis/course-padic.pdf

The aim of this ticket is to propose a rough implementation of these ideas.

You can play with the latest version of this by clicking on launch binder here.

Below is a small demo (extracted from the doctest).


   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 =  rac {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
   "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
   "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 "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 "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 "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) 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 "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}

CC: @roed314 @sagetrac-TristanVaccon @saraedum @mezzarobba @sagetrac-swewers

Component: padics

Keywords: sd87

Author: Xavier Caruso

Branch: f29502e

Reviewer: David Roe, Julian Rüth

Issue created by migration from https://trac.sagemath.org/ticket/23505

xcaruso commented 7 years ago

Branch: u/caruso/lattice_precision

xcaruso commented 7 years ago

New commits:

b46b474 Pseudo-code for lattice precision

xcaruso commented 7 years ago

Description changed:

--- 
+++ 
@@ -2,126 +2,141 @@

 For much more details, I refer to my lecture notes http://xavier.toonywood.org/papers/publis/course-padic.pdf

-The aim of this ticket is to propose to implement these ideas in [SageMath](../wiki/SageMath). 
-More precisely, I propose to create a new parent `QpLP` (LP for "lattice precision") for which the precision is tracked using lattices. This leads to some difficulties:
-- Since a lattice does not in general split properly, the precision datum is a global object which must be handled by the parent (and not by its elements). For this reason, the parent has to manage the precision itself for all its elements. In particular, he has to maintain an up-to-date list of all its elements (and should be then informed whenever one of its elements is collected by the garbage collector).
-- If f is not surjective, df_x(H) is no longer a lattice. However, dealing with general Zp-submodules (and not only lattices) is more complicated for several reasons (the first one is that they are not exact objects whereas lattices are). For this reason, we introduce a `working_precision_cap` and feel free to add to our precision lattice any multiple of `p^working_precision_cap` if needed.
+The aim of this ticket is to propose a rough implementation of these ideas. 

-## About implementation
-
-The pseudo-code below gives some hints about the implementation I have in mind. Comments are welcome! 
+Below is a small demo (extracted from the doctest).

-# The parent -############

sage: R = ZpLP(3)
sage: x = R(1,10)

-class pAdicFieldLattice(pAdicRingBaseGeneric):

Internal variables:
. self._working_precision_cap
a cap for the working precision
meaning that the precision lattice always contains p^(self._working_precision_cap)
. self._elements
list of weak references of elements in this parent
. self._precision_lattice
a matrix over ZZ representing the lattice of precision
(its columns are indexed by self._elements)

+Of course, when we multiply by 3, we gain one digit of absolute +precision::
def init(self, p, working_precision_cap, print_mode):
Initialize variables here
sage: 3*x
3 + O(3^11)
def _echelonize(self):
Echelonize the matrix giving the precision lattice

+The lattice precision machinery sees this even if we decompose +the computation into several steps::
def _add_element(self, elt, prec):
A new element in the parent has just been created
We should:
. add a weak reference to it to the list self._elements
. add a column to self._precision_lattice
and update this matrix according to prec
(and possibly the working precision cap)
NOTE: prec be either an integer or a formal linear combinaison
of the precision on the other elements of this parent
sage: y = x+x
sage: y
2 + O(3^10)
sage: x+y
3 + O(3^11)
def _del_element(self, elt):
The element elt has just been garbage collected
We should:
. remove it to the list self._elements
. remove the corresponding column of self._precision_lattice

+The same works for the multiplication::
def precision_absolute(self, elt):
Return the (optimal) absolute precision of the element elt
This precision can be read off on self._precision_lattice:
it is the smallest valuation of an entry of the column of
self._precision_lattice corresponding to the element elt
sage: z = x^2
sage: z
1 + O(3^10)
sage: x*z
1 + O(3^11)
def working_precision(self, elt):
Return the working precision of the element elt
This precision can be read off on the precision lattice:
it is the smallest integer n for which the precision
lattice contains p^n*[elt] where [elt] denotes the
basis vector corresponding to elt

+This comes more funny when we are working with elements given +at different precisions::
def _elementconstructor(self, x, prec):
We ask for the creation of an element in this parent
We should:
. create this element (called elt hereafter) from x
Note: elt is an instance of the class pAdicLatticeElement below
. call the method _new_element(elt, prec)
. install a callback so that when elt will be collected
by the garbage collector, the method _del_element will
be called
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)

-QpLP = pAdicFieldLattice

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::

-# The elements -############## +..MATH::

-class pAdicLatticeElement(Element):

Internal variable:
. self._approximation
an approximation of this p-adic number
it is defined modulo p^(working_precision)
un = \frac {u{n-1} u{n-3} + u{n-2}^2} {u_{n-4}}
def working_precision(self):
return self.parent().working_precision(self) +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::
def approximation(self):
We should:
. reduce self._approximation modulo p^(self.working_precision())
(and update self._approximation accordingly)
. return the result
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,(bd+cc)/a; print d
4 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
13 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
55 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
21975 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
6639 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
7186 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
569 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
253 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
4149 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
2899 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
3072 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
349 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
619 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
243 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
3 + O(2^2)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
2 + O(2^2)
def precision_absolute(self):
return self.parent().precision_absolute(self) +If instead, we use the lattice precision, everything goes well::
def valuation(self, secure=False):
We should:
. compute the valuation of self.approximation()
. compare it to the self.precision_absolute()
. if the former is less than the latter, we return the former
. otherwise, if secure is False, we return the latter
if secure is True, we raise an error
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,(bd+cc)/a; print d
4 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
13 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
55 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
21975 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
23023 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
31762 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
16953 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
16637 + O(2^15)
def precision_relative(self, secure=False):
return self.precision_absolute() - self.valuation(secure=secure)
def repr(self):
Print like this:
self.approximation() + O(p^self.precision_absolute())
def add(self, other):
We should:
. compute app = self.approximation() + other.approximation()
. create a new element from app and the precision given by the differential:
dapp = dself + dother
def mul(self, other):
We should:
. compute app = self.approximation() * other.approximation()
. create a new element from app and the precision given by the differential:
dapp = selfdother + otherdself
sage: for _ in range(100):
....: a,b,c,d = b,c,d,(bd+cc)/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)
+PS: For now, the code is far from being ready for review. Any help will be of course quite appreciated.

xcaruso commented 7 years ago

Commit: b46b474

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

4b2af06 First more-or-less working implementation

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Changed commit from b46b474 to 4b2af06

xcaruso commented 7 years ago

comment:5

I've worked more on my implementation and, at least, it seems now to be usable.

It still requires a lot of work (convert to Cython, write templates, write doctests, rewrite completely the class pRational which is a hack, etc). I however post it because I think that some people might have fun playing with it and discovering what lattice precision can do. Enjoy :-)

@mezzarobba: I add your name in Cc because I think that you could be interested. If you're not, feel free to remove it.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Changed commit from 4b2af06 to 858492b

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

858492b Second rough implementation of lattice precision

xcaruso commented 7 years ago

Description changed:

--- 
+++ 
@@ -8,135 +8,297 @@

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::
un = \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,(bd+cc)/a; print d
4 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
13 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
55 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
21975 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
6639 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
7186 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
569 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
253 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
4149 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
2899 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
3072 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
349 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
619 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
243 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
3 + O(2^2)
sage: a,b,c,d = b,c,d,(bd+cc)/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,(bd+cc)/a; print d
4 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
13 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
55 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
21975 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
23023 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
31762 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
16953 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
16637 + O(2^15)
sage: for _ in range(100):
....: a,b,c,d = b,c,d,(bd+cc)/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)
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::
un = rac {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,(bd+cc)/a; print d
4 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
13 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
55 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
21975 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
6639 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
7186 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
569 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
253 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
4149 + O(2^13)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
2899 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
3072 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
349 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
619 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
243 + O(2^12)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
3 + O(2^2)
sage: a,b,c,d = b,c,d,(bd+cc)/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,(bd+cc)/a; print d
4 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
13 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
55 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
21975 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
23023 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
31762 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
16953 + O(2^15)
sage: a,b,c,d = b,c,d,(bd+cc)/a; print d
16637 + O(2^15)
sage: for _ in range(100):
....: a,b,c,d = b,c,d,(bd+cc)/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
"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
"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 "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 "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 "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 = 3x+2y, 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 "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 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
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}

-PS: For now, the code is far from being ready for review. Any help will be of course quite appreciated.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`418cff6`	`Typos`
`399e953`	`Fix element_class`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Changed commit from 858492b to 399e953

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Changed commit from 399e953 to 1951205

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

1951205 Fix convert_multiple

roed314 commented 7 years ago

Changed branch from u/caruso/lattice_precision to u/roed/lattice_precision

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Changed commit from 1951205 to 90a79f7

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

90a79f7 Merge branch 'u/roed/lattice_precision' of git://trac.sagemath.org/sage into t/23505/lattice_precision

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

aa8ae62 Fix small problem in QpLP in factory

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 7 years ago

Changed commit from 90a79f7 to aa8ae62

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from aa8ae62 to 0bf65f1

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

0bf65f1 Merge branch 'u/roed/lattice_precision' of git://trac.sagemath.org/sage into t/23505/lattice_precision

xcaruso commented 6 years ago

Changed branch from u/roed/lattice_precision to u/caruso/lattice_precision

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 0bf65f1 to 62c5d56

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

62c5d56 Some doctests in lattice_precision.py

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 62c5d56 to f06603b

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`b70227b`	`More doctests in lattice_precision.py`
`f06603b`	`Doctest in padic_base_leaves.py`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from f06603b to 5ef1735

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

5ef1735 Doctest in padic_lattice_element.py

roed314 commented 6 years ago

Changed branch from u/caruso/lattice_precision to u/roed/lattice_precision

xcaruso commented 6 years ago

Changed branch from u/roed/lattice_precision to u/caruso/lattice_precision

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 5ef1735 to e9bb90b

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

e9bb90b ZpLP -> ZpLC

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

087eb33 Fix small bug

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from e9bb90b to 087eb33

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 087eb33 to 8d68a69

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

8d68a69 Implementation of ZpLF

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`b0a8b0b`	`Doctest for ZpLF`
`7bb677c`	`Add a thresold for column deletion`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 8d68a69 to 7bb677c

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

63eee81 Merge branch 'develop' into lattice_precision

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 7bb677c to 63eee81

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 63eee81 to c04c370

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`009bbfd`	`Merge branch 'develop' into lattice_precision`
`c04c370`	`Doctest for abstract methods`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from c04c370 to c26c01e

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

c26c01e Fix bug in is_precision_capped

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from c26c01e to 9780da5

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`c3c8692`	`Doctest for the class pRational`
`9780da5`	`More doctests`

xcaruso commented 6 years ago

comment:29

I guess that this ticket is now more or less ready for a first review.

saraedum commented 6 years ago

comment:30

I would be happy to help with reviewing this. Since this is quite a lot of code, would you mind to push the code to github/gitlab where it is much easier to comment on lines of code? I can also push it to a repository there myself if you don't mind.

saraedum commented 6 years ago

comment:31

Do you document somewhere why this does not follow the usual linkage pattern used elsewhere in the p-adics code?

roed314 commented 6 years ago

comment:32

Replying to @saraedum:

Do you document somewhere why this does not follow the usual linkage pattern used elsewhere in the p-adics code?

I think the main answer is that it just hasn't been switched to the linkage pattern, but it should be. The pRational will take some work to make this happen. I'll be busy in the next few weeks working on an ANTS paper, but we can talk about this in Rennes. In the mean time, I think if you want to push the code to github/gitlab that should be fine.

sagemath / sage

Lattice precision for p-adics #23505

Internal variables:

. self._working_precision_cap

a cap for the working precision

meaning that the precision lattice always contains p^(self._working_precision_cap)

. self._elements

list of weak references of elements in this parent

. self._precision_lattice

a matrix over ZZ representing the lattice of precision

(its columns are indexed by self._elements)

Initialize variables here

Echelonize the matrix giving the precision lattice

A new element in the parent has just been created

We should:

. add a weak reference to it to the list self._elements

. add a column to self._precision_lattice

and update this matrix according to prec

(and possibly the working precision cap)

NOTE: prec be either an integer or a formal linear combinaison

of the precision on the other elements of this parent

The element elt has just been garbage collected

We should:

. remove it to the list self._elements

. remove the corresponding column of self._precision_lattice

Return the (optimal) absolute precision of the element elt

This precision can be read off on self._precision_lattice:

it is the smallest valuation of an entry of the column of

self._precision_lattice corresponding to the element elt

Return the working precision of the element elt

This precision can be read off on the precision lattice:

it is the smallest integer n for which the precision

lattice contains p^n*[elt] where [elt] denotes the

basis vector corresponding to elt

We ask for the creation of an element in this parent

We should:

. create this element (called elt hereafter) from x

Note: elt is an instance of the class pAdicLatticeElement below

. call the method _new_element(elt, prec)

. install a callback so that when elt will be collected

by the garbage collector, the method _del_element will

be called

Internal variable:

. self._approximation

an approximation of this p-adic number

it is defined modulo p^(working_precision)

We should:

. reduce self._approximation modulo p^(self.working_precision())

(and update self._approximation accordingly)

. return the result

We should:

. compute the valuation of self.approximation()

. compare it to the self.precision_absolute()

. if the former is less than the latter, we return the former

. otherwise, if secure is False, we return the latter

if secure is True, we raise an error

Print like this:

self.approximation() + O(p^self.precision_absolute())

We should:

. compute app = self.approximation() + other.approximation()

. create a new element from app and the precision given by the differential:

dapp = dself + dother

We should:

. compute app = self.approximation() * other.approximation()

. create a new element from app and the precision given by the differential:

dapp = selfdother + otherdself