Changes between Initial Version and Version 2 of Ticket #23505
 Timestamp:
 07/22/17 04:09:02 (9 months ago)
Legend:
 Unmodified
 Added
 Removed
 Modified

Ticket #23505

Property
Commit
changed from
to
b46b474adaa8d0f54a2bfa039939e9185ee79b4e

Property
Branch
changed from
to
u/caruso/lattice_precision

Property
Commit
changed from

Ticket #23505 – Description
initial v2 3 3 For much more details, I refer to my lecture notes http://xavier.toonywood.org/papers/publis/coursepadic.pdf 4 4 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 uptodate 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 Zpsubmodules (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. 5 The aim of this ticket is to propose a rough implementation of these ideas. 9 6 10 == About implementation 11 12 The pseudocode below gives some hints about the implementation I have in mind. Comments are welcome! 7 Below is a small demo (extracted from the doctest). 13 8 14 9 {{{ 15 10 16 # The parent 17 ############ 11 sage: R = ZpLP(3) 12 sage: x = R(1,10) 18 13 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) 14 Of course, when we multiply by 3, we gain one digit of absolute 15 precision:: 29 16 30 def __init__(self, p, working_precision_cap, print_mode):31 # Initialize variables here17 sage: 3*x 18 3 + O(3^11) 32 19 33 def _echelonize(self): 34 # Echelonize the matrix giving the precision lattice 20 The lattice precision machinery sees this even if we decompose 21 the computation into several steps:: 35 22 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) 45 28 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 29 The same works for the multiplication:: 51 30 52 def precision_absolute(self, elt):53 # Return the (optimal) absolute precision of the element elt54 # This precision can be read off on self._precision_lattice:55 # it is the smallest valuation of an entry of the column of56 # self._precision_lattice corresponding to the element elt31 sage: z = x^2 32 sage: z 33 1 + O(3^10) 34 sage: x*z 35 1 + O(3^11) 57 36 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 37 This comes more funny when we are working with elements given 38 at different precisions:: 64 39 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 = xy; t 46 0 + O(2^5) 47 sage: z+t # observe that z+t = 2*x 48 2 + O(2^11) 49 sage: zt # observe that zt = 2*y 50 2 + O(2^6) 74 51 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) 76 58 59 The SOMOS sequence is the sequence defined by the recurrence:: 77 60 78 # The elements 79 ############## 61 ..MATH:: 80 62 81 class pAdicLatticeElement(Element): 82 # Internal variable: 83 # . self._approximation 84 # an approximation of this padic number 85 # it is defined modulo p^(working_precision) 63 u_n = \frac {u_{n1} u_{n3} + u_{n2}^2} {u_{n4}} 86 64 87 def working_precision(self): 88 return self.parent().working_precision(self) 65 It is known for its numerical instability. 66 On the one hand, one can show that if the initial values are 67 invertible in `\mathbb{Z}_p` and known at precision `O(p^N)` 68 then all the next terms of the SOMOS sequence will be known 69 at the same precision as well. 70 On the other hand, because of the division, when we unroll 71 the recurrence, we loose a lot of precision. Observe:: 89 72 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) 95 107 96 def precision_absolute(self): 97 return self.parent().precision_absolute(self) 108 If instead, we use the lattice precision, everything goes well:: 98 109 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) 106 128 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) 126 139 127 140 }}} 141 142 PS: For now, the code is far from being ready for review. Any help will be of course quite appreciated.