Correct karatsuba multiplication of univariate polynomials for different degree polynomials

[lftabera]

In the generic implementation of univariate polynomials over exact rings, plain karatsuba is used.

However, for different degree polynomials, the karatsuba code makes more products and additions than the classical multiplication code. Moreover, for equally sized polynomials, the degree in which karatsuba starts to be better than plain multiplication looks too high for me.

See attachments comparison_product_50_400.png and comparison_addition_50_400.png for the number of operations multiplying degree 50 and 400 polynomials, in yellow, the number of operations using _mul_generic, in red using current _mul_karatsuba as of 4.6.1 and in blue with the patch proposed.

sage: K=QQ[I][x]   
sage: f=K.random_element(1500)
sage: g=K.random_element(1500)
sage: %time _ = f._mul_generic(g)
CPU times: user 6.03 s, sys: 0.08 s, total: 6.11 s
Wall time: 6.43 s
sage: %time _ = f._mul_karatsuba(g)
CPU times: user 6.95 s, sys: 0.06 s, total: 7.01 s
Wall time: 7.76 s

See comment:13 for some benchmarks

Apply: trac_10255_karatsuba_improbement.patch

Component: basic arithmetic

Keywords: karatsuba, multiplication, polynomial

Author: Luis Felipe Tabera Alonso, Marc Mezzarobba

Branch/Commit: u/mmezzarobba/10255-karatsuba @ 5977470

Reviewer: Marc Mezzarobba, Luis Felipe Tabera Alonso

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

[lftabera]

comment:1

Same example with my patch. Pure karatsuba:

sage: K=QQ[I][x]   
sage: f=K.random_element(1500)
sage: g=K.random_element(1500)
sage: %time _ = f._mul_karatsuba(g)
CPU times: user 1.53 s, sys: 0.01 s, total: 1.54 s
Wall time: 1.73 s

TODO:

-Take a more careful look to the correctness of code.
-Check that the changes do not introduce performance penalties on a wider set of rings, degrees and size of polynomials.
-Introduce a user-definable threshold (with a sensible default for most common rings) for which karatsuba falls back to classic multiplication.

[lftabera]

comment:2

the code is too naive for different sized polynomials. Even with my improvements I get this:

Two degree 849 polynomials, no threshold:

274750 additions

87729 products

One degree 849 polynomial and a 449 polynomial:

472812 additions

218649 products

Experimental data

sage: f=ZZ[x].random_element(849)
sage: g=ZZ[x].random_element(849)
sage: %timeit f._mul_karatsuba(g)
5 loops, best of 3: 95.9 ms per loop
sage: f=ZZ[x].random_element(849)
sage: g=ZZ[x].random_element(449)
sage: %timeit f._mul_karatsuba(g)
5 loops, best of 3: 218 ms per loop

The shorter product is more expensive!

I have a patch right now to solve this, have to clean it up though.

[lftabera]

comment:3

Further improvements for polynomials of of different sizes:

if f is of degree m-1 and g is of degree n-1 with m<n, consider g as:

g0 + x^m*g1 + x^(2*m)g2 + ... + x^(q*m)*gq

Compute the products fgi with karatsuba and reconstruct fg

With this strategy, for polynomials of degree 849 and 449 I get:

214507 additions

39942 products

Far better from previous method. Anyway this is not optimal and it is hard to messure the recursion overhead. I have taken this strategy from old versions of gmp documentation.

Some timmings (Debian 64bits). All polyomials have been generated with random_element without further parameters. Pure karatsuba without threshold for better comparison with old code. With an appropriate threshold the timming whould slightly better.

Note that for rigns like QQ[x] and GF(2)[x] sage uses faster special libraries (FLINT, NTL and the like).

Degree f:1000,g:1000

QQ[x] 
old code: 828 ms
new code: 498 ms
GF(2)[x] 
old code: 284 ms
new code: 99.1 ms
QQ[sqrt(2)][x]
old code: 2.61 s
new code: 480 ms # 5 times faster
ZZ[x]['y']['z']['t']
old code: 86.95 s
new code: 71.10 s

degree f:849, g:449

QQ[x] 
old code: 988 ms
new code: 343 ms
GF(2)[x] 
old code: 148 ms
new code: 67.6 ms
QQ[sqrt(2)][x]
old code: 1.93 s
new code: 313 ms #6 times faster
ZZ[x]['y']['z']['t']
old code: 342.16 s
new code: 48.11 s # 7 times faster

degree 3, 5

QQ[x] 
old code: 135 µs
new code: 83.6 µs
GF(2)[x] 
old code: 137 µs
new code: 100 µs
QQ[sqrt(2)][x]
old code: 323 µs
new code: 86.4 µs
ZZ[x]['y']['z']['t']
old code: 19.5 ms
new code: 13.2 ms 

