#23043 closed enhancement (fixed)

Fast p-adic logarithm

Reported by: caruso Owned by:
Priority: major Milestone: sage-8.0
Component: padics Keywords: logarithm, sd86.5
Cc: roed, saraedum Merged in:
Authors: Xavier Caruso Reviewers: David Roe
Report Upstream: N/A Work issues:
Branch: 16d8757 (Commits) Commit: 16d8757c92dccab7cb8164700217c6710905cc95
Dependencies: Stopgaps:

Description (last modified by caruso)

This ticket provides a fast C implementation of the logarithm for elements of \mathbb Z_p.

The gain of speed is significant:

    sage: from sage.rings.padics.padic_generic_element import pAdicGenericElement
    sage: generic_log = pAdicGenericElement.log

    sage: for i in range(2,11):
    ....:     print 'At precision 2^%s:' % i
    ....:     R = Zp(17, prec=2^i)
    ....:     a = 1 + 17*R.random_element()
    ....:     timeit('b = generic_log(a)')
    ....:
    At precision 2^2:
    625 loops, best of 3: 243 µs per loop
    At precision 2^3:
    625 loops, best of 3: 377 µs per loop
    At precision 2^4:
    625 loops, best of 3: 709 µs per loop
    At precision 2^5:
    625 loops, best of 3: 1.26 ms per loop
    At precision 2^6:
    125 loops, best of 3: 2.47 ms per loop
    At precision 2^7:
    125 loops, best of 3: 4.87 ms per loop
    At precision 2^8:
    25 loops, best of 3: 10.3 ms per loop
    At precision 2^9:
    25 loops, best of 3: 22.3 ms per loop
    At precision 2^10:
    5 loops, best of 3: 53.6 ms per loop

    sage: for i in range(2,16):
    ....:     print 'At precision 2^%s:' % i
    ....:     R = Zp(17, prec=2^i)
    ....:     a = 1 + 17*R.random_element()
    ....:     timeit('b = a.log()')
    ....:
    At precision 2^2:
    625 loops, best of 3: 4.36 µs per loop
    At precision 2^3:
    625 loops, best of 3: 5.62 µs per loop
    At precision 2^4:
    625 loops, best of 3: 8.92 µs per loop
    At precision 2^5:
    625 loops, best of 3: 16.8 µs per loop
    At precision 2^6:
    625 loops, best of 3: 32.2 µs per loop
    At precision 2^7:
    625 loops, best of 3: 62.8 µs per loop
    At precision 2^8:
    625 loops, best of 3: 143 µs per loop
    At precision 2^9:
    625 loops, best of 3: 333 µs per loop
    At precision 2^10:
    625 loops, best of 3: 975 µs per loop
    At precision 2^11:
    125 loops, best of 3: 2.49 ms per loop
    At precision 2^12:
    125 loops, best of 3: 7.12 ms per loop
    At precision 2^13:
    25 loops, best of 3: 19.3 ms per loop
    At precision 2^14:
    5 loops, best of 3: 54.6 ms per loop
    At precision 2^15:
    5 loops, best of 3: 158 ms per loop

Change History (26)

comment:1 Changed 15 months ago by caruso

  • Branch set to u/caruso/padiclog

comment:2 Changed 15 months ago by git

  • Commit set to e83902f25c3fa2dd75f4b7180cd793f4245fb194

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

e83902fSmall fixes

comment:3 Changed 15 months ago by caruso

  • Description modified (diff)

comment:4 Changed 15 months ago by chapoton

I have just launched my patchbot on this ticket.

Maybe you would be interested in having a look at #12560 ? I have done my best, but it still needs some care from an expert in padics and C/cython.

comment:5 Changed 15 months ago by caruso

Ok, I'll have a look.

comment:6 Changed 15 months ago by git

  • Commit changed from e83902f25c3fa2dd75f4b7180cd793f4245fb194 to a7d40cc4d4e59c6c51befbf6285aeea504c4d141

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

