cleaning of linbox for dense integer matrices

[videlec]

Various cleaning to the linbox interface (sage.libs.linbox) and Sage dense integer matrices (sage.matrix.matrix_integer_dense).

• implement a more direct linbox/flint interface in sage/libs/linbox/linbox_flint_interface.pyx (supersedes part of linbox-sage.C inside linbox)

• remove all mpz pointers from Matrix_integer_dense (that is the attributes bint _initialized_mpz, mpz_t * _entrie and mpz_t ** _rows) as well as the associated constructors/destructors (int _init_mpz(self), int _init_mpz_impl(self), int _init_linbox(self), void _dealloc_mpz(self))

• fix the following error

sage: matrix(ZZ, 2).det(algorithm='padic')
Traceback (most recent call last):
...
ImportError: No module named matrix_integer_dense_hnf

• add more examples and tests to the documentation

(see #22872 for the linbox task ticket)

CC: @ClementPernet @sagetrac-dlucas

Component: packages: standard

Keywords: linbox

Author: Vincent Delecroix

Branch/Commit: 0b082c2

Reviewer: Clément Pernet

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

[videlec]

Branch: u/vdelecroix/22924

[videlec]

Commit: 38d9da3

[videlec]

New commits:

cca8edf	fix headers
38d9da3	Cleaning linbox interface of Matrix_integer_dense

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

86d8af3

Cleaning linbox interface of Matrix_integer_dense

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 38d9da3 to 86d8af3

[videlec]

Description changed:

--- 
+++ 
@@ -17,3 +17,4 @@

 - add more examples and tests to the documentation

+(see #22872 for the linbox task ticket)

[videlec]

Description changed:

--- 
+++ 
@@ -1,6 +1,6 @@
 Various cleaning to the linbox interface (`sage.libs.linbox`) and Sage dense integer matrices (`sage.matrix.matrix_integer_dense`).

-- implement in pure Cython the interface to linbox in a new file `sage/libs/linbox/linbox_flint_interface.pyx`
+- implement a more direct linbox/flint interface in `sage/libs/linbox/linbox_flint_interface.pyx` (supersedes part of `linbox-sage.C` inside linbox)

 - remove all mpz pointers from `Matrix_integer_dense` (that is the attributes
   `bint _initialized_mpz`, `mpz_t * _entrie` and `mpz_t ** _rows`) as

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 86d8af3 to 4cae968

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

cccbe30	Cleaning linbox interface of Matrix_integer_dense
4cae968	22924: remove (instead of deprecate) code from linbox.pyx

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

09f6280

22924: get rid of spies

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 4cae968 to 09f6280

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 09f6280 to 69a140e

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

48cd83b	fix headers
3f5ec83	Cleaning linbox interface of Matrix_integer_dense
ef07ed5	22924: remove (instead of deprecate) code from linbox.pyx
69a140e	22924: get rid of spies

[videlec]

comment:9

(rebased on 8.0.b4)

[ClementPernet]

comment:10

Great. Let's start with a few remarks:

1. Why do need the dependency on #22886 ? That ticket is stalled and is therefore also blocking this one.
2. While we are changing a lot of stuff, I wouldn't mind getting rid of these ugly poly_linbox method, with a string parameter typ that determines whether to compute minpoly or charpoly. We should rather only have the 2 methods _minpoly_linbox_ and _charpoly_linbox_. There is not much duplicated code between these two, especially now after this clean-up.

[ClementPernet]

Reviewer: Clément Pernet

[videlec]

comment:12

Replying to @ClementPernet:

Great.

Thanks for having a look!

Let's start with a few remarks:

1. Why do need the dependency on #22886 ? That ticket is stalled and is therefore also blocking this one.

I think that Cython compilation fails otherwise (I am trying again right now on top of 8.0.beta4... will take some time since I need to recompile the whole Sage after playing with #22493).

2. While we are changing a lot of stuff, I wouldn't mind getting rid of these ugly poly_linbox method, with a string parameter typ that determines whether to compute minpoly or charpoly. We should rather only have the 2 methods _minpoly_linbox_ and _charpoly_linbox_. There is not much duplicated code between these two, especially now after this clean-up.

The ticket only intends to clean matrix_integer_dense.pyx. The two functions you mention are actually removed from matrix/matrix_integer_dense.pyx but kept in libs/linbox/linbox_flint_interface.pyx (but note that the argument there is an integer). I can get rid of the following function in the interface

# set p to the minpoly or the charpoly of A
cdef inline void linbox_fmpz_mat_poly(fmpz_poly_t p, fmpz_mat_t A, int minpoly)

Is that what you want?

[ClementPernet]

comment:13

Replying to @videlec:

I can get rid of the following function in the interface

# set p to the minpoly or the charpoly of A
cdef inline void linbox_fmpz_mat_poly(fmpz_poly_t p, fmpz_mat_t A, int minpoly)

Is that what you want?

Ok sorry for the confusion. My remark indeed applies to the new libs/linbox_flint_interface.pyx. Sure, it should definitely not be exposed through the new interface as it is an internal trick to factor out code. Then I also suggest to avoid this trick and simply write the two methods for the sake of clarity.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

55bdc0a	Cleaning linbox interface of Matrix_integer_dense
c8b95de	22924: remove (instead of deprecate) code from linbox.pyx
fcafa60	22924: get rid of spies
637fd87	22924: cleaning in linbox_flint_interface.pyx

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 69a140e to 637fd87

[videlec]

comment:15

No dependency on #22886 anymore (it is needed for cylinbox but it is fine for this ticket).

I cleaned the flint interface a bit more than what you asked for.

[videlec]

Changed dependencies from #22886 to none

[videlec]

comment:16

I was wrong (and you were right). There is still def _poly_linbox(self, typ, var='x'): as a method of Matrix_integer_dense.

[videlec]

comment:17

By the way, I am adding tests and LinBox has trouble with the characteristic polynomial of the 0 x 0 matrix

sage: matrix([]).charpoly('y', 'linbox')
Traceback (most recent call last):
...
MemoryError: failed to allocate 2305843008945258512 bytes

Does it qualify as a bug? For now, I will add a special case in the interface LinBox/flint.

EDIT: and it hangs forever for minpoly.

[videlec]

comment:18

For minpoly it is written

        .. NOTE::

           Linbox charpoly disabled on 64-bit machines, since it hangs
           in many cases.

However, I do not see where this actually happens in the code. Do you know what this is?

[videlec]

comment:20

Un petit timing

sage: m = random_matrix(ZZ, 50)
sage: %timeit m._clear_cache(); p = m.charpoly(algorithm='linbox')
100 loops, best of 3: 16.8 ms per loop
sage: %timeit m._clear_cache(); p = m.charpoly(algorithm='flint')
10 loops, best of 3: 44.2 ms per loop
sage: %timeit m._clear_cache(); p = m.charpoly(algorithm='generic')
10 loops, best of 3: 152 ms per loop

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 637fd87 to fde1353

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

fde1353

22924: clean minpoly/charpoly in Matrix_integer_dense

[videlec]

comment:22

All right: Matrix_integer_dense cleaned and interface fixed on corner cases with commit fde1353.

[videlec]

comment:23

And patchbot is happy (excepted some timeout I am not responsible for).

[ClementPernet]

comment:24

Replying to @videlec:

For minpoly it is written

        .. NOTE::

           Linbox charpoly disabled on 64-bit machines, since it hangs
           in many cases.

However, I do not see where this actually happens in the code. Do you know what this is?

It's a super old comment. Some archeology on trac tickets should tell where it comes from but it's safe to remove it.

[ClementPernet]

comment:25

There are a couple of warnings that one could avoid:

[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:313:49: local variable 'res' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:389:46: local variable 'res' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:394:42: local variable 'res' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:398:42: local variable 'res' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:498:34: local variable 'g' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:532:32: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:534:16: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:562:30: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:564:16: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:717:33: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:736:33: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:756:35: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1110:39: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1127:39: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1150:39: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1167:39: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1184:39: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1217:35: local variable 'r' referenced before assignment
[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:1218:59: local variable 'r' referenced before assignment
[

Basically, the variables are passed as "out" arguments, and cython complains that they are not initialized. It is safe to force their value to 0 at declaration (e.g. int res = 0;) and this would silence out the warning.

[ClementPernet]

comment:26

Replying to @videlec:

By the way, I am adding tests and LinBox has trouble with the characteristic polynomial of the 0 x 0 matrix

sage: matrix([]).charpoly('y', 'linbox')
Traceback (most recent call last):
...
MemoryError: failed to allocate 2305843008945258512 bytes

Does it qualify as a bug? For now, I will add a special case in the interface LinBox/flint.

EDIT: and it hangs forever for minpoly.

This is definitely a bug. Reported upstream: https://github.com/linbox-team/linbox/issues/51

[ClementPernet]

comment:27

In linbox_flint_interface.pyx lines 84 and 157, you seem to use the method size of givvector to get the degree of the corresponding polynomial. While this is very likely to be correct in the current code, the specification do not require that the size of the vector storing the polynomial is ajusted to its degree. One should rather use Givaro::Poly1Dom::degree() instead.

[videlec]

comment:28

Replying to @ClementPernet:

There are a couple of warnings that one could avoid:

[sagelib-8.0.beta4] warning: sage/rings/finite_rings/element_givaro.pyx:313:49: local variable 'res' referenced before assignment
...

Basically, the variables are passed as "out" arguments, and cython complains that they are not initialized. It is safe to force their value to 0 at declaration (e.g. int res = 0;) and this would silence out the warning.

Yes, but these warnings have nothing to do with this ticket. Feel free to open another one to deal with rings/finite_rings (and put me in cc if you do).

[videlec]

comment:29

Replying to @ClementPernet:

In linbox_flint_interface.pyx lines 84 and 157, you seem to use the method size of givvector to get the degree of the corresponding polynomial. While this is very likely to be correct in the current code, the specification do not require that the size of the vector storing the polynomial is ajusted to its degree. One should rather use Givaro::Poly1Dom::degree() instead.

In LinBox code there are two occurrences of v.size() in examples/charpoly.C for a const Polynomial &v (lines 47 and 57). Idem in examples/minpoly.C (line 36). Are these also wrong?

How am I supposed to use degree? I was able to do it 7 lines in Cython... it looks extremely complicated to get such elementary quantity.

cdef GivaroIntegerRing ZZ
cdef GivaroDegree deg
cdef LinBoxIntegerPolynomialRing * R
R = new LinBoxIntegerPolynomialRing(ZZ)
R.degree(deg, q)
del R
cdef size_t length = deg.value() + 1

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from fde1353 to 6b45893

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

a75ec91	22924: better FIXME message
6b45893	22924: compute length from degree

[videlec]

comment:31

I implemented your remarks. Though, I do not think that 6b45893 is reasonable.

[ClementPernet]

comment:32

Replying to @videlec:

I implemented your remarks. Though, I do not think that 6b45893 is reasonable.

Ok this seems indeed quite tedious. After browsing through the current use of size() for degree in LinBox, I would say it's rather safe to assume that the vector has been resized and therefore that the size method would return the degree. I was probably being too pedantic, sorry. I suggest to drop 6b45893.

[ClementPernet]

comment:33

Replying to @videlec:

Yes, but these warnings have nothing to do with this ticket. Feel free to open another one to deal with rings/finite_rings (and put me in cc if you do).

My bad, I thought they were introduced in your new code. Will open a new ticket

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 6b45893 to a75ec91

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

[ClementPernet]

comment:35

Also, I see some new code not directly related to the subject of the ticket (e.g. in _rational_kernel_flint, rational_kernel_iml), but this is rather ok for me as it is documentation improvement on routines related to the ticket).

[videlec]

comment:36

Replying to @ClementPernet:

Replying to @videlec:

I implemented your remarks. Though, I do not think that 6b45893 is reasonable.

Ok this seems indeed quite tedious. After browsing through the current use of size() for degree in LinBox, I would say it's rather safe to assume that the vector has been resized and therefore that the size method would return the degree. I was probably being too pedantic, sorry. I suggest to drop 6b45893.

I dropped it.

It would be fine for me to use degree however it is horribly complicated. Why is there a class Degree that is a trivial wrapper around a int64_t? Why could we not use degree directly as a method of a polynomial?

[videlec]

comment:37

Replying to @ClementPernet:

Also, I see some new code not directly related to the subject of the ticket (e.g. in _rational_kernel_flint, rational_kernel_iml), but this is rather ok for me as it is documentation improvement on routines related to the ticket).

It is. _rational_kernel_iml now uses fmpz_mat_to_mpz_array that was introduced in the ticket.

[ClementPernet]

comment:38

Replying to @videlec:

Feel free to open another one to deal with rings/finite_rings (and put me in cc if you do).

Done in #22966.

[ClementPernet]

comment:39

Let me try an explanation for this mess (which I agree should be simplified):

Replying to @videlec:

Why is there a class Degree that is a trivial wrapper around a int64_t?

I would have said: so as to deal with the degree of the 0 polynomial which is -1 and not size()-1=0

Why could we not use degree directly as a method of a polynomial?

I assume because in Givaro (and the whole LinBox ecosystem), the container (here givvector does not know about the semantic). It is the domain of computation (here GivPoly1Dom who does).

[videlec]

comment:40

Replying to @ClementPernet:

Let me try an explanation for this mess (which I agree should be simplified):

Replying to @videlec:

Why is there a class Degree that is a trivial wrapper around a int64_t?

I would have said: so as to deal with the degree of the 0 polynomial which is -1 and not size()-1=0

In LinBox you have size(0) = 1!? Why would you store the zero at all?

Why could we not use degree directly as a method of a polynomial?

I assume because in Givaro (and the whole LinBox ecosystem), the container (here givvector does not know about the semantic). It is the domain of computation (here GivPoly1Dom who does).

I see. This way is indeed more flexible.

[videlec]

comment:41

Tell me if there are other improvements to be done on this ticket.