[very positive review] charpoly of matrices over number fields is ridiculously slow right now (easy fixes exist)

[williamstein]

See the below examples that illustrate that charpoly(A), where A is over a number field, currently totally sucks. This has very bad implications for modular forms computations. Even switching to use PARI for charpoly would give a significant speedup (which is a lot faster than Magma, interestingly).

K.<a> = NumberField(x^2 + 17)
a^2
///
-17

n = 40
m = matrix(K, n, [(a+1)^randint(0,3) for _ in xrange(n^2)])
m
///
40 x 40 dense matrix over Number Field in a with defining polynomial x^2 + 17

m[0] # first row
///
(2*a - 16, 1, -14*a - 50, a + 1, 2*a - 16, 2*a - 16, a + 1, -14*a - 50, 1, -14*a - 50, 2*a - 16, a + 1, 1, a + 1, a + 1, a + 1, a + 1, -14*a - 50, -14*a - 50, 2*a - 16, 1, a + 1, 1, a + 1, 2*a - 16, -14*a - 50, 1, 2*a - 16, 2*a - 16, 1, a + 1, a + 1, 1, 2*a - 16, -14*a - 50, 2*a - 16, 2*a - 16, -14*a - 50, 2*a - 16, -14*a - 50)

time k = m*m
///
CPU time: 0.14 s,  Wall time: 0.27 s

time f=m.charpoly()
///
CPU time: 23.93 s,  Wall time: 26.22 s

NOTE: Sage should use PARI for charpoly's over number fields, but currently it doesn't. Notice how much faster PARI is at the same computation!

m._clear_cache()
g = pari(m)
time h = g.charpoly()
///
Time: CPU 2.52 s, Wall: 2.76 s

But a multimodular algorithm should blow away all of them. Student project? Basically, charpoly mod p is extremely fast in Sage, and one could reduce modulo primes of a number field, do the computation mod p, then lift, and get the correct result. I know of no implementations of this. A student project could be to implement this and tune it to be the fastest program in the world for charpoly over number fields.

time m.echelonize()
///
CPU time: 0.35 s,  Wall time: 0.35 s

%magma
R<x> := PolynomialRing(RationalField());
K<a> := NumberField(x^2 + 17);
n := 40;
m := MatrixAlgebra(K, n)![(1+a)^Random(0, 3) : i in [1..n^2]];
time k := m*m;
time f := CharacteristicPolynomial(m);
time e := EchelonForm(m);
///
Time: 0.000
Time: 15.220
Time: 0.150

Component: linear algebra

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

[85eec1a4-3d04-4b4d-b711-d4db03337c41]

Description changed:

--- 
+++ 
@@ -36,7 +36,7 @@
 CPU time: 23.93 s,  Wall time: 26.22 s

-NOTE: Sage should use PARI for charpoly's over number +NOTE: Sage should use PARI for charpoly's over number fields, but currently it doesn't. Notice how much faster PARI is at the same computation!

@@ -48,13 +48,12 @@ Time: CPU 2.52 s, Wall: 2.76 s

-But a <b>multimodular algorithm</b> should blow away all of them.  Student project?
+But a **multimodular algorithm** should blow away all of them.  Student project?
 Basically, charpoly mod p is extremely fast in Sage, and one could reduce modulo
-<i>primes of a number field</i>, do the computation mod p, then lift, and get the
+*primes of a number field*, do the computation mod p, then lift, and get the
 correct result.  I know of no implementations of this.  A student project could be
 to implement this and tune it to be the fastest program in the world for charpoly
 over number fields. 
-<br>

time m.echelonize()

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

comment:2

The following IRC transcript is relevant.

ncalexan-827: wstein-1606: is there a good way to find degree 1 primes in a number field.
ncalexan-827: ?
wstein-1606: theoretically or in sage right now?
ncalexan-827: wstein-1606: right now, but also in theory.
wstein-1606: for all but finitely many primes you factor the defining poly of the field and see if it splits completely.
dmharvey left the chat room.
wstein-1606: you can do that more quickly by taking the gcd mod p with x^p - x
ncalexan-827: Right.
wstein-1606: or something like that.
wstein-1606: I bet you get the idea.
ncalexan-827: Doesn't this only work if your ring of integers is generated by a single element?
wstein-1606: nick -- that's why I said "all but finitely many".
wstein-1606: For all but a couple primes the ring of integers is monogenic.
wstein-1606: just avoid primes dividing the discriminant.
wstein-1606: for those do the usual factorization algorithm.
ncalexan-827: wstein-1606: sure.  It sounds like a degree_one_primes() generator would be handy for multi modular algorithms, then.
wstein-1606: yep!
wstein-1606: But beware -- it is very very very hard to find any degree one primes if the field isn't cyclotomic.
wstein-1606: For a galois S_n extension of degree n, only 1/n! of the primes split completely (by cebo.)
ncalexan-827: Hmm.  So it's not even worth it?
wstein-1606: Well it's very worth it in the cyclotomic case.
wstein-1606: more generally, in the abelian case.
wstein-1606: E.g., for quadratic fields.
wstein-1606: and even 1/n! isn't bad if n is small, which it often is.

[JohnCremona]

comment:3

