Opened 19 months ago

Last modified 4 months ago

#23505 closed enhancement

Lattice precision for p-adics — at Version 2

Reported by: caruso Owned by:
Priority: major Milestone: sage-8.1
Component: padics Keywords: sd87
Cc: roed, TristanVaccon, saraedum, mmezzarobba, swewers Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: u/caruso/lattice_precision (Commits) Commit: b46b474adaa8d0f54a2bfa039939e9185ee79b4e
Dependencies: Stopgaps:

Description (last modified by caruso)

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.

Change History (2)

comment:1 Changed 19 months ago by caruso

  • Branch set to u/caruso/lattice_precision

comment:2 Changed 19 months ago by caruso

  • Commit set to b46b474adaa8d0f54a2bfa039939e9185ee79b4e
  • Description modified (diff)

New commits:

b46b474Pseudo-code for lattice precision
Note: See TracTickets for help on using tickets.