11 | | sage: R = ZpLP(3) |
12 | | sage: x = R(1,10) |
13 | | |
14 | | Of course, when we multiply by 3, we gain one digit of absolute |
15 | | precision:: |
16 | | |
17 | | sage: 3*x |
18 | | 3 + O(3^11) |
19 | | |
20 | | The lattice precision machinery sees this even if we decompose |
21 | | the computation into several steps:: |
22 | | |
23 | | sage: y = x+x |
24 | | sage: y |
25 | | 2 + O(3^10) |
26 | | sage: x+y |
27 | | 3 + O(3^11) |
28 | | |
29 | | The same works for the multiplication:: |
30 | | |
31 | | sage: z = x^2 |
32 | | sage: z |
33 | | 1 + O(3^10) |
34 | | sage: x*z |
35 | | 1 + O(3^11) |
36 | | |
37 | | This comes more funny when we are working with elements given |
38 | | at different precisions:: |
39 | | |
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) |
51 | | |
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) |
58 | | |
59 | | The SOMOS sequence is the sequence defined by the recurrence:: |
60 | | |
61 | | ..MATH:: |
62 | | |
63 | | u_n = \frac {u_{n-1} u_{n-3} + u_{n-2}^2} {u_{n-4}} |
64 | | |
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:: |
72 | | |
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) |
107 | | |
108 | | If instead, we use the lattice precision, everything goes well:: |
109 | | |
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) |
128 | | |
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) |
| 11 | Below is a small demo of the features by this model of precision: |
| 12 | |
| 13 | sage: R = ZpLP(3, print_mode='terse') |
| 14 | sage: x = R(1,10) |
| 15 | |
| 16 | Of course, when we multiply by 3, we gain one digit of absolute |
| 17 | precision: |
| 18 | |
| 19 | sage: 3*x |
| 20 | 3 + O(3^11) |
| 21 | |
| 22 | The lattice precision machinery sees this even if we decompose the |
| 23 | computation into several steps: |
| 24 | |
| 25 | sage: y = x+x |
| 26 | sage: y |
| 27 | 2 + O(3^10) |
| 28 | sage: x + y |
| 29 | 3 + O(3^11) |
| 30 | |
| 31 | The same works for the multiplication: |
| 32 | |
| 33 | sage: z = x^2 |
| 34 | sage: z |
| 35 | 1 + O(3^10) |
| 36 | sage: x*z |
| 37 | 1 + O(3^11) |
| 38 | |
| 39 | This comes more funny when we are working with elements given at |
| 40 | different precisions: |
| 41 | |
| 42 | sage: R = ZpLP(2, print_mode='terse') |
| 43 | sage: x = R(1,10) |
| 44 | sage: y = R(1,5) |
| 45 | sage: z = x+y; z |
| 46 | 2 + O(2^5) |
| 47 | sage: t = x-y; t |
| 48 | 0 + O(2^5) |
| 49 | sage: z+t # observe that z+t = 2*x |
| 50 | 2 + O(2^11) |
| 51 | sage: z-t # observe that z-t = 2*y |
| 52 | 2 + O(2^6) |
| 53 | |
| 54 | sage: x = R(28888,15) |
| 55 | sage: y = R(204,10) |
| 56 | sage: z = x/y; z |
| 57 | 242 + O(2^9) |
| 58 | sage: z*y # which is x |
| 59 | 28888 + O(2^15) |
| 60 | |
| 61 | The SOMOS sequence is the sequence defined by the recurrence: |
| 62 | |
| 63 | ..MATH:: |
| 64 | |
| 65 | u_n = rac {u_{n-1} u_{n-3} + u_{n-2}^2} {u_{n-4}} |
| 66 | |
| 67 | It is known for its numerical instability. On the one hand, one can |
| 68 | show that if the initial values are invertible in mathbb{Z}_p and |
| 69 | known at precision O(p^N) then all the next terms of the SOMOS |
| 70 | sequence will be known at the same precision as well. On the other |
| 71 | hand, because of the division, when we unroll the recurrence, we |
| 72 | loose a lot of precision. Observe: |
| 73 | |
| 74 | sage: R = Zp(2, 30, print_mode='terse') |
| 75 | sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15) |
| 76 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 77 | 4 + O(2^15) |
| 78 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 79 | 13 + O(2^15) |
| 80 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 81 | 55 + O(2^15) |
| 82 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 83 | 21975 + O(2^15) |
| 84 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 85 | 6639 + O(2^13) |
| 86 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 87 | 7186 + O(2^13) |
| 88 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 89 | 569 + O(2^13) |
| 90 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 91 | 253 + O(2^13) |
| 92 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 93 | 4149 + O(2^13) |
| 94 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 95 | 2899 + O(2^12) |
| 96 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 97 | 3072 + O(2^12) |
| 98 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 99 | 349 + O(2^12) |
| 100 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 101 | 619 + O(2^12) |
| 102 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 103 | 243 + O(2^12) |
| 104 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 105 | 3 + O(2^2) |
| 106 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 107 | 2 + O(2^2) |
| 108 | |
| 109 | If instead, we use the lattice precision, everything goes well: |
| 110 | |
| 111 | sage: R = ZpLP(2, 30, print_mode='terse') |
| 112 | sage: a,b,c,d = R(1,15), R(1,15), R(1,15), R(3,15) |
| 113 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 114 | 4 + O(2^15) |
| 115 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 116 | 13 + O(2^15) |
| 117 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 118 | 55 + O(2^15) |
| 119 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 120 | 21975 + O(2^15) |
| 121 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 122 | 23023 + O(2^15) |
| 123 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 124 | 31762 + O(2^15) |
| 125 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 126 | 16953 + O(2^15) |
| 127 | sage: a,b,c,d = b,c,d,(b*d+c*c)/a; print d |
| 128 | 16637 + O(2^15) |
| 129 | |
| 130 | sage: for _ in range(100): |
| 131 | ....: a,b,c,d = b,c,d,(b*d+c*c)/a |
| 132 | sage: a |
| 133 | 15519 + O(2^15) |
| 134 | sage: b |
| 135 | 32042 + O(2^15) |
| 136 | sage: c |
| 137 | 17769 + O(2^15) |
| 138 | sage: d |
| 139 | 20949 + O(2^15) |
| 140 | |
| 141 | BEHIND THE SCENE: |
| 142 | |
| 143 | The precision is global. It is encoded by a lattice in a huge |
| 144 | vector space whose dimension is the number of elements having this |
| 145 | parent. |
| 146 | |
| 147 | Concretely, this precision datum is an instance of the class |
| 148 | "sage.rings.padic.lattice_precision.PrecisionLattice". It is |
| 149 | attached to the parent and is created at the same time as the |
| 150 | parent. (It is actually a bit more subtle because two different |
| 151 | parents may share the same instance; this happens for instance for |
| 152 | a p-adic ring and its field of fractions.) |
| 153 | |
| 154 | This precision datum is accessible through the method |
| 155 | "precision()": |
| 156 | |
| 157 | sage: R = ZpLP(5, print_mode='terse') |
| 158 | sage: prec = R.precision() |
| 159 | sage: prec |
| 160 | Precision Lattice on 0 object |
| 161 | |
| 162 | This instance knows about all elements of the parent, it is |
| 163 | automatically updated when a new element (of this parent) is |
| 164 | created: |
| 165 | |
| 166 | sage: x = R(3513,10) |
| 167 | sage: prec |
| 168 | Precision Lattice on 1 object |
| 169 | sage: y = R(176,5) |
| 170 | sage: prec |
| 171 | Precision Lattice on 2 objects |
| 172 | sage: z = R.random_element() |
| 173 | sage: prec |
| 174 | Precision Lattice on 3 objects |
| 175 | |
| 176 | The method "tracked_elements()" provides the list of all tracked |
| 177 | elements: |
| 178 | |
| 179 | sage: prec.tracked_elements() |
| 180 | [3513 + O(5^10), 176 + O(5^5), ...] |
| 181 | |
| 182 | Similarly, when a variable is collected by the garbage collector, |
| 183 | the precision lattice is updated. Note however that the update |
| 184 | might be delayed. We can force it with the method "del_elements()": |
| 185 | |
| 186 | sage: z = 0 |
| 187 | sage: prec |
| 188 | Precision Lattice on 3 objects |
| 189 | sage: prec.del_elements() |
| 190 | sage: prec |
| 191 | Precision Lattice on 2 objects |
| 192 | |
| 193 | The method "precision_lattice()" returns (a matrix defining) the |
| 194 | lattice that models the precision. Here we have: |
| 195 | |
| 196 | sage: prec.precision_lattice() |
| 197 | [9765625 0] |
| 198 | [ 0 3125] |
| 199 | |
| 200 | Observe that 5^10 = 9765625 and 5^5 = 3125. The above matrix then |
| 201 | reflects the precision on x and y. |
| 202 | |
| 203 | Now, observe how the precision lattice changes while performing |
| 204 | computations: |
| 205 | |
| 206 | sage: x, y = 3*x+2*y, 2*(x-y) |
| 207 | sage: prec.del_elements() |
| 208 | sage: prec.precision_lattice() |
| 209 | [ 3125 48825000] |
| 210 | [ 0 48828125] |
| 211 | |
| 212 | The matrix we get is no longer diagonal, meaning that some digits |
| 213 | of precision are diffused among the two new elements x and y. They |
| 214 | nevertheless show up when we compute for instance x+y: |
| 215 | |
| 216 | sage: x |
| 217 | 1516 + O(5^5) |
| 218 | sage: y |
| 219 | 424 + O(5^5) |
| 220 | sage: x+y |
| 221 | 17565 + O(5^11) |
| 222 | |
| 223 | It is these diffused digits of precision (which are tracked but do |
| 224 | not appear on the printing) that allow to be always sharp on |
| 225 | precision. |
| 226 | |
| 227 | PERFORMANCES: |
| 228 | |
| 229 | Each elementary operation requires significant manipulations on the |
| 230 | lattice precision and then is costly. Precisely: |
| 231 | |
| 232 | * The creation of a new element has a cost O(n) when n is the |
| 233 | number of tracked elements. |
| 234 | |
| 235 | * The destruction of one element has a cost O(m^2) when m is the |
| 236 | distance between the destroyed element and the last one. |
| 237 | Fortunately, it seems that m tends to be small in general (the |
| 238 | dynamics of the list of tracked elements is rather close to that |
| 239 | of a stack). |
| 240 | |
| 241 | It is nevertheless still possible to manipulate several hundred |
| 242 | variables (e.g. squares matrices of size 5 or polynomials of degree |
| 243 | 20 are accessible). |
| 244 | |
| 245 | The class "PrecisionLattice" provides several features for |
| 246 | introspection (especially concerning timings). If enables, it |
| 247 | maintains an history of all actions and stores the wall time of |
| 248 | each of them: |
| 249 | |
| 250 | sage: R = ZpLP(3) |
| 251 | sage: prec = R.precision() |
| 252 | sage: prec.history_enable() |
| 253 | sage: M = random_matrix(R, 5) |
| 254 | sage: d = M.determinant() |
| 255 | sage: print prec.history() # somewhat random |
| 256 | --- |
| 257 | 0.004212s oooooooooooooooooooooooooooooooooooo |
| 258 | 0.000003s oooooooooooooooooooooooooooooooooo~~ |
| 259 | 0.000010s oooooooooooooooooooooooooooooooooo |
| 260 | 0.001560s ooooooooooooooooooooooooooooooooooooooooo |
| 261 | 0.000004s ooooooooooooooooooooooooooooo~oooo~oooo~o |
| 262 | 0.002168s oooooooooooooooooooooooooooooooooooooo |
| 263 | 0.001787s ooooooooooooooooooooooooooooooooooooooooo |
| 264 | 0.000004s oooooooooooooooooooooooooooooooooooooo~~o |
| 265 | 0.000198s ooooooooooooooooooooooooooooooooooooooo |
| 266 | 0.001152s ooooooooooooooooooooooooooooooooooooooooo |
| 267 | 0.000005s ooooooooooooooooooooooooooooooooo~oooo~~o |
| 268 | 0.000853s oooooooooooooooooooooooooooooooooooooo |
| 269 | 0.000610s ooooooooooooooooooooooooooooooooooooooo |
| 270 | ... |
| 271 | 0.003879s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 272 | 0.000006s oooooooooooooooooooooooooooooooooooooooooooooooooooo~~~~~ |
| 273 | 0.000036s oooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 274 | 0.006737s oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 275 | 0.000005s oooooooooooooooooooooooooooooooooooooooooooooooooooo~~~~~ooooo |
| 276 | 0.002637s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 277 | 0.007118s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 278 | 0.000008s oooooooooooooooooooooooooooooooooooooooooooooooooooo~~~~o~~~~oooo |
| 279 | 0.003504s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 280 | 0.005371s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 281 | 0.000006s ooooooooooooooooooooooooooooooooooooooooooooooooooooo~~~o~~~ooo |
| 282 | 0.001858s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 283 | 0.003584s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 284 | 0.000004s oooooooooooooooooooooooooooooooooooooooooooooooooooooo~~o~~oo |
| 285 | 0.000801s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 286 | 0.001916s ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo |
| 287 | 0.000022s ooooooooooooooooooooooooooooo~~~~~~~~~~~~~~~~~~~~~~oooo~o~o |
| 288 | 0.014705s ooooooooooooooooooooooooooooooooooo |
| 289 | 0.001292s ooooooooooooooooooooooooooooooooooooo |
| 290 | 0.000002s ooooooooooooooooooooooooooooooooooo~o |
| 291 | |
| 292 | The symbol o symbolized a tracked element. The symbol ~ means that |
| 293 | the element is marked for deletion. |
| 294 | |
| 295 | The global timings are also accessible as follows: |
| 296 | |
| 297 | sage: prec.timings() # somewhat random |
| 298 | {'add': 0.25049376487731934, |
| 299 | 'del': 0.11911273002624512, |
| 300 | 'mark': 0.0004909038543701172, |
| 301 | 'partial reduce': 0.0917658805847168} |