search

sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.4k stars 473 forks source link

Relaxed p-adics #31108

Closed xcaruso closed 3 years ago

xcaruso commented 3 years ago

We propose an implementation of exact p-adic numbers relaying on relaxed arithmetics proposed by van der Hoeven and al.

Here is a small demo:

    The model for relaxed `p`-adics is quite different from any of the
    other types of `p`-adics. In addition to storing a finite
    approximation, one also stores a method for increasing the
    precision.

    Relaxed `p`-adic rings are created by the constructor :func:`ZpER`::

        sage: R = ZpER(5, print_mode="digits")
        sage: R
        5-adic Ring handled with relaxed arithmetics

    The precision is not capped in `R`::

        sage: R.precision_cap()
        +Infinity

    However, a default precision is fixed. This is the precision
    at which the elements will be printed::

        sage: R.default_prec()
        20

    A default halting precision is also set. It is the default absolute
    precision at which the elements will be compared. By default, it is
    twice the default precision::

        sage: R.halting_prec()
        40

    However, both the default precision and the halting precision can be
    customized at the creation of the parent as follows:

        sage: S = ZpER(5, prec=10, halt=100)
        sage: S.default_prec()
        10
        sage: S.halting_prec()
        100

    One creates elements as usual::

        sage: a = R(17/42)
        sage: a
        ...00244200244200244201

        sage: R.random_element()  # random
        ...21013213133412431402

    Here we notice that 20 digits (that is the default precision) are printed.
    However, the computation model is designed in order to guarantee that more
    digits of `a` will be available on demand.
    This feature is reflected by the fact that, when we ask for the precision
    of `a`, the software answers `+\infty`::

        sage: a.precision_absolute()
        +Infinity

    Asking for more digits is achieved by the methods :meth:`at_precision_absolute`
    and :meth:`at_precision_relative`::

        sage: a.at_precision_absolute(30)
        ...?244200244200244200244200244201

    As a shortcut, one can use the bracket operator::

        sage: a[:30]
        ...?244200244200244200244200244201

    Of course, standard operations are supported::

        sage: b = R(42/17)
        sage: a + b
        ...03232011214322140002
        sage: a - b
        ...42311334324023403400
        sage: a * b
        ...00000000000000000001
        sage: a / b
        ...12442142113021233401
        sage: sqrt(a)
        ...20042333114021142101

    We observe again that only 20 digits are printed but, as before,
    more digits are available on demand::

        sage: sqrt(a)[:30]
        ...?142443342120042333114021142101

    .. RUBRIC:: Equality tests

    Checking equalities between relaxed `p`-adics is a bit subtle and can
    sometimes be puzzling at first glance.

    When the parent is created with ``secure=False`` (which is the
    default), elements are compared at the current precision, or at the
    default halting precision if it is higher::

        sage: a == b
        False

        sage: a == sqrt(a)^2
        True
        sage: a == sqrt(a)^2 + 5^50
        True

    In the above example, the halting precision is `40`; it is the
    reason why a congruence modulo `5^50` is considered as an equality.
    However, if both sides of the equalities have been previously
    computed with more digits, those digits are taken into account.
    Hence comparing two elements at different times can produce
    different results::

        sage: aa = sqrt(a)^2 + 5^50
        sage: a == aa
        True
        sage: a[:60]
        ...?244200244200244200244200244200244200244200244200244200244201
        sage: aa[:60]
        ...?244200244300244200244200244200244200244200244200244200244201
        sage: a == aa
        False

    This annoying situation, where the output of `a == aa` may change
    depending on previous computations, cannot occur when the parent is
    created with ``secure=True``.
    Indeed, in this case, if the equality cannot be decided, an error
    is raised::

        sage: S = ZpER(5, secure=True)
        sage: u = S.random_element()
        sage: uu = u + 5^50
        sage: u == uu
        Traceback (most recent call last):
        ...
        PrecisionError: unable to decide equality; try to bound precision

        sage: u[:60] == uu
        False

    .. RUBRIC:: Self-referent numbers

    A quite interesting feature with relaxed `p`-adics is the possibility to
    create (in some cases) self-referent numbers. Here is an example.
    We first declare a new variable as follows::

        sage: x = R.unknown()
        sage: x
        ...?.0

    We then use the method :meth:`set` to define `x` by writing down an equation
    it satisfies::

        sage: x.set(1 + 5*x^2)
        True

    The variable `x` now contains the unique solution of the equation
    `x = 1 + 5 x^2`::

        sage: x
        ...04222412141121000211

    This works because the `n`-th digit of the right hand size of the
    defining equation only involves the `i`-th digits of `x` with `i < n`
    (this is due to the factor `5`).

    As a comparison, the following does not work::

        sage: y = R.unknown()
        sage: y.set(1 + 3*y^2)
        True
        sage: y
        ...?.0
        sage: y[:20]
        Traceback (most recent call last):
        ...
        RecursionError: definition looks circular

    Self-referent definitions also work with systems of equations::

        sage: u = R.unknown()
        sage: v = R.unknown()
        sage: w = R.unknown()

        sage: u.set(1 + 2*v + 3*w^2 + 5*u*v*w)
        True
        sage: v.set(2 + 4*w + sqrt(1 + 5*u + 10*v + 15*w))
        True
        sage: w.set(3 + 25*(u*v + v*w + u*w))
        True

        sage: u
        ...31203130103131131433
        sage: v
        ...33441043031103114240
        sage: w
        ...30212422041102444403

