Changes between Version 2 and Version 7 of Ticket #23505


Ignore:
Timestamp:
08/02/17 00:00:07 (16 months ago)
Author:
caruso
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • Ticket #23505

    • Property Cc saraedum mmezzarobba added
    • Property Commit changed from b46b474adaa8d0f54a2bfa039939e9185ee79b4e to 858492b4b0cd87dfb63fba0f3ad239801b98cb62
    • Property Type changed from task to enhancement
  • Ticket #23505 – Description

    v2 v7  
    99{{{
    1010
    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}
    139302
    140303}}}
    141 
    142 PS: For now, the code is far from being ready for review. Any help will be of course quite appreciated.