Opened 2 years ago
Closed 2 years ago
#23043 closed enhancement (fixed)
Fast padic logarithm
Reported by:  caruso  Owned by:  

Priority:  major  Milestone:  sage8.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 )
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 2 years ago by
 Branch set to u/caruso/padiclog
comment:2 Changed 2 years ago by
 Commit set to e83902f25c3fa2dd75f4b7180cd793f4245fb194
comment:3 Changed 2 years ago by
 Description modified (diff)
comment:4 Changed 2 years ago by
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 2 years ago by
Ok, I'll have a look.
comment:6 Changed 2 years ago by
 Commit changed from e83902f25c3fa2dd75f4b7180cd793f4245fb194 to a7d40cc4d4e59c6c51befbf6285aeea504c4d141
Branch pushed to git repo; I updated commit sha1. New commits:
a7d40cc  Incorporate fast logarithm in ZpCA and ZpFM

comment:7 Changed 2 years ago by
 Commit changed from a7d40cc4d4e59c6c51befbf6285aeea504c4d141 to 33817875c41e8e182da9a7e911d7fad8562cd2ae
Branch pushed to git repo; I updated commit sha1. New commits:
3381787  Small optimization

comment:8 Changed 2 years ago by
 Description modified (diff)
 Status changed from new to needs_review
Patch ready for review.
Note that I've changed the behaviour of log for fixedmodulus 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 nonunits (while the absolute precision is fixed in the parent in that case). Please argue if you don't agree :).
comment:9 Changed 2 years ago by
 Commit changed from 33817875c41e8e182da9a7e911d7fad8562cd2ae to f1c7c2e0cba402c402f81001a357c0c986a975d1
Branch pushed to git repo; I updated commit sha1. New commits:
f1c7c2e  Memory leak problems fixed

comment:10 Changed 2 years ago by
 Commit changed from f1c7c2e0cba402c402f81001a357c0c986a975d1 to af9c989bdc87c833ca2963745092e3600264699a
Branch pushed to git repo; I updated commit sha1. New commits:
af9c989  Doctest issue fixed

comment:11 Changed 2 years ago by
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 2 years ago by
 Commit changed from af9c989bdc87c833ca2963745092e3600264699a to 47ce6c0dd0cad23ad777fa855e37abb6a0d8de74
Branch pushed to git repo; I updated commit sha1. New commits:
47ce6c0  Faster algorithm (hopefully)

comment:13 Changed 2 years ago by
 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 nonoptimized 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 witha//p
working.
 It would be good to use
sig_malloc
andsig_free
fromcysignals/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 padic logs and compare. Adding someTestSuite
tests (by implementing_test_log
for example; seepadic_generic.pyx
) would also be good, checking thatlog(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 padics. It would be cool to get this positively reviewed soon!
comment:14 Changed 2 years ago by
 Commit changed from 47ce6c0dd0cad23ad777fa855e37abb6a0d8de74 to 630ad0290ce6b6438b365ee9ade2b5bbd81d7105
Branch pushed to git repo; I updated commit sha1. New commits:
4e3a751  Added keywork 'algorithm'

13b1316  Switch to generic algorithm when p does not fit in a long

7472fc1  Switch to generic algorithm when p does not fit in a long (second try) + clean doctests

145da3c  More doctests

630ad02  Added TestSuite for padic logarithm

comment:15 Changed 2 years ago by
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:
4e3a751  Added keywork 'algorithm'

13b1316  Switch to generic algorithm when p does not fit in a long

7472fc1  Switch to generic algorithm when p does not fit in a long (second try) + clean doctests

145da3c  More doctests

630ad02  Added TestSuite for padic logarithm

comment:16 Changed 2 years ago by
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, callsig_unblock()
. This will duplicate the effect ofsig_malloc
incysignals/memory.pxd
.
comment:17 Changed 2 years ago by
 Commit changed from 630ad0290ce6b6438b365ee9ade2b5bbd81d7105 to b47597412e0db06fd2a21d70736efef8c97866c9
Branch pushed to git repo; I updated commit sha1. New commits:
b475974  Disable signal handling while memory allocation

comment:18 followup: ↓ 19 Changed 2 years ago by
 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?
comment:19 in reply to: ↑ 18 Changed 2 years ago by
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 followup: ↓ 21 Changed 2 years ago by
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 2 years ago by
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:23 Changed 2 years ago by
 Branch changed from u/caruso/padiclog to u/roed/padiclog
comment:24 Changed 2 years ago by
 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:
16d8757  Adding limits to long tests in padic_extension_leaves, fix typo

comment:26 Changed 2 years ago by
 Branch changed from u/roed/padiclog to 16d8757c92dccab7cb8164700217c6710905cc95
 Resolution set to fixed
 Status changed from positive_review to closed
Branch pushed to git repo; I updated commit sha1. New commits:
Small fixes