CC: @sagetrac-TristanVaccon @ThibautVerron @roed314 @saraedum

Component: padics

Author: Xavier Caruso

Branch: ad06299

Reviewer: David Roe

Issue created by migration from https://trac.sagemath.org/ticket/31108

xcaruso commented 3 years ago

Branch: u/caruso/relaxed

xcaruso commented 3 years ago

Commit: 7229cef

xcaruso commented 3 years ago

Description changed:

--- 
+++ 
@@ -95,13 +95,14 @@
 (this is due to the factor `5`).

 As a comparison, the following produces an error::
+
     sage: y = R()
     sage: y == 1 + 3*y^2
     Traceback (most recent call last):
     ...
     RecursionError: definition looks circular

-Previous self-referrent definitions also work with system of equations::
+Previous self-referent definitions also work with system of equations::

     sage: u = R(); v = R(); w = R()
xcaruso commented 3 years ago

New commits:

7229ceffirst implementation of lazy p-adics
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

1fb0d59rough implementation of square roots
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 7229cef to 1fb0d59

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 1fb0d59 to b69ba5a

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

b69ba5amethod selfref
xcaruso commented 3 years ago

Description changed:

--- 
+++ 
@@ -1,4 +1,4 @@
-We propose an implement of lazy p-adic numbers relaying on relaxed arithmetics proposed by van der Hoeven and al.
+We propose an implementation of lazy p-adic numbers relaying on relaxed arithmetics proposed by van der Hoeven and al.

 Here is a small demo:

@@ -104,7 +104,9 @@

 Previous self-referent definitions also work with system of equations::

-    sage: u = R(); v = R(); w = R()
+    sage: u = R.selfref()
+    sage: v = R.selfref()
+    sage: w = R.selfref()

     sage: u == 1 + 2*v + 3*w^2 + 5*u*v*w
     True
xcaruso commented 3 years ago

New commits:

b69ba5amethod selfref
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from b69ba5a to 4044cad

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

4044cadrewrite in cython+flint (not finished)
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

5649641cdef jump + improve shift
688061eimplement Qp and fraction_field
f339732more shared variables
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 4044cad to f339732

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

be24dadimplement division
a35d5d4coercion in equality test
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from f339732 to a35d5d4

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from a35d5d4 to 707287b

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

707287bc helper file
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

b130f14fix bugs in element_constructor
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 707287b to b130f14

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from b130f14 to ac5dbbc

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

d119381implement square root
ac5dbbcimplement lift method, so print_mode='terse' now works
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

4f2f8e6explicit error codes
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from ac5dbbc to 4f2f8e6

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

e23dccfbetter with files added
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 4f2f8e6 to e23dccf

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

f2ccf72fix bugs
596280eupdate tutorial
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from e23dccf to 596280e

xcaruso commented 3 years ago

