fast PowerSeries_poly multiplication

[sagetrac-pernici]

In this patch truncated multiplication of dense polynomials is
used in PowerSeries_poly multiplication.

in Sage-4.6 on my computer with Intel Core i7 2.8GHz

sage: R.<a,b> = QQ[]
sage: K.<t> = PowerSeriesRing(R)
sage: time p1 = (1 + a*t + b*t^2 + O(t^50))^-40
Wall time: 7.62s

with this patch it takes 0.12s
The speed-up increases with the number of variables and with the precision.

Apply: trac_10480_fast_PowerSeries_poly_multiplication-alternative.patch

Depends on #10255

CC: @nilesjohnson @zimmermann6

Component: commutative algebra

Keywords: power series

Author: Mario Pernici, Luis Felipe Tabera Alonso

Branch/Commit: public/ticket/10480 @ e2cdc54

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

[sagetrac-pernici]

Attachment: trac_10480_fast_PowerSeries_poly_multiplication.patch.gz

[lftabera]

comment:1

Could you please check if the changes of #10255 in karatsuba multiplication provide further improvement?

PowerSeries are constructed over commutative rings, but polynomial_elements do not. Hence mul_trunc needs to work also over noncommutative rings. It does not right now (al least lines 5342 and 5345 do not differenciate between left and right multiplication).

I see several improvements:

I think that

cdef int i, j, a1len, a2len, n1, n2, ih

in line 5334 should be Py_ssize_t, but maybe a cython expert should look at this.

c = [0]*h in line 5346 of the patch is inefficient in several cases. One of the most contributions in 10480 is to avoid python 0 'int' in the code, that can be (relatively) expensive to coerce. I guess that this code will be slow if the base ring is a number field (that has slow coercion).

The last loop is classical multiplication up to the precision required. I wonder if part of that multiplication could use do_karatsuba or a truncated_do_karatsuba.

[lftabera]

comment:2

Some more numbers with your patch showing what I said on my previous post.
Patch aginst a clean sage 4.6

with patch:

sage: K.<t>=QQ['a'][]
sage: p1 = K.random_element(800) + O(t^990)
sage: p2 = K.random_element(800) + O(t^990)
sage: %time _=p1*p2
CPU times: user 6.09 s, sys: 0.08 s, total: 6.17 s
Wall time: 6.49 s

without patch

sage: sage: K.<t>=QQ['a'][]
sage: sage: p1 = K.random_element(800) + O(t^990)
sage: sage: p2 = K.random_element(800) + O(t^990)
sage: sage: %time _=p1*p2
CPU times: user 3.24 s, sys: 0.03 s, total: 3.27 s
Wall time: 3.39 s

I think that we can define an algorithm that works like karatsuba but that has a prec parameter. If a portion will have order greater than the prec. That part is discarded.
The prec is updated on recursive calls.

[lftabera]

comment:4

I will try to implement the algorithm at "A long note on Mulders’ short product" of Hanrot and Zimmermann that looks not harder than Karatsuba. The authors claim that heuristically the gain time is on the average 0.7 x Karatsuba time.

Also, the following cases have to be correctly covered:

sage: O(x)*O(x)
0
sage: O(x^2)*O(x^3)
0

[sagetrac-pernici]

Attachment: trac_10480_fast_PowerSeries_poly_multiplication2.patch.gz

Attachment: un.sage.gz

benchmark file

[sagetrac-pernici]

comment:5

Hello,

lftabera wrote:

I think that we can define an algorithm that works like karatsuba but that has a prec parameter. If a portion will have order greater than the prec. That part is discarded. The prec is updated on recursive calls.

I followed your suggestion to perform truncated multiplication recursively,
keeping track of the precision, calling Karatsuba multiplication
when there are no truncations.

Also, the following cases have to be correctly covered:

sage: O(x)*O(x)
0
sage: O(x^2)*O(x^3)
0

The new patch trac_10480_fast_PowerSeries_poly_multiplication2.patch passes all tests, included these.
It performs generally better than the unpached Sage4.6,
except the case in which the series has degree lower than the precision,
like your example

p1 = K.random_element(800) + O(t^990)

on which it is typically up to 50% slower, while with the previous
patch it was 2x slower.

The first patch was also slow in other cases, see below.

The performance of truncated multiplication vs Karatsuba multiplication
depends on the distribution of the terms in the univariate series;
truncated multiplication is fast when the high powers in the univariate
variable t have many more terms than the low powers, like in this example

using this patch:

sage: from sage.rings.power_series_poly import PowerSeries_poly
sage: R.<a,b> = QQ[]
sage: K.<t> = R[]
sage: p = (1 + a*t + b*t^2 + O(t^50))^-1
sage: v = [len(x.coefficients()) for x in p.coefficients()]
sage: v[:6], v[-6:]
([1, 1, 2, 2, 3, 3], [23, 23, 24, 24, 25, 25])
sage: %timeit p*p
25 loops, best of 3: 26 ms per loop
sage: %timeit PowerSeries_poly(p.parent(),p.polynomial().mul_trunc(p.polynomial(),50,1),50)
125 loops, best of 3: 4.89 ms per loop
sage: p*p - PowerSeries_poly(p.parent(),p.polynomial().mul_trunc(p.polynomial(),50,1),50)
O(t^50)

unpached Sage (main):

sage: %timeit p*p
5 loops, best of 3: 833 ms per loop

mul_trunc(self, right, h, generic=1) calls do_mul_trunc_generic, which
is the algorithm used in the first patch; in this case it is faster,
but in other cases it is slower, so I think it is better to use do_mul_trunc.
However in multivariate series do_mul_trunc_generic seems to perform always
better; I will post a patch about this in ticket #1956

p = K.random_element(50) + O(t^50) 

has uniform distribution; for p*p the patched version performs as the main version.

In un.sage one takes N subsequent p=p*(p+1) starting with a
random univariate series on a polynomial ring with n variables;
the ratio of the number of terms in the upper half of the univariate series
and the lower half is given.
The ratio is close to 1 in a random polynomial, and increases with n and N.
There is a correlation between the speedup
using truncated multiplication and this ratio.
The old patch (patch1) in some cases performs badly; the worst case is given.

n variables in polynomial ring on QQ; h precision;
N = number of times p = p*(p+1) is taken
n=1
h=20  N=1     2    3     4    5    6   7
ratio  0.86  0.98  1     1    1    1   1
main   0.005 0.008 0.022 0.09 0.41 3.5 53.5
patch  0.004 0.006 0.016 0.06 0.31 2.5 41.5
h=100
ratio  0.93  1     1     1    1
main   0.049 0.14  0.47  1.6  7.9
patch  0.043 0.10  0.31  1.1  5.5
h=500
ratio  0.96  1.0
main   1.5   20.9
patch1 2.0   49.1
patch  1.3   17.0

n=2
h=20
ratio  0.81  1.06  1.10  1.19 1.34
main   0.009 0.04  0.30  4.3  134
patch  0.007 0.03  0.22  2.6  61
h=100
ratio  1.05  1.01  1.0   1.0
main   0.11  0.80  7.8   171
patch  0.09  0.59  5.5   103

n=4
h=20
ratio  1.08  1.53  1.64
main   0.015 0.35  21.0
patch  0.010 0.19  10.9

n=8
h=10
ratio  1.28  2.67  5.43
main   0.008 0.38  127
patch  0.006 0.046 4.4

n=16
h=10
ratio  1.00  2.77  8.90
main   0.013 1.95  > 19m(6GB)
patch  0.009 0.14  29.0(0.5GB)

In the last example I stopped the computation with main because too much
memory was used.

I will try to implement the algorithm at "A long note on Mulders’ short product" of Hanrot and Zimmermann that looks not harder than Karatsuba. The authors claim that heuristically the gain time is on the average 0.7 x Karatsuba time.

If I understand correctly, that paper deals with univariate series
on numeric rings. Maybe in the case of univariate series over polynomial
rings one could use the list of the number of elements to guide the
algorithm.

[lftabera]

comment:6

Oops, we have been duplicating work,

I have written a proof of concept of a truncated Karatsuba to check that it derives a correct algorithm. It is much much slower than plain sage. It needs to eliminate unnecesary slicing and calls of type K(f).list()

I will take a look at your patch.

My patch depends on #10255