a7d40ccIncorporate fast logarithm in ZpCA and ZpFM

comment:7 Changed 15 months ago by git

  • Commit changed from a7d40cc4d4e59c6c51befbf6285aeea504c4d141 to 33817875c41e8e182da9a7e911d7fad8562cd2ae

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

3381787Small optimization

comment:8 Changed 15 months ago by caruso

  • Description modified (diff)
  • Status changed from new to needs_review

Patch ready for review.

Note that I've changed the behaviour of log for fixed-modulus padics: when a is not invertible, log(a) now raises an error. The reason for this is that taking logarithm decreases the absolute precision for non-units (while the absolute precision is fixed in the parent in that case). Please argue if you don't agree :-).

comment:9 Changed 15 months ago by git

  • Commit changed from 33817875c41e8e182da9a7e911d7fad8562cd2ae to f1c7c2e0cba402c402f81001a357c0c986a975d1

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

f1c7c2eMemory leak problems fixed

comment:10 Changed 15 months ago by git

  • Commit changed from f1c7c2e0cba402c402f81001a357c0c986a975d1 to af9c989bdc87c833ca2963745092e3600264699a

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

af9c989Doctest issue fixed

comment:11 Changed 15 months ago by caruso

I've improved a bit the algorithm (by making first the argument closer to 1 by raising - possibly several times - to the power p). New timings :

sage: for i in range(2,16):
....:     print 'At precision 2^%s:' % i
....:     R = Zp(17, prec=2^i)
....:     a = 1 + 17*R.random_element()
....:     timeit('b = a.log()')
....:     
At precision 2^2:
625 loops, best of 3: 5.51 µs per loop
At precision 2^3:
625 loops, best of 3: 7.83 µs per loop
At precision 2^4:
625 loops, best of 3: 11 µs per loop
At precision 2^5:
625 loops, best of 3: 13.6 µs per loop
At precision 2^6:
625 loops, best of 3: 21.9 µs per loop
At precision 2^7:
625 loops, best of 3: 40.5 µs per loop
At precision 2^8:
625 loops, best of 3: 88.7 µs per loop
At precision 2^9:
625 loops, best of 3: 184 µs per loop
At precision 2^10:
625 loops, best of 3: 526 µs per loop
At precision 2^11:
625 loops, best of 3: 1.57 ms per loop
At precision 2^12:
125 loops, best of 3: 4.59 ms per loop
At precision 2^13:
25 loops, best of 3: 13.2 ms per loop
At precision 2^14:
25 loops, best of 3: 38.3 ms per loop
At precision 2^15:
5 loops, best of 3: 123 ms per loop

comment:12 Changed 15 months ago by git

  • Commit changed from af9c989bdc87c833ca2963745092e3600264699a to 47ce6c0dd0cad23ad777fa855e37abb6a0d8de74

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

47ce6c0Faster algorithm (hopefully)

comment:13 Changed 15 months ago by roed

  • Status changed from needs_review to needs_work

Looks great!

  • p is allowed to be larger than an unsigned long. Presumably you should fall back on the non-optimized implementation then (the precision won't be large in this case, so it shouldn't really matter much).
  • I think that fixed modulus logs should still work even when a is not invertible, in analogy with a//p working.
  • It would be good to use sig_malloc and sig_free from cysignals/memory.pxd, but I don't know exactly how to achieve this in a C file.
  • It would be good to have additional tests. In particular, maybe have an algorithm keyword that allows falling back to the default implementation of p-adic logs and compare. Adding some TestSuite tests (by implementing _test_log for example; see padic_generic.pyx) would also be good, checking that log(xy) = log(x) + log(y) for example.
  • I should be fairly responsive for the next 9 days or so, since I'm at a Sage Days then working with Julian on p-adics. It would be cool to get this positively reviewed soon!