Description changed:

--- 
+++ 
@@ -5,9 +5,9 @@

This module provides basic support for lazy p-adic integers.

@@ -18,20 +18,17 @@ sage: R.random_element() # random ...21013213133412431402

-By default, 20 digits of precision are computed (and printed). -If more (or less) digits are needed, one can specify it as follows::

-We observe again that only 20 digits are computed. -If we need more, we have to create a new variable:: +We observe again that only 20 digits are printed, even if the precision +on the operands is higher:: +

-Equality test works but equality is only checked up to the -minimal current precision of the elements::

@@ -96,13 +91,13 @@

As a comparison, the following produces an error::

-Previous self-referent definitions also work with system of equations:: +Self-referent definitions also work with systems of equations::

 sage: u = R.selfref()
 sage: v = R.selfref()

@@ -110,15 +105,15 @@

 sage: u == 1 + 2*v + 3*w^2 + 5*u*v*w
 True
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 596280e to 202e160

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

202e160square roots when p=2
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 202e160 to 7e6e696

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

d45dffdimplement Teichmuller lifts
7e7fb0achange construction of Teichmuller lifts
7e6e696implement method residue
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

1799646write templates
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 7e6e696 to 1799646

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 1799646 to b2b377a

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

90abf2fmake API and implementation coherent
d9aa5d7cdef inline classes
2687c79move element_constructor
b2b377areorganize classes
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from b2b377a to 2f51842

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

56001b4slices
4818fefinit jump
30d36dabounded and unbounded elements
b692710replace jump_* by at_precision_*
6ea0a8cchange behaviour of `_repr_`, precision_absolute and precision_absolute
bca6f7bdo not jump at initialization for unbounded elements
2f51842fix some bugs
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

357721dallow to set first digits in self referent numbers
70c752cremove last occurrence of _zeroprec
0777dd7implement add_bigoh and lift_to_precision
1959203fix teichmuller and residue
f395104update tutorial
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 2f51842 to f395104

xcaruso commented 3 years ago

Description changed:

--- 
+++ 
@@ -3,11 +3,27 @@
 Here is a small demo:

-This module provides basic support for lazy p-adic integers. +The model for lazy elements is quite different from any of the +other types of p-adics. In addition to storing a finite +approximation, one also stores a method for increasing the +precision. + +Lazy p-adic rings are created by the constructor :func:ZpL::

 sage: R = ZpL(5, print_mode="digits")
 sage: R
 5-adic Ring with lazy precision

+ +Observe that the precision is not capped on R:: +

@@ -18,52 +34,79 @@ sage: R.random_element() # random ...21013213133412431402

-By default, 20 digits of precision are printed. -If more digits are needed, one can increase the precision by using the -meth:jump:: +Here we notice that 20 digits (that is the default precision) are printed. +However, the computation model is designed in order to guarantee that more +digits of a are available on demand. +This feature is reflected by the fact that, when we ask for the precision +of a, the software answers \infty::

-Standard operations are implemented:: +Asking for more digits is achieved by the methods :meth:at_precision_absolute +and :meth:at_precision_relative:: +

-We observe again that only 20 digits are printed, even if the precision -on the operands is higher:: +We observe again that only 20 digits are printed but, as before, +more digits are available on demand::

-If more digits are needed, we have to create a new variable:: +Checking equalities between lazy p-adics is a bit subtle can could +sometimes be puzzling at first glance. +Actually, when it is obvious (from the previous computations) that +the two sides of the equality are different, everything works well::

-Note that this automatically increases the precision on a and b:: +On the contrary, when the two numbers we want to compare are indeed +equal, it is not possible to conclude after a finite amount of +computations. In this case, an error is raised::

@@ -75,12 +118,13 @@ sage: x ...?.0

-We then write down an equation satisfied by x:: +We then use the method :meth:set to define x by writing down an equation +it satisfies::

-The variable x now contains the unique solution of the above equation:: +The variable x now contains the unique solution of the equation +x = 1 + 5 x^2::

 sage: x
 ...04222412141121000211

@@ -89,10 +133,13 @@ defining equation only involves the i-th digits of x with i < n (this is due to the factor 5).