I needed exactly this for my work on Unimodular Integer Circulants (To appear in: Mathematics of Computation) or see my web page. At the time I used pari and have been meaning to try it in Sage too. Rather than take gcd(f,x^p-x) -- all mof p of course -- , I took Mod(x,fMod(1,p))^p and that worked pretty well. My worst example: deg(f) = 18, the first degree 1 prime was about 250 million and did not work for me so I sed the second one which was around 2 billion. This poly has a Galois group of order only 2(n/2)! ; if it been the generic n! I would probably still be waiting.

This suggests to me that for anything other than small degrees the implementation will need to use primes of higher degree.

[JohnCremona]

comment:4

Replying to @JohnCremona:

Sorry about bad formatting.

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

comment:5

See #2835 for a possible primes of degree one iterator implementation.

[aghitza]

Attachment: 903-charpoly_pari.patch.gz

[aghitza]

comment:6

I suggest that we implement the easy fix and get that merged, and open an enhancement ticket along the lines "implement multimodular algorithm for charpoly over number fields".

I'm attaching a patch with the easy fix. On my machine, William's example runs 15 times faster after the patch than before.

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

comment:7

I am not happy with this patch. I want to see an explicit doctest (not random) that has a polynomial defined over the number field, not over QQ.

I also want to see doctests for stacked relative fields -- I.e.

x = QQ['x'].gen()
K.<a> = NumberField(x^2 - 2)
L.<b> = K.extension(x^3 - a)
M.<c> = L.extension(x^2 - a*x + b)

and see it work for matrices over all those fields.

This is hard -- we need polynomials over arbitrary pari number fields, and it doesn't work for #2329 yet, so it shouldn't work here yet.

[aghitza]

Attachment: 903-charpoly_pari_doctests.patch.gz

apply after 903-charpoly_pari.patch

[aghitza]

comment:8

The point about the doctests is correct. I put up a new patch that adds the appropriate doctests and fixes a small bug (due to some unrelated weirdness that I don't want to deal with right now).

The doctests pass and I checked the results. It looks to me like it's working, so I'm not sure I agree with the last part of the last comment.

[mwhansen]

comment:9

The patches apply and pass tests for me. Nick, what are your thoughts?

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

Attachment: 903-ncalexan-charpoly-over-number-field-1.patch.gz

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

comment:10

I prefer the attached patch because it does not require string parsing, at least not in the code. It will be more robust to changes to PARI and uses the Sage library. It's shorter, too, and I prefer the factored-out function.

[mwhansen]

comment:11

Nick's patch applies and passes the tests. I think it is preferable since it avoids string parsing. I'm going to give it a positive review. Alex, if you disagree, feel free to comment and change it.

[aghitza]

comment:12

Yes, I much prefer Nick's code over mine. I was struggling to get things from Pari into Sage and the way I was doing things was obscure and rather fragile.

There is one problem I run into with this: 2 doctests in number_field_element.pyx fail. Probably an easy fix?

[aghitza]

Attachment: 903-small_fix.patch.gz

apply after Nick's patch

[aghitza]

comment:13

It was indeed an easy fix. Pari treats the t_POLMOD 0 in a different way than the others, so the code chocked on that. I've put up a patch that fixes this and adds a direct doctest for this behavior.

We should apply 903-ncalexan-charpoly-over-number-field-1.patch and 903-small_fix.patch. When this ticket gets closed we should make a new enhancement ticket containing the above discussion of a multimodular algorithm.

BTW, in 3.0.alpha5, before the patches:

sage: K.<a> = NumberField(x^2 + 17)
sage: n = 40
sage: m = matrix(K, n, [(a+1)^randint(0,3) for _ in xrange(n^2)])
sage: time f = m.charpoly('Z')
CPU times: user 59.69 s, sys: 0.01 s, total: 59.70 s
Wall time: 59.70

and after the patches:

sage: K.<a> = NumberField(x^2 + 17)
sage: n = 40
sage: m = matrix(K, n, [(a+1)^randint(0,3) for _ in xrange(n^2)])
sage: time f = m.charpoly('Z')
CPU times: user 3.96 s, sys: 0.02 s, total: 3.98 s
Wall time: 3.98

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

Attachment: 903-ncalexan-charpoly-over-number-field-2.patch.gz

[93749b3a-c0a4-47a7-b178-004ba08b0b97]

comment:14

903-ncalexan-charpoly-over-number-field-2.patch stands alone, and doesn't muck about with PARI at all: it turns out that polynomials over number fields work to the extent that special parsing is not necessary.

BTW, great team effort guys!

[aghitza]

comment:15

I love the end result, it's so clean and simple (I wish there was a way to delete my old patches from this ticket, I'm somewhat embarrassed :)

Anyway, it works great and fast and it should be merged.

[85eec1a4-3d04-4b4d-b711-d4db03337c41]

comment:16

Nick, Alex,

am I correct to understand that I must only apply 903-ncalexan-charpoly-over-number-field-2.patch?

it applies by itself and does build, and doctests fine.

Cheers,

Michael

[85eec1a4-3d04-4b4d-b711-d4db03337c41]

comment:17

Merged 903-ncalexan-charpoly-over-number-field-2.patch in Sage 3.0.alpha6 only. In case that is not what is intended please comment here and I will merge those patches.

Cheers,

Michael

[aghitza]

comment:18

Yes, that patch was the only one that needed merging.