Improve polynomial powering in positive characteristic

[bgrenet]

Before:

sage: R.<x> = PolynomialRing(GF(5), sparse=True)
sage: (1+x)^(5^10)
... (hangs for ever)

After:

sage: R.<x> = PolynomialRing(GF(5), sparse=True)
sage: (1+x)^(5^10)
x^9765625 + 1

Notes.

1. I've tried to be careful with quite extensive tests not to slow down the computations for "easy" cases. This explains the threshold values in the code.

2. This is related to inefficient polynomial powering for positive characteristic #7253 and make raising polynomials in characteristic p to large powers (and printing them) more efficient #12660. Though, there is another problem for non-sparse polynomials: Simply assigning the value and printing the polynomial take a very long time. Compare the timings for

   sage: R.<x> = ZZ[]
   sage: 1+x^5^10
   x^9765625 + 1

   and

   sage: R.<x> = GF(5)[]
   sage: 1+x^5^10
   x^9765625 + 1

   Thus for tickets inefficient polynomial powering for positive characteristic #7253 and make raising polynomials in characteristic p to large powers (and printing them) more efficient #12660 to be closed, it remains to handle this other problem, but this is not the purpose of this ticket.

Component: commutative algebra

Keywords: polynomial, powering

Author: Bruno Grenet

Branch/Commit: 41dd905

Reviewer: Jeroen Demeyer

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

[sagetrac-git]

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

a8ecb3e

Special code for positive characteristic

[sagetrac-git]

Commit: a8ecb3e

[sagetrac-git]

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

8f4125b

Typo in the doc

[sagetrac-git]

Changed commit from a8ecb3e to 8f4125b

[jdemeyer]

comment:4

Your added doctest still takes # long time (28s on my computer). You should shorten the test to take less than 1s: 5^8 instead of 5^10 would be better.

[jdemeyer]

comment:5

I don't get it, it seems that the doctest doesn't even use the new code. It's equally slow with or without this patch...

[sagetrac-fwclarke]

comment:6

I think it does use the code. But what's happening is that in the line

                return self.subs({x:x**(q*p)}) * self**r

x**(q*p) is taking (almost) all the time. And this comes from the section of the code dealing with powers of the polynomial generator, in particular the line

            return self.parent()(v, check=False)

where v is defined in the previous line as

                v = [R.zero()]*right + [R.one()]

with R being the base ring. It just takes a long time to create a non-sparse polynomial of this size.

[sagetrac-git]

Changed commit from 8f4125b to 139ac60

[sagetrac-git]

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

139ac60

Correct doctest

[sagetrac-git]

Changed commit from 139ac60 to 3088973

[sagetrac-git]

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

3088973

Special code only works for field... of course!

[bgrenet]

comment:9

Replying to @jdemeyer:

I don't get it, it seems that the doctest doesn't even use the new code. It's equally slow with or without this patch...

Right, because my doctest was wrong: The code is intended for sparse polynomials. I have not really been able yet to determine what is precisely going on for dense polynomial. Though, as indicated by fwclarke and mentioned in my note 2. in the description, there is another problem for dense polynomial, namely that creating a large polynomial over a finite field takes a very long time. It is not the purpose of this ticket to correct this, even though it should be corrected!

This patch thus improves the situation for sparse polynomials only. The code is indeed not used for dense polynomials over Zmod(k) for any k since specific code appears in src/sage/rings/polynomial/polynomial_template.pxi for this case.

Important note: I forgot to test whether the base ring is a field, which is of course mandatory. I correct this mistake in my latest commit.

[jdemeyer]

comment:10

You added 3 cases here:

if self.parent().base_ring().is_finite():
    return ...
if self.parent().is_sparse():
    ...
else:
    ...

You should add a doctest for each of these 3 cases (and perhaps replace the second if by elif for clarity).

[jdemeyer]

comment:11

Replying to @bgrenet:

Important note: I forgot to test whether the base ring is a field, which is of course mandatory.

Why is this "of course" mandatory? You need that the characteristic is a prime number, not that the base ring is a field.

[sagetrac-git]

Changed commit from 3088973 to 7760a43

[sagetrac-git]

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

7760a43

Correct algorithm + adapted doctests

[bgrenet]

comment:13

The algorithm I proposed was completely wrong: If p is the characteristic and n=p * q + r, I used the wrong identity f(x)^n = f(x^{p*q}) * f(x)^r instead of the correct f(x)^n = f(x<sup>p)</sup>q * f(x)^r.

In the commit, I correct this mistake, and in the same time I take jdemeyer's comments into account:

• I add doctests for the different possibilities. (Note that I wrote code for the case where the base ring is a finite field and the polynomials are in dense representation, though as I wrote above this code is not used since there is specific code in polynomial_template.pxi. Thus there is no doctest for this case.)
• "Of course", it is not "of course mandatory" for the base ring to be a finite field...

