Relaxed p-adics

[caruso]

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

[caruso]

New commits:

​7229cef

first implementation of lazy p-adics

[git]

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

​1fb0d59

rough implementation of square roots

[git]

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

​b69ba5a

method selfref

[caruso]

New commits:

​b69ba5a

method selfref

[git]

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

​4044cad

rewrite in cython+flint (not finished)

[git]

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

​5649641	cdef jump + improve shift
​688061e	implement Qp and fraction_field
​f339732	more shared variables

[git]

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

​be24dad	implement division
​a35d5d4	coercion in equality test

[git]

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

​707287b

c helper file

[git]

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

​b130f14

fix bugs in element_constructor

[git]

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

​d119381	implement square root
​ac5dbbc	implement lift method, so print_mode='terse' now works

[git]

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

​4f2f8e6

explicit error codes

[git]

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

​e23dccf

better with files added

[git]

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

​f2ccf72	fix bugs
​596280e	update tutorial

[git]

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

​202e160

square roots when p=2

[git]

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

​d45dffd	implement Teichmuller lifts
​7e7fb0a	change construction of Teichmuller lifts
​7e6e696	implement method residue

[git]

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

​1799646

write templates

[git]

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

​90abf2f	make API and implementation coherent
​d9aa5d7	cdef inline classes
​2687c79	move element_constructor
​b2b377a	reorganize classes

[git]

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

​56001b4	slices
​4818fef	init jump
​30d36da	bounded and unbounded elements
​b692710	replace jump_* by at_precision_*
​6ea0a8c	change behaviour of _repr_, precision_absolute and precision_absolute
​bca6f7b	do not jump at initialization for unbounded elements
​2f51842	fix some bugs

[git]

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

​357721d	allow to set first digits in self referent numbers
​70c752c	remove last occurrence of _zeroprec
​0777dd7	implement add_bigoh and lift_to_precision
​1959203	fix teichmuller and residue
​f395104	update tutorial

[nbruin]

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.

[caruso]

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

[caruso]

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

• the underlying arithmetics looses a factor log(precision) in the complexity
• I'm working on digits, which means that I basically only manipulate integers between 0 and p (here p=5) and encode them on longs.

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)

[git]

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

​72f2f70	documentation in linkage files
​a632ec9	move tutorial to the documentation of ZpL
​3199f73	printer is always defined
​5aa0af9	try to accelerate addition and subtraction
​baeeb10	documentation in generic_nodes
​5316bee	fix doctest
​7fc3fb3	make TestSuite pass (except pickling)

[git]

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

​64a41c2	use prefix for element_classes
​47d3921	implement pickling
​25f2879	improve random element
​92758cb	simplify a bit pickling
​bac0fd1	make TestSuite pass for elements

[git]

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

​a995e30	some documentation
​cfb242d	more documentation
​ab8111d	remove the @ operator

[git]

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

​ccbdfee	more doctest + fix bugs
​710dead	more doctests
​501c26c	doctest for all classes handling lazy p-adics
​30089ba	fix doctests

[git]

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

​7b5b2c3	Merge branch 'develop' of https://github.com/sagemath/sage into relaxed
​e0acbd9	implement expansion in smallest mode

[git]

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

​0cf1d59

implement expansion in Teichmuller mode

[caruso]

I guess that the ticket is ready for a first review...

[git]

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

​9398e97

fix pyflakes

[roed]

Some issues that we talked about on video:

• Add some halting parameters to the parent
• Inconsistency with valuation and set method. Add a check parameter to set to ensure that the element satisfies the given equation (even with provided valuation and/or digits) or have set return True or False depending on whether the equation is satisfied.
• Figure out segfaults in gcd computations, conversions
• Change selfref to unknown

For another ticket: add more transcendental functions like exp and log

[roed]

More comments as I read through the code.