TODO:

• Sensible default threshold.

• Documentation cleanup. (almost done)

• More testing for correctness rather than speed. This is a very basic feature, so we cannot allow correctness bugs here.

[lftabera]

comment:4

All doctest passes and I am trying to make doctest independent of future changes of the default threshold. But I get the following difference that puzzles me.

sage -t -long -force_lib "devel/sage/sage/rings/qqbar.py"   
**********************************************************************
File "/home/rosamari/sage/devel/sage/sage/rings/qqbar.py", line 373:
    sage: sage_input(one, verify=True)
Expected:
    # Verified
    R.<x> = QQbar[]
    v1 = AA(2)
    v2 = QQbar(sqrt(v1))
    v3 = QQbar(3)
    v4 = sqrt(v3)
    v5 = v2*v4
    v6 = (1 - v2)*(1 - v4) - 1 - v5
    v7 = QQbar(sqrt(v1))
    v8 = sqrt(v3)
    si1 = v7*v8
    cp = AA.common_polynomial(x^2 + ((1 - v7)*(1 + v8) - 1 + si1)*x - si1)
    v9 = QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
    v10 = 1 - v9
    v11 = v6 + (v10 - 1)
    v12 = -1 - v4 - QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
    v13 = 1 + v12
    v14 = v10*(v6 + v5) - (v6 - v5*v9)
    si2 = v5*v9
    AA.polynomial_root(AA.common_polynomial(x^4 + (v11 + (v13 - 1))*x^3 + (v14 + (v13*v11 - v11))*x^2 + (v13*(v14 - si2) - (v14 - si2*v12))*x - si2*v12), RIF(RR(0.99999999999999989), RR(1.0000000000000002)))
Got:
    # Verified
    R.<x> = QQbar[]
    v1 = QQbar(3)
    v2 = sqrt(v1)
    v3 = AA(2)
    v4 = QQbar(sqrt(v3))
    v5 = sqrt(v1)
    si1 = v4*v5
    cp = AA.common_polynomial(x^2 + ((1 - v4)*(1 + v5) - 1 + si1)*x - si1)
    v6 = -1 - v2 - QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
    v7 = QQbar.polynomial_root(cp, RIF(-RR(1.7320508075688774), -RR(1.7320508075688772)))
    v8 = QQbar(sqrt(v3))
    v9 = v8*v2
    v10 = (1 - v8)*(1 - v2) - 1 - v9
    v11 = -v7 + v10
    v12 = v6*v11
    si2 = v7*v9
    v13 = (1 - v7)*(v10 + v9) - v10 + si2
    si3 = v6*si2
    AA.polynomial_root(AA.common_polynomial(x^4 + ((1 + v6)*(1 + v11) - 1 - v12)*x^3 + (v12 + v13)*x^2 + ((1 + v6)*(v13 - si2) - v13 + si3)*x - si3), RIF(RR(0.99999999999999989), RR(1.0000000000000002)))
**********************************************************************
1 items had failures:
   1 of 152 in __main__.example_0
***Test Failed*** 1 failures.

[lftabera]

comment:5

I have eliminated karatsuba_sumand karatsuba_diff so there are less function calls and less if branches in the coding (we know an can predict at each step the length of every list.

I have also eliminated most slicing of lists so less objects are creating. The recursion does not operates over slices but over the orginial list plus some pointers to the slice we are working on.

With this improvements the code is about 10% faster. I also think that improves readability.

For number fields of degree 2 and 3 the code is about three times slower than pari. For number fields of degree >=6 sage seems faster than pari.

For dense bivariate polynomials in QQ[x][y] the multiplication starts to be faster than in QQ[x,y] for degree 16, for dense polynomials of degree 100 QQ[x][y] is 45x faster than QQ[x,y]. Of course, for non-dense polynomials multivariate multiplication is much much faster.

Some doctests added to compare karatsuba with generic multiplication and doctest for the non-commutative ring case.

Still has the QQbar fail.

[sagetrac-spancratz]

comment:6

Hello,

First of all, I want to say that it's great that you are looking at this part of the Sage library and are trying to improve it!

From reading your ticket progress, it seems that perhaps you might benefit from some help with this. Have you considered advertising this problem on the developers' list at sage-devel? Perhaps someone else is interested in this, too.

In any case, here are a couple of comments:

1)  What is your main motivation for wanting to improve this code?  If you want to achieve faster computations in polynomial rings over specific rings, QQ[i], you should most definitely either write some Cython code or even a separate C library.  You should *not* be tinkering with the generic Python multiplication code.  To keep with the example, if you are interested in QQ[i], then (just as a first idea) you should represent your polynomials as sums of two polynomials internally, one for the real part and one for the imaginary code.  This approach also applies to other number fields of small degree.