Furthermore, I improve the code not using self.subs(...), but rather doing the computations myself.

[jdemeyer]

comment:18

Please justify in the code why

# characteristic 2 is special

[jdemeyer]

comment:19

I suggest to remove the special case for finite fields. The generic code should work just fine for finite fields.

[jdemeyer]

comment:20

And I'm not too happy with this either:

self.base_ring().is_field()

Checking whether some complicated ring is a field might be an expensive operation. I prefer

self.base_ring() in IntegralDomains

[sagetrac-git]

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

7788e84	Merge branch 'u/bruno/polynomial_powering_positive_characteristic' of git://trac.sagemath.org/sage into t18547_poly_power_positive_char
3a452c7	trac 18547: characteristic 2 is not special
5a4965c	trac 18547: is_field() replaced by in IntegralDomains

[sagetrac-git]

Changed commit from 959947d to 5a4965c

[jdemeyer]

comment:22

I didn't say you had to remove this:

if ... or p.is_prime()

[bgrenet]

comment:23

Replying to @jdemeyer:

Please justify in the code why

# characteristic 2 is special

I had run some tests and my conclusion that was using generic_power was faster in characteristic 2. (Note that binary exponentiation is well suited to this characteristic.) I am not able to reproduce this behavior though, so I removed the special case.

Replying to @jdemeyer:

And I'm not too happy with this either:

self.base_ring().is_field()

Checking whether some complicated ring is a field might be an expensive operation. I prefer

self.base_ring() in IntegralDomains

Done.

Replying to @jdemeyer:

I suggest to remove the special case for finite fields. The generic code should work just fine for finite fields.

Actually, the current code is faster for finite fields. With the special case:

sage: R.<x> = PolynomialRing(GF(17),sparse=True)
sage: p = R.random_element(100)
sage: %timeit p^(17^5)
The slowest run took 11.13 times longer than the fastest. This could mean that an intermediate result is being cached 
1000 loops, best of 3: 159 µs per loop

Without the special case, with the same p:

%timeit p^(17^5)
The slowest run took 26.59 times longer than the fastest. This could mean that an intermediate result is being cached 
1000 loops, best of 3: 249 µs per loop

I may refactor the code though so that it is less redundant. Tell me what you think.

Replying to @jdemeyer:

I think this test is much too weak

I am happy if you find a good way to test this function!

[sagetrac-git]

Changed commit from 5a4965c to ba484df

[sagetrac-git]

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

ba484df

trac 18547: restore is_prime(), previously removed as an error

[bgrenet]

comment:25

Replying to @jdemeyer:

I didn't say you had to remove this:

if ... or p.is_prime()

Right, my mistake!

[jdemeyer]

comment:26

Replying to @bgrenet:

I am happy if you find a good way to test this function!

You have many cases and they should all be tested:

sage: R1 = PolynomialRing(GF(8, 'a'), 'x')
sage: R2 = PolynomialRing(GF(9, 'b'), 'x', sparse=True)
sage: R3 = PolynomialRing(R2, 'y')
sage: R4 = PolynomialRing(R1, 'y', sparse=True)
sage: for d in range(40):  # long time
....:     for R in [R1, R2, R3, R4]:
....:         a = R.random_element()
....:         assert a^d == generic_power(a, d)

[sagetrac-git]

Changed commit from ba484df to 41dd905

[sagetrac-git]

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

e8d11a6	Trac 18547: better doctests
41dd905	Trac 18547: Refactor and correct code

[bgrenet]

comment:28

Replying to @jdemeyer:

Replying to @bgrenet:

I am happy if you find a good way to test this function!

You have many cases and they should all be tested:

sage: R1 = PolynomialRing(GF(8, 'a'), 'x')
sage: R2 = PolynomialRing(GF(9, 'b'), 'x', sparse=True)
sage: R3 = PolynomialRing(R2, 'y')
sage: R4 = PolynomialRing(R1, 'y', sparse=True)
sage: for d in range(40):  # long time
....:     for R in [R1, R2, R3, R4]:
....:         a = R.random_element()
....:         assert a^d == generic_power(a, d)

I introduced your proposed doctest, with the following change: Since I put a threshold 20 on the exponent before using the new code (for efficiency reasons), the loop is for d in range(20,40).

Replying to @bgrenet:

Replying to @jdemeyer:

I suggest to remove the special case for finite fields. The generic code should work just fine for finite fields.

Actually, the current code is faster for finite fields.

Actually, the code was not correct for non-prime finite fields... (The special case should have been for prime finite fields only.) So I finally followed your advice and removed the special case. In the same time, I refactored a bit the code.

[jdemeyer]

Reviewer: Jeroen Demeyer

[vbraun]

Changed branch from u/bruno/polynomial_powering_positive_characteristic to 41dd905