• In digit_smallest you call fmpz_init but not fmpz_clear. I think this is a memory leak. I may have missed other leaks.... I also think that the documentation of digit_smallest should clarify that the input should already be in the range 0..p-1.
• For digit_sqrt it may be worth using FLINT's n_sqrtmod if p is word-size. You could also do some other easy cases in-line rather than calling out to Sage's integer_mod code. See square_root_mod_prime in sage/rings/finite_rings/integer_mod.pyx.
• You should add documentation to sage/misc/persist.pyx on how already_pickled and already_unpickled are intended to work, and some tests to check that it's working correctly.
• When converting an element between a lazy p-adic ring and its field of fractions it looks in _element_constructor_ that you create a bounded element. Why is this necessary? It also looks like you don't check that the primes are the same....
• The error codes in padic_lazy_error.pyx and next_c() aren't the same.
• In _repr_ it looks like some print modes aren't included (e.g. terse and val-unit). I would guess it to be easier to use code from padic_printing.pyx here.
• Last week we talked about adding some parameters about how halting is handled (e.g. a strict mode that raises an error whenever unbounded elements are compared). __nonzero__ should also be updated. I think the precision behavior for __nonzero__ should be the same as for == 0; currently it will return True rather than raising a precision error if it looks like 0 at known precision.
• I found the use of permissive=None a bit strange in at_precision_absolute. I think if a user specifies infinite precision, permissive should be treated as False (so, replace if permissive is False with if not permissive.
• Shouldn't at_precision_relative also support prec=infinity?
• You can probably remove some code by changing if error: raise_error(error) to just raise_error(error) since this has no effect if error is 0.
• _sub_ has a bug: currently 0-1 will give 1.

More comments to come....

[roed]

• In next_c() for random elements, do you need to digit_clear(r) in both branches of the if statement?
• It would be good to have some tests for slices of slices.
• In next_c() for LazyElement_div, if definition is None and you set it, don't you still want to execute the else clause? If not, you can just return self._bootstrap_c() rather than if error: return error. Same for LazyElement_sqrt.
• In the __init__ method for LazyElement_sqrt, shouldn't you check that the valuation is even?
• You should probably delete the next method on LazyElement_sqrt (I assume it was for debugging; it's not documented).

[caruso]

Replying to roed and roed:

• In digit_smallest you call fmpz_init but not fmpz_clear. I think this is a memory leak. I may have missed other leaks.... I also think that the documentation of digit_smallest should clarify that the input should already be in the range 0..p-1.
• You should add documentation to sage/misc/persist.pyx on how already_pickled and already_unpickled are intended to work, and some tests to check that it's working correctly.
• The error codes in padic_lazy_error.pyx and next_c() aren't the same.
• I found the use of permissive=None a bit strange in at_precision_absolute. I think if a user specifies infinite precision, permissive should be treated as False (so, replace if permissive is False with if not permissive.
• Shouldn't at_precision_relative also support prec=infinity?
• You can probably remove some code by changing if error: raise_error(error) to just raise_error(error) since this has no effect if error is 0.
• _sub_ has a bug: currently 0-1 will give 1.
• In next_c() for random elements, do you need to digit_clear(r) in both branches of the if statement?
• It would be good to have some tests for slices of slices.
• You should probably delete the next method on LazyElement_sqrt (I assume it was for debugging; it's not documented).

Fixed.

• In _repr_ it looks like some print modes aren't included (e.g. terse and val-unit). I would guess it to be easier to use code from padic_printing.pyx here.

They are included. I agree that we should include this code in padic_printing.pyx but it does not sounds easy.

• When converting an element between a lazy p-adic ring and its field of fractions it looks in _element_constructor_ that you create a bounded element. Why is this necessary? It also looks like you don't check that the primes are the same....

It's true than I call LazyElement_bound but I pass in boundprec=maxordp so that the resulting element remains unbounded. It's a convenient way to do copies. (I cannot return the same element because the parent changes.)

• In next_c() for LazyElement_div, if definition is None and you set it, don't you still want to execute the else clause? If not, you can just return self._bootstrap_c() rather than if error: return error. Same for LazyElement_sqrt.

Well, _bootstrap_c already computes the first significant digit, so I think it's fine like this. I changed if error: return error to return self._bootstrap_c().

• In the __init__ method for LazyElement_sqrt, shouldn't you check that the valuation is even?

It's checked in the bootstrap function. I think we cannot check it in the __init__ method because the valuation might be not known for self-referent numbers.

---

New commits:

​54c46e2

Fixing typos, grammar and whitespace

[caruso]

And here the comments I have not addressed yet:

• Last week we talked about adding some parameters about how halting is handled (e.g. a strict mode that raises an error whenever unbounded elements are compared). __nonzero__ should also be updated. I think the precision behavior for __nonzero__ should be the same as for == 0; currently it will return True rather than raising a precision error if it looks like 0 at known precision.
• Inconsistency with valuation and set method. Add a check parameter to set to ensure that the element satisfies the given equation (even with provided valuation and/or digits) or have set return True or False depending on whether the equation is satisfied.
• Figure out segfaults in gcd computations, conversions
• Change selfref to unknown
• For digit_sqrt it may be worth using FLINT's n_sqrtmod if p is word-size. You could also do some other easy cases in-line rather than calling out to Sage's integer_mod code. See square_root_mod_prime in sage/rings/finite_rings/integer_mod.pyx.

[git]

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

​5a5bdd9

Address David's remarks

[mkoeppe]

Sage development has entered the release candidate phase for 9.3. Setting a new milestone for this ticket based on a cursory review of ticket status, priority, and last modification date.

[git]

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

​30c5641	Merge branch 'develop' into relaxed
​429c83e	Merge branch 'develop' into relaxed
​9b0d5d1	selfref -> unknown
​b520df3	halting precision

[git]

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

​cfad0f2	improve equality test and computation of valuation
​a4a72a5	fix conversion and other bugs

[git]

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

​b21c36d

use fmpz_sqrtmod

[caruso]

I think that all your comments have been addressed now.

We still need to discuss the choice of the name for the classes implemented in this ticket (cf discussion on zulip).

[git]

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

​7f49e18

add keyword secure in the parent

[git]

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

​24ef236

lazy -> relaxed

[caruso]

I changed the naming (Relaxed instead of Lazy, ZpER instead of ZpL).

Ticket ready for a second review.

[git]

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

​2b14c83	update README.txt
​df13f1e	fix pyflakes

[git]

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

​50c9719

simplify exposition of relaxed hierarchy

[roed]

Looks good to me! I think any additional issues can be addressed in followup tickets.

[git]

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

​e30393b

minor typos

[caruso]

Great, thanks! I set again the ticket to positive review, I've just fixed two typos.

[git]

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

​ad06299

typo

[mkoeppe]

Looks like this broke the build --> #31903