2)  Stating what you definitely already know:  Karatsuba multiplication trades additions and multiplications in the base ring.  Thus, the threshold highly depends on the relative speeds of addition and multiplication, and it will be very difficult for you to make a default choice that is useful in general.  Even worse, what about base rings where multiplication is faster than addition?

3)  I haven't looked at the details, but I believe you might simply not be able to apply this code to polynomial rings with inexact base rings.

All that said, while this patch looks very interesting, I think there are a lot of things that need to be sorted out. By the way, I am not quite sure whether any one else is currently reading this ticket, but at least during the next week I should have a lot of time to look at this.

Best wishes,

Sebastian

[lftabera]

comment:7

Sebastian,

Thanks for looking at this. I feel that the ticket is ready for review, but there is a doctest fail and I wanted to understand why. My impression is that the doctest should be changed, but I did not have enough time to dig it.

Replying to @sagetrac-spancratz:

1)  What is your main motivation for wanting to improve this code?  If you want to achieve faster computations in polynomial rings over specific rings, QQ[i], you should most definitely either write some Cython code or even a separate C library.  You should *not* be tinkering with the generic Python multiplication code.  To keep with the example, if you are interested in QQ[i], then (just as a first idea) you should represent your polynomials as sums of two polynomials internally, one for the real part and one for the imaginary code.  This approach also applies to other number fields of small degree.

My primary interest was univariate polynomial rings over number fields of degree around 30-40. I was considering to use the specific class I have added in #8558, because, for number field base rings, the slowness comes really from adding number field elements by python 0. But then I realised that the karatsuba multiplication is not well enough implemented if both polynomials have different degree, so I changed the cython code improving here and there.

2)  Stating what you definitely already know:  Karatsuba multiplication trades additions and multiplications in the base ring.  Thus, the threshold highly depends on the relative speeds of addition and multiplication, and it will be very difficult for you to make a default choice that is useful in general.  Even worse, what about base rings where multiplication is faster than addition?

That is why the default threshold is zero. No classical multiplication at all. But I have added a ring-dependent parameter to change this value. Univariate polynomial rings have now a method set_karatsuba_threshold to deal with this issue. Rings in which addition is faster than multiplication also benefits from karatsuba, since the number of additions is also subquadratic in the degree of the polynomials.

3)  I haven't looked at the details, but I believe you might simply not be able to apply this code to polynomial rings with inexact base rings.

By default, unless the user calls explicitelly _mul_karatsuba, mutiplication for inexact rings is the classical multiplication algorithm. I have not changed this. Also,if a class (ZZ[x], GF(5)[x], etc.) has its own multiplication code, my patch does not change this.

All that said, while this patch looks very interesting, I think there are a lot of things that need to be sorted out. By the way, I am not quite sure whether any one else is currently reading this ticket, but at least during the next week I should have a lot of time to look at this.

Thanks!, if you have any question on why I changed this or that or have further suggestions/improvements, do not hesitate to ask.

Luis

[jdemeyer]

comment:8

One can actually a speedup of more than a factor 10 by using PARI:

The following function converts a polynomial over a number field to a PARI object:

def nfpol_to_pari(f):
    return pari([c._pari_('a') for c in f.list()]).Polrev()

Now let's try your example:

sage: K=QQ[I][x]   
sage: f=K.random_element(1500)
sage: g=K.random_element(1500)
sage: %time _ = f._mul_generic(g)
CPU times: user 3.57 s, sys: 0.00 s, total: 3.57 s
Wall time: 3.59 s
sage: %time _ = f._mul_karatsuba(g)
CPU times: user 5.06 s, sys: 0.05 s, total: 5.11 s
Wall time: 5.12 s
sage: fpari = nfpol_to_pari(f)
sage: gpari = nfpol_to_pari(g)
sage: %time _ = fpari*gpari
CPU times: user 0.26 s, sys: 0.00 s, total: 0.26 s
Wall time: 0.26 s

Since PARI is written completely in C, it's obvious why one could expect such a speedup. Obviously, writing specialized code to deal with number field polynomials would even be better.

[lftabera]

comment:9

Obviously you did not try the computations with my patch, it is still much slower than PARI for small number fields (as I acknowledge in the ticket, PARI is around three times faster for degree two extensions):

