Changes between Initial Version and Version 2 of Ticket #23505


Ignore:
Timestamp:
07/22/17 04:09:02 (16 months ago)
Author:
caruso
Comment:

New commits:

b46b474Pseudo-code for lattice precision

Legend:

Unmodified
Added
Removed
Modified
  • Ticket #23505

    • Property Commit changed from to b46b474adaa8d0f54a2bfa039939e9185ee79b4e
    • Property Branch changed from to u/caruso/lattice_precision
  • Ticket #23505 – Description

    initial v2  
    33For much more details, I refer to my lecture notes http://xavier.toonywood.org/papers/publis/course-padic.pdf
    44
    5 The aim of this ticket is to propose to implement these ideas in SageMath.
    6 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:
    7 - 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).
    8 - 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.
     5The aim of this ticket is to propose a rough implementation of these ideas.
    96
    10 == About implementation
    11 
    12 The pseudo-code below gives some hints about the implementation I have in mind. Comments are welcome!
     7Below is a small demo (extracted from the doctest).
    138
    149{{{
    1510
    16 # The parent
    17 ############
     11    sage: R = ZpLP(3)
     12    sage: x = R(1,10)
    1813
    19 class pAdicFieldLattice(pAdicRingBaseGeneric):
    20     # Internal variables:
    21     #  . self._working_precision_cap
    22     #    a cap for the working precision
    23     #    meaning that the precision lattice always contains p^(self._working_precision_cap)
    24     #  . self._elements
    25     #    list of weak references of elements in this parent
    26     #  . self._precision_lattice
    27     #    a matrix over ZZ representing the lattice of precision
    28     #    (its columns are indexed by self._elements)
     14Of course, when we multiply by 3, we gain one digit of absolute
     15precision::
    2916
    30     def __init__(self, p, working_precision_cap, print_mode):
    31         # Initialize variables here
     17    sage: 3*x
     18    3 + O(3^11)
    3219
    33     def _echelonize(self):
    34         # Echelonize the matrix giving the precision lattice
     20The lattice precision machinery sees this even if we decompose
     21the computation into several steps::
    3522
    36     def _add_element(self, elt, prec):
    37         # A new element in the parent has just been created
    38         # We should:
    39         #  . add a weak reference to it to the list self._elements
    40         #  . add a column to self._precision_lattice
    41         #    and update this matrix according to prec
    42         #    (and possibly the working precision cap)
    43         # NOTE: prec be either an integer or a formal linear combinaison
    44         # of the precision on the other elements of this parent
     23    sage: y = x+x
     24    sage: y
     25    2 + O(3^10)
     26    sage: x+y
     27    3 + O(3^11)
    4528
    46     def _del_element(self, elt):
    47         # The element elt has just been garbage collected
    48         # We should:
    49         #  . remove it to the list self._elements
    50         #  . remove the corresponding column of self._precision_lattice
     29The same works for the multiplication::
    5130
    52     def precision_absolute(self, elt):
    53         # Return the (optimal) absolute precision of the element elt
    54         # This precision can be read off on self._precision_lattice:
    55         # it is the smallest valuation of an entry of the column of
    56         # self._precision_lattice corresponding to the element elt
     31    sage: z = x^2
     32    sage: z
     33    1 + O(3^10)
     34    sage: x*z
     35    1 + O(3^11)
    5736
    58     def working_precision(self, elt):
    59         # Return the working precision of the element elt
    60         # This precision can be read off on the precision lattice:
    61         # it is the smallest integer n for which the precision
    62         # lattice contains p^n*[elt] where [elt] denotes the
    63         # basis vector corresponding to elt
     37This comes more funny when we are working with elements given
     38at different precisions::
    6439
    65     def _element_constructor_(self, x, prec):
    66         # We ask for the creation of an element in this parent
    67         # We should:
    68         #  . create this element (called elt hereafter) from x
    69         #    Note: elt is an instance of the class pAdicLatticeElement below
    70         #  . call the method _new_element(elt, prec)
    71         #  . install a callback so that when elt will be collected
    72         #    by the garbage collector, the method _del_element will
    73         #    be called
     40    sage: R = ZpLP(2)
     41    sage: x = R(1,10)
     42    sage: y = R(1,5)
     43    sage: z = x+y; z
     44    2 + O(2^5)
     45    sage: t = x-y; t
     46    0 + O(2^5)
     47    sage: z+t  # observe that z+t = 2*x
     48    2 + O(2^11)
     49    sage: z-t  # observe that z-t = 2*y
     50    2 + O(2^6)
    7451
    75 QpLP = pAdicFieldLattice
     52    sage: x = R(28888,15)
     53    sage: y = R(204,10)
     54    sage: z = x/y; z
     55    242 + O(2^9)
     56    sage: z*y  # which is x
     57    28888 + O(2^15)
    7658
     59The SOMOS sequence is the sequence defined by the recurrence::
    7760
    78 # The elements
    79 ##############
     61..MATH::
    8062
    81 class pAdicLatticeElement(Element):
    82     # Internal variable:
    83     #  . self._approximation
    84     #    an approximation of this p-adic number
    85     #    it is defined modulo p^(working_precision)
     63    u_n = \frac {u_{n-1} u_{n-3} + u_{n-2}^2} {u_{n-4}}
    8664
    87     def working_precision(self):
    88         return self.parent().working_precision(self)
     65It is known for its numerical instability.
     66On the one hand, one can show that if the initial values are
     67invertible in `\mathbb{Z}_p` and known at precision `O(p^N)`
     68then all the next terms of the SOMOS sequence will be known
     69at the same precision as well.
     70On the other hand, because of the division, when we unroll
     71the recurrence, we loose a lot of precision. Observe::
    8972
    90     def approximation(self):
    91         # We should:
    92         #  . reduce self._approximation modulo p^(self.working_precision())
    93         #    (and update self._approximation accordingly)
    94         #  . return the result
     73    sage: R = Zp(2, print_mode='terse')
     74    sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15)
     75    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     76    4 + O(2^15)
     77    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     78    13 + O(2^15)
     79    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     80    55 + O(2^15)
     81    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     82    21975 + O(2^15)
     83    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     84    6639 + O(2^13)
     85    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     86    7186 + O(2^13)
     87    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     88    569 + O(2^13)
     89    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     90    253 + O(2^13)
     91    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     92    4149 + O(2^13)
     93    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     94    2899 + O(2^12)
     95    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     96    3072 + O(2^12)
     97    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     98    349 + O(2^12)
     99    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     100    619 + O(2^12)
     101    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     102    243 + O(2^12)
     103    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     104    3 + O(2^2)
     105    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     106    2 + O(2^2)
    95107
    96     def precision_absolute(self):
    97         return self.parent().precision_absolute(self)
     108If instead, we use the lattice precision, everything goes well::
    98109
    99     def valuation(self, secure=False):
    100         # We should:
    101         #  . compute the valuation of self.approximation()
    102         #  . compare it to the self.precision_absolute()
    103         #  . if the former is less than the latter, we return the former
    104         #  . otherwise, if secure is False, we return the latter
    105         #               if secure is True, we raise an error
     110    sage: R = ZpLP(2)
     111    sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15)
     112    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     113    4 + O(2^15)
     114    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     115    13 + O(2^15)
     116    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     117    55 + O(2^15)
     118    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     119    21975 + O(2^15)
     120    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     121    23023 + O(2^15)
     122    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     123    31762 + O(2^15)
     124    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     125    16953 + O(2^15)
     126    sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d
     127    16637 + O(2^15)
    106128
    107     def precision_relative(self, secure=False):
    108         return self.precision_absolute() - self.valuation(secure=secure)
    109 
    110     def _repr_(self):
    111         # Print like this:
    112         #   self.approximation() + O(p^self.precision_absolute())
    113 
    114 
    115     def _add_(self, other):
    116         # We should:
    117         #  . compute app = self.approximation() + other.approximation()
    118         #  . create a new element from app and the precision given by the differential:
    119         #    dapp = dself + dother
    120 
    121     def _mul_(self, other):
    122         # We should:
    123         #  . compute app = self.approximation() * other.approximation()
    124         #  . create a new element from app and the precision given by the differential:
    125         #    dapp = self*dother + other*dself
     129    sage: for _ in range(100):
     130    ....:     a,b,c,d = b,c,d,(b*d+c*c)/a
     131    sage: a
     132    15519 + O(2^15)
     133    sage: b
     134    32042 + O(2^15)
     135    sage: c
     136    17769 + O(2^15)
     137    sage: d
     138    20949 + O(2^15)
    126139
    127140}}}
     141
     142PS: For now, the code is far from being ready for review. Any help will be of course quite appreciated.