-As a comparison, the following produces an error:: +As a comparison, the following does not work::

 sage: y = R.selfref()
nbruin commented 3 years ago
comment:23

A note about operator choice: While the @ notation is definite informative and cute when considered stand-alone, it doesn't fit very well with other places where the operator is used in python. Specifically, the operator was introduced upon request of the numpy people, to signify matrix multiplication. With that in mind @ has a rather strong association with (function) composition (which, in a specific way, is also how you could read @ decorators to function definitions).

There's a reasonable chance that @ will at some point be integrated in the coercion framework. The use for @ proposed here will cause problems. I'd suggest just using a named method to cap precision.

Otherwise, this package looks really cool! It would be great to see some performance specs at some point.

xcaruso commented 3 years ago
comment:24

I'm a bit disappointed but I note your argument for @ and will probably disallow this construction.

About performances, it's of course not as fast as other types of p-adics but I would say that it's already reasonable (although optimizations are still possible):

    sage: R = ZpL(5)
    sage: a = R.random_element()
    sage: b = R.random_element()
    sage: c = 1 + 5*R.random_element()
    sage: %timeit _ = (a + b)@10000
    848 µs ± 614 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
    sage: %timeit _ = (a * b)@10000
    9.59 ms ± 116 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    sage: %timeit _ = (a / b)@10000
    11.7 ms ± 215 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    sage: %timeit _ = c.sqrt()@10000
    13.7 ms ± 240 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Compare with:

    sage: R = Zp(5, prec=10000)
    sage: a = R.random_element()
    sage: b = R.random_element()
    sage: c = 1 + 5*R.random_element()
    sage: %timeit _ = a + b
    1.21 µs ± 9.2 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
    sage: %timeit _ = a * b
    87.8 µs ± 130 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
    sage: %timeit _ = a / b
    667 µs ± 4.08 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
    sage: %timeit _ = c.sqrt()
    27.4 ms ± 125 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

EDIT: use %timeit instead of %time

xcaruso commented 3 years ago
comment:25

Maybe some comments on the timings. The "bad" performances of lazy p-adic come from two different sources:

Although I cannot do much to get rid of the first problem, I can certainly do something for the second one: instead of working with digits, I can work with blocks of digits (limbs) of the correct size. Actually, I want to do this but this requires many changes in my current code and I planned to delay this for another ticket.

In any case, when p becomes larger, the differences of performances between lazy p-adics and classical p-adics tend to get smaller.

    sage: p = next_prime(2^63)
    sage: R = ZpL(p)

    sage: a = R.random_element()
    sage: b = R.random_element()
    sage: c = 1 + p*R.random_element()
    sage: %timeit _ = (a + b)@10000
    3.51 ms ± 407 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    sage: %timeit _ = (a * b)@10000
    73.7 ms ± 189 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    sage: %timeit _ = (a / b)@10000
    83.4 ms ± 471 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    sage: %timeit _ = c.sqrt()@10000
    160 ms ± 3.31 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

and:

    sage: p = next_prime(2^63)
    sage: R = Zp(p, prec=10000)

    sage: a = R.random_element()
    sage: b = R.random_element()
    sage: c = 1 + p*R.random_element()
    sage: %timeit _ = a + b
    86.4 µs ± 578 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
    sage: %timeit _ = a * b
    9.04 ms ± 63.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    sage: %timeit _ = a / b
    78.8 ms ± 415 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    sage: %timeit _ = c.sqrt()
    672 ms ± 1.23 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from f395104 to 7fc3fb3

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

72f2f70documentation in linkage files
a632ec9move tutorial to the documentation of ZpL
3199f73printer is always defined
5aa0af9try to accelerate addition and subtraction
baeeb10documentation in generic_nodes
5316beefix doctest
7fc3fb3make TestSuite pass (except pickling)
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Changed commit from 7fc3fb3 to bac0fd1

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

64a41c2use prefix for element_classes
47d3921implement pickling
25f2879improve random element
92758cbsimplify a bit pickling
bac0fd1make TestSuite pass for elements
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 3 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

a995e30some documentation
cfb242dmore documentation
ab8111dremove the @ operator