comment:14 Changed 15 months ago by git

  • Commit changed from 47ce6c0dd0cad23ad777fa855e37abb6a0d8de74 to 630ad0290ce6b6438b365ee9ade2b5bbd81d7105

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

4e3a751Added keywork 'algorithm'
13b1316Switch to generic algorithm when p does not fit in a long
7472fc1Switch to generic algorithm when p does not fit in a long (second try) + clean doctests
145da3cMore doctests
630ad02Added TestSuite for p-adic logarithm

comment:15 Changed 15 months ago by caruso

Great! Thanks for your review!

For sig_malloc and sig_free, I also do not know how to use them in a C code. Apart from that, I addressed all your requests.


New commits:

4e3a751Added keywork 'algorithm'
13b1316Switch to generic algorithm when p does not fit in a long
7472fc1Switch to generic algorithm when p does not fit in a long (second try) + clean doctests
145da3cMore doctests
630ad02Added TestSuite for p-adic logarithm

comment:16 Changed 15 months ago by roed

Great. Peter Bruin suggests the following for sig_malloc and sig_free:

  • Include cysignals/macros.h in your transcendentals.c file. This might be as simple as #include <cysignals/macros.h>, but we're not sure.
  • Before your mallocs, call sig_block() and afterward, call sig_unblock(). This will duplicate the effect of sig_malloc in cysignals/memory.pxd.

comment:17 Changed 15 months ago by git

  • Commit changed from 630ad0290ce6b6438b365ee9ade2b5bbd81d7105 to b47597412e0db06fd2a21d70736efef8c97866c9

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

b475974Disable signal handling while memory allocation

comment:18 follow-up: Changed 15 months ago by caruso

  • Status changed from needs_work to needs_review

Ok, done. But I am then wondering about other memory allocations (required by gmp). Are there correctly handled?

Last edited 15 months ago by caruso (previous) (diff)

comment:19 in reply to: ↑ 18 Changed 15 months ago by pbruin

Replying to caruso:

Ok, done. But I am then wondering about other memory allocations (required by gmp). Are there correctly handled?

I don't think there is a problem here; Sage sets the malloc/realloc/free functions to be used by GMP to sage_sig_malloc etc. (see sage/ext/memory.pyx). Actually, it could be an option for your code to use the same malloc as GMP does.

comment:20 follow-up: Changed 15 months ago by caruso

Fine.

Btw, what happens exactly when the computation is interrupted by KeyboardInterrupt. Is the memory allocated by the function automatically free? Or should we handle this?

comment:21 in reply to: ↑ 20 Changed 15 months ago by roed

Replying to caruso:

Fine.

Btw, what happens exactly when the computation is interrupted by KeyboardInterrupt. Is the memory allocated by the function automatically free? Or should we handle this?

I don't think it's automatically freed, but I don't know how to handle it in C. In Cython, you can use a try/finally block and free inside the finally. I think in this case we should just live with the possibility of a memory leak.

comment:22 Changed 15 months ago by roed

  • Keywords sd86.5 added

I'm building it now to run tests.

comment:23 Changed 15 months ago by roed

  • Branch changed from u/caruso/padiclog to u/roed/padiclog

comment:24 Changed 15 months ago by roed

  • Commit changed from b47597412e0db06fd2a21d70736efef8c97866c9 to 16d8757c92dccab7cb8164700217c6710905cc95
  • Reviewers set to David Roe

I've run all tests. Xavier, if you're happy with my changes (fixing a typo and changing some TestSuites? to not take forever), go ahead and set this to positive review.


New commits:

16d8757Adding limits to long tests in padic_extension_leaves, fix typo

comment:25 Changed 15 months ago by caruso

  • Status changed from needs_review to positive_review

Ok, Thanks!

comment:26 Changed 15 months ago by vbraun

  • Branch changed from u/roed/padiclog to 16d8757c92dccab7cb8164700217c6710905cc95
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.