sage: K=QQ[I][x]
sage: f=K.random_element(1500)
sage: g=K.random_element(1500)
sage: %time _ = f._mul_generic(g)
CPU times: user 5.94 s, sys: 0.02 s, total: 5.96 s
Wall time: 6.01 s
sage: %time _ = f._mul_karatsuba(g)
CPU times: user 1.31 s, sys: 0.00 s, total: 1.32 s
Wall time: 1.33 s
sage: K.set_karatsuba_threshold(8)
sage: %time _ = f._mul_karatsuba(g)  
CPU times: user 0.97 s, sys: 0.00 s, total: 0.97 s
Wall time: 1.00 s
sage: %time fpari = nfpol_to_pari(f)
CPU times: user 0.24 s, sys: 0.00 s, total: 0.24 s
Wall time: 0.25 s
sage: %time gpari = nfpol_to_pari(g)
CPU times: user 0.24 s, sys: 0.00 s, total: 0.24 s
Wall time: 0.25 s
sage: %time hpari = fpari*gpari
CPU times: user 0.42 s, sys: 0.00 s, total: 0.42 s
Wall time: 0.44 s
sage: %time h = K(hpari)       
CPU times: user 0.70 s, sys: 0.02 s, total: 0.72 s
Wall time: 0.73 s

If you take into account a naive transformation from sage data to pari and back, the total time is worse than the mul_karatsuba time. So if a specialiczed class would use PARI, this transformation has to be taken into account.

This pari advantage disapears for bigger number fields

sage: K=CyclotomicField(7)[x]
sage: f=K.random_element(1500)
sage: g=K.random_element(1500)
sage: %time _ = f._mul_karatsuba(g) 
CPU times: user 15.33 s, sys: 0.02 s, total: 15.35 s
Wall time: 15.44 s
sage: %time _ = f._mul_karatsuba(g,8)
CPU times: user 12.96 s, sys: 0.08 s, total: 13.03 s
Wall time: 13.10 s
sage: %time fpari = nfpol_to_pari(f)
CPU times: user 0.30 s, sys: 0.00 s, total: 0.30 s
Wall time: 0.32 s
sage: %time gpari = nfpol_to_pari(g) 
CPU times: user 0.30 s, sys: 0.00 s, total: 0.31 s
Wall time: 0.32 s
sage: %time hpari = fpari*gpari      
CPU times: user 24.26 s, sys: 0.02 s, total: 24.28 s
Wall time: 25.23 s
sage: %time h = K(hpari)             
CPU times: user 1.48 s, sys: 0.03 s, total: 1.50 s
Wall time: 1.52 s

Note: for absolute number fields, I am used to set a karatsuba threshold in 8.

[sagetrac-spancratz]

comment:10

Hi Luis,

I'm sorry for the delay in my reply. Also, thank you Jeroen for adding the explicit timings. Let me just go through a couple of things:

• "I feel that the ticket is ready for review"

I don't really agree. It is certainly in a rather useful state, but it requires more people looking at it, using it, and providing feedback.

• For one, as you say, this "is a very basic feature, so we cannot allow correctness bugs here." There should be a lot more tests. Rather than having a single test case for illustrative purposes (as is the case in many places in the Sage library), there should be more bulk tests here. You could have that in the Sage framework e.g. by creating a list with f._mul_karatsuba(g) - f*g for many random pairs (f,g), then remove duplicates from the list, and check that the list has length one and only contains 0.

• Your code is still noticeably slower than PARI for the base ring being a number field of small degree. I guess strongly that this will be a much more typical use case than large degree number fields and as such this is not acceptable.

• Another one of your examples is ZZ[x]['y']['z']['t']. I think this is a very bad example. I don't think one should ever expect a user to work with a multivariate polynomial ring in this way just to get some obscure speed-up. Of course, dense and sparse multivariate polynomial rings offer massively different performance characteristic. Perhaps there should be a dense option to MPolynomialRing, but that's not for this ticket. In any case, this example should not be relevant here. Finally, if you really care about dense multivariate polynomials over the integers, I am sure this can be done much, much faster using C or Cython and arrays of mpz_t's for the coefficients. The speed-up of this will be much, much larger than anything you get from nesting a bunch univariate polynomial rings in Python.

• Finally, the statement

"If you take into account a naive transformation from sage data to pari and back, the total time is worse than the mul_karatsuba time. So if a specialiczed class would use PARI, this transformation has to be taken into account."

is not true. In short, the code for polynomials over number fields should then be changed so that the internal representation simply was the PARI object --- no converting back and forth in between the computations whatsoever!

I have looked at the currently only remaining failing test case in QQbar. Unfortunately, I am very puzzled by what could possibly cause this and don't have any help to offer right now.

By the way, I am sorry if this might sound a little harsh in places. Of course, as said, I am still very happy to keep looking at this etc.

Best wishes,

Sebastian