Opened 19 months ago

# Lattice precision for p-adics — at Version 2

Reported by: Owned by: caruso major sage-8.1 padics sd87 roed, TristanVaccon, saraedum, mmezzarobba, swewers N/A u/caruso/lattice_precision (Commits) b46b474adaa8d0f54a2bfa039939e9185ee79b4e

In several recent papers, David Roe, Tristan Vaccon and I explain that lattices allow a sharp track of precision: if f is a function we want to evaluate and x is an input given with some uncertainty modeled by a lattice H, then the uncertainty on the output f(x) is exactly df_x(H).

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 a rough implementation of these ideas.

Below is a small demo (extracted from the doctest).

    sage: R = ZpLP(3)
sage: x = R(1,10)

Of course, when we multiply by 3, we gain one digit of absolute
precision::

sage: 3*x
3 + O(3^11)

The lattice precision machinery sees this even if we decompose
the computation into several steps::

sage: y = x+x
sage: y
2 + O(3^10)
sage: x+y
3 + O(3^11)

The same works for the multiplication::

sage: z = x^2
sage: z
1 + O(3^10)
sage: x*z
1 + O(3^11)

This comes more funny when we are working with elements given
at different precisions::

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)

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::

..MATH::

u_n = \frac {u_{n-1} u_{n-3} + u_{n-2}^2} {u_{n-4}}

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::

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)

If instead, we use the lattice precision, everything goes well::

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)

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.

### comment:1 Changed 19 months ago by caruso

• Branch set to u/caruso/lattice_precision

### comment:2 Changed 19 months ago by caruso

 ​b46b474 Pseudo-code for lattice precision