Opened 5 years ago

Closed 5 years ago

## #23043 closed enhancement (fixed)

Reported by: Owned by: caruso major sage-8.0 padics logarithm, sd86.5 roed, saraedum Xavier Caruso David Roe N/A 16d8757 16d8757c92dccab7cb8164700217c6710905cc95

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: 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
```

### comment:2 Changed 5 years ago by git

• Commit set to e83902f25c3fa2dd75f4b7180cd793f4245fb194

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

 ​e83902f `Small fixes`

### comment:3 Changed 5 years ago by caruso

• Description modified (diff)

### comment:4 Changed 5 years 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 5 years ago by caruso

Ok, I'll have a look.

### comment:6 Changed 5 years ago by git

• 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 5 years ago by git

• Commit changed from a7d40cc4d4e59c6c51befbf6285aeea504c4d141 to 33817875c41e8e182da9a7e911d7fad8562cd2ae

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

 ​3381787 `Small optimization`

### comment:8 Changed 5 years ago by caruso

• Description modified (diff)
• Status changed from new to needs_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 5 years ago by git

• 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 5 years ago by git

• Commit changed from f1c7c2e0cba402c402f81001a357c0c986a975d1 to af9c989bdc87c833ca2963745092e3600264699a

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

 ​af9c989 `Doctest issue fixed`

### comment:11 Changed 5 years 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 5 years ago by git

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

 ​47ce6c0 `Faster algorithm (hopefully)`

### comment:13 Changed 5 years 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 5 years ago by git

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 p-adic logarithm`

### comment:15 Changed 5 years ago by caruso

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 p-adic logarithm`

### comment:16 Changed 5 years 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 5 years ago by git

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

 ​b475974 `Disable signal handling while memory allocation`

### comment:18 follow-up: ↓ 19 Changed 5 years 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 5 years ago by caruso (previous) (diff)

### comment:19 in reply to: ↑ 18 Changed 5 years ago by pbruin

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: ↓ 21 Changed 5 years 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 5 years ago by roed

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 5 years ago by roed

I'm building it now to run tests.

### comment:24 Changed 5 years 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:

 ​16d8757 `Adding limits to long tests in padic_extension_leaves, fix typo`

### comment:25 Changed 5 years ago by caruso

• Status changed from needs_review to positive_review

Ok, Thanks!

### comment:26 Changed 5 years 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.