# Changes between Initial Version and Version 2 of Ticket #23505

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

New commits:

 ​b46b474 Pseudo-code for lattice precision

### Legend:

Unmodified
• Property Commit changed from to b46b474adaa8d0f54a2bfa039939e9185ee79b4e
• Property Branch changed from to u/caruso/lattice_precision
 initial 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. 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 _element_constructor_(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) u_n = \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,(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) 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,(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) 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 = self*dother + other*dself 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) }}} PS: For now, the code is far from being ready for review. Any help will be of course quite appreciated.