speed of multiplying polynomials

[dmha...@math.harvard.edu]

Hi,

I'm trying to figure out how fast mathematica can multiply polynomials
with integer coefficients. I don't know mathematica well at all, but I
had a friend write me a mathematica program to test this. It look like
on a regular desktop machine it can multiply e.g. degree 1000
polynomials with coefficients in the range 0 <= x <= 1000 in about 1.3
seconds.

This is ludicrously slow compared to some other computer algebra
systems, which can do this multiplication in about 0.0003
seconds. I can't believe mathematica is in the order of 10000 times
slower for this simple task. I think perhaps we are doing something
wrong. Can anyone suggest a way of coaxing mathematica into doing this
kind of arithmetic at a comparable pace?

Many thanks

David

[Roman Pearce]

The speed will depend heavily on the representation of the polynomials and
the algorithms used. If Mathematica is using a sparse representation
then any algorithm will be at least O(d^2). I know another system
that uses a uses a dense representation for univariate polynomials,
with Karatsuba and FFT multiplication. These methods are an order of
magnitude faster if your polynomials are dense.

[dimitris]

Don't judge Mathematica for a specific performance. It is the overal
that makes her better for other CAS. For example compare the obtained
results by Symbolic Integration of Definite Integrals and you will be
quite impressed about its skills.
Anyway, it just happens some CAS to be better on one task, some on
another.

To tell you the truth I believed (and I still believe somehow!)
Mathematica can achieve better performance here contrary to my
findings...
May be someone will suggest a better way!

In[64]:=
Plus @@ Table[Random[Integer, {0, 1000}]*x^i, {i, 0, 1000}]
Plus @@ Table[Random[Integer, {0, 1000}]*x^i, {i, 0, 1000}]

In[69]:=
Timing[Collect[%*%%, x]][[1]]
Out[69]=
2.3279999999999745*Second

(CPU 2.8 GHz)

[mcmc...@unca.edu]

On Jan 6, 11:28 pm, dmhar...@math.harvard.edu wrote:
> I'm trying to figure out how fast mathematica can multiply polynomials
> with integer coefficients.

Depends how you do it. Here's very simple way to
multiply polynomials.

poly := Sum[Random[Integer, {1, 1000}]*x^i, {i, 0, 1000}];
SeedRandom[1];
Expand[poly*poly]; // Timing

{1.20313 Second, Null}

On the other hand, suppose we represent a polynomial
as a list; then, we can use ListConvolve.

listpoly := Table[Random[Integer, {1, 1000}], {1001}];
SeedRandom[1];
ListConvolve[PadRight[listpoly, 2001],
PadRight[listpoly, 2001], 1]; // Timing

{0.002973 Second, Null}

About 400 times faster.

You can remove the semi-colons to examine the output and
see that the results are equivalent.

Lists are a fundamental data structure and in this case
we are able to express the desired object using a list and
we can express the operation using a command which acts
on Lists. Symbolic manipulation, by contrast, is likely
to be much slower.

Mark

[Carl Woll]

dmha...@math.harvard.edu wrote:

>Hi,

>
>I'm trying to figure out how fast mathematica can multiply polynomials

>with integer coefficients. I don't know mathematica well at all, but I
>had a friend write me a mathematica program to test this. It look like
>on a regular desktop machine it can multiply e.g. degree 1000
>polynomials with coefficients in the range 0 <= x <= 1000 in about 1.3
>seconds.
>
>This is ludicrously slow compared to some other computer algebra
>systems, which can do this multiplication in about 0.0003
>seconds. I can't believe mathematica is in the order of 10000 times
>slower for this simple task. I think perhaps we are doing something
>wrong. Can anyone suggest a way of coaxing mathematica into doing this
>kind of arithmetic at a comparable pace?
>
>Many thanks
>
>David
>
>

I think the problem is with the representation of a polynomial. Assuming
we are dealing with univariate polynomials, we can represent the
polynomial as a list of coefficients:

c = {900, 801, etc}

or as a polynomial with explicit multiplications and additions:

p = 900 + 801*x + etc.

If we work with lists, we can use ListConvolve to do the multiplication.
If we work with polynomials, we will probably use Expand. Here is an
example multiplying two polynomials in both ways. Here are the coefficients:

In[1]:=
c1 = Table[Random[Integer, {0, 1000}], {1001}];
c2 = Table[Random[Integer, {0, 1000}], {1001}];

Here are the equivalent polynomials:

In[3]:=
p1 = x^Range[0, 1000] . c1;
p2 = x^Range[0, 1000] . c2;

Here we time multiplication using ListConvolve:

In[5]:=
c = ListConvolve[c1, c2, {1, -1}, 0]; // Timing
Out[5]=
{0. Second, Null}

Here we time multiplication using Expand:

In[6]:=
p = Expand[p1 p2]; // Timing
Out[6]=
{2.25 Second, Null}

Finally, we check that the two approaches yield the same result:

In[7]:=
c === CoefficientList[p, x]
Out[7]=
True

So, if one is interested in multiplying large, dense univariate
polynomials in Mathematica, the fastest approach is to use ListConvolve.

Carl Woll
Wolfram Research

[Giovanni Resta]

dmha...@math.harvard.edu wrote:

> This is ludicrously slow compared to some other computer algebra
> systems, which can do this multiplication in about 0.0003
> seconds. I can't believe mathematica is in the order of 10000 times
> slower for this simple task. I think perhaps we are doing something
> wrong. Can anyone suggest a way of coaxing mathematica into doing this
> kind of arithmetic at a comparable pace?

I think that the point is that Mathematica multiplies the
two expression without regarding them as polynomials, but as general
expressions.
Indeed, if I have two polynomials:

a = Sum[Random[Integer,{1,1000}] x^k,{k,0,1000}];
b = Sum[Random[Integer,{1,1000}] x^k,{k,0,1000}];

And I compute the explicit product:

Timing[c = Expand[a b];]
{1.28408 Second, Null}

If I have something that is not a polynomial:
Timing[c = Expand[(a+7/x)(b+Sin[x])];]
{1.36808 Second, Null}

I got a time which is very near.
Hence, probably Mathematica is not taking advantage of the fact
that a and b are polynomials, in the first case.
In other symbolic programs you have to declare explicitly that something
is a polynomial (sometimes you can manipulate only polynomials...) so
easily the computation is faster.
I do not know if there is a way to tell mathematica that that
expressions are polynomials. Usually to perform such a simple
tasks, if I need speed, I write a specific C program
(unless coefficients are symbolic...).
g.

[Paul Abbott]

In article <ent4s4$etq$1...@smc.vnet.net>, Carl Woll <ca...@wolfram.com>
wrote:

All this is true but the essential question -- one that David Harvey
alluded to in his original post -- is why does Mathematica not do this
automatically? For example, multiplication of large integers and
high-precision approximate numbers is done using interleaved schoolbook,
Karatsuba, three-way Toom-Cook and number-theoretic transform algorithms.

In other words if a computation involves multiplying large, dense
univariate polynomials then _internally_ Mathematica should use
ListConvolve. Now, given two polynomials, it takes time to check that
each _is_ a polynomial,

PolynomialQ[p1, x]; // Timing
PolynomialQ[p2, x]; // Timing

time to extract the coefficients,

c1 = CoefficientList[p1, x]; // Timing
c2 = CoefficientList[p2, x]; // Timing

time to multiplication using ListConvolve:

c = ListConvolve[c1, c2, {1, -1}, 0]; // Timing

and time to reconstruct the resulting polynomial

p = x^Range[0, Length[c] - 1] . c; // Timing

The total time of these operations is still less than that of direct
multiplication, but there are now significant time penalties -- which
makes one wonder why there is not an explicit Polynomial type or
representation involving the coefficients and variables (perhaps
invoking other tools such as sparse and/or packed arrays, as required
for sparse polynmomials) so that multiplication and other polynomial
operations is as fast as possible automatically?

Cheers,
Paul

_______________________________________________________________________
Paul Abbott Phone: 61 8 6488 2734
School of Physics, M013 Fax: +61 8 6488 1014
The University of Western Australia (CRICOS Provider No 00126G)
AUSTRALIA http://physics.uwa.edu.au/~paul

[mcmc...@unca.edu]

On Jan 12, 5:16 am, Paul Abbott <p...@physics.uwa.edu.au> wrote:

> ... why does Mathematica not do this

> automatically? For example, multiplication of large integers and
> high-precision approximate numbers is done using interleaved schoolbook,
> Karatsuba, three-way Toom-Cook and number-theoretic transform algorithms.

I'm not sure that your comparison to integer multiplication
(perhaps as fundamental an operation on as fundamental a
data type as one can imagine) is a fair one. Your code
suggests that we simply test if the input expressions are
polynomials in a variable x; this is already far more
complicated than checking for an integer. Furthermore, we
really need to know if the expression is a polynomial in
*some* variable. This seems quite a bit harder,
particularly if the expression does not look like
a polynomial.

What if the polynomial is not pre-expanded?
For example:

CoefficientList[(x+a)(x-b),x]
{-a b,a-b,1}

This has serious consequences for the time complexity of
of an algorithm using List Convolve, even if a and b are
integers.

More fundamentally, Expand and Collect are primarily
symbolic operations while ListConvolve is primarily a
numerical operation; each command works broadly, but is
strongest in its own domain. I think it is reasonable,
even necessary, to expect users to understand this - in
much the same way that users need to understand when to
use Solve vs NSolve. Although, this is an admittedly
subtle situation.

Here's an example to illustrate the time complexity in
various situation. First, define a PolynomialMultiply
function based on ListConvolve:

PolynomialMultiply[Times[p1_, p2_], x_] := Module[
{c1, c2, c},

c1 = CoefficientList[p1, x];

c2 = CoefficientList[p2, x];

c = ListConvolve[c1, c2, {1, -1}, 0];

x^Range[0, Length[c] - 1].c] /;
PolynomialQ[p1, x] && PolynomialQ[p2, x];

It is indeed quite fast when multiplying pre-expanded
polynomials with integer coefficients.

poly := Sum[Random[Integer, {1, 1000}]*x^i, {i, 0, 1000}];
SeedRandom[1];

PolynomialMultiply[poly*poly]; // Timing
{0.009521 Second, Null}

vs. Expand

SeedRandom[1];
Expand[poly*poly];//Timing
{1.21226 Second,Null}

But the Expand version works much faster with symbolic
coefficients - even very simple ones

p1 = Sum[a[n]x^n,{n,0,1000}];
p2 = Sum[b[n]x^n,{n,0,1000}];
PolynomialMultiply[p1*p2,x];//Timing
{45.9369 Second,Null}

Expand[p1*p2];//Timing
{1.68533 Second,Null}

Another observation: the different commands may return
answers in different forms. For example:

p1 = (a+b)x+c;
p2 = d*x+(e+f);

Expand[p1*p2] //InputForm
c*e + c*f + c*d*x + a*e*x + b*e*x + a*f*x + b*f*x + a*d*x^2 + b*d*x^2

Expand[p1*p2,x] //InputForm
c*e + c*f + c*d*x + (a + b)*e*x + (a + b)*f*x + (a + b)*d*x^2

Collect[p1*p2,x] //InputForm
c*e + c*f + (c*d + (a + b)*e + (a + b)*f)*x + (a + b)*d*x^2

PolynomialMultiply[p1*p2,x] //InputForm
c*(e + f) + (c*d + (a + b)*(e + f))*x + (a + b)*d*x^2

Which is the correct answer?

MM

[dmha...@math.harvard.edu]

Paul Abbott wrote:
> All this is true but the essential question -- one that David Harvey

> alluded to in his original post -- is why does Mathematica not do this

> automatically? For example, multiplication of large integers and
> high-precision approximate numbers is done using interleaved schoolbook,
> Karatsuba, three-way Toom-Cook and number-theoretic transform algorithms.

Exactly, thanks Paul.

Actually, there is another issue. I have tried ListConvolve for a few
cases, and compared its performance to polynomial multiplication in
certain other computer algebra systems.

For large degree, small coefficient problems (say degree 10000,
coefficients in [0, 1000]), Mathematica is definitely in the same ball
park, maybe still a little slower, but not hugely slower.

But for small degree, larger coefficients (say degree 100, coefficients
in [0, 10^500]), Mathematica was something like 20 times slower.

Is ListConvolve really the best Mathematica can do? I mean, if someone
really desperately needs to multiply polynomials in Mathematica, I find
it difficult to believe that the best they can do is still twenty times
slower than other systems. Of course they could go out and write the
multiplication code in C or something themselves, but that's missing
the point :-)

David

[Paul Abbott]

In article <eoa6fh$mfs$1...@smc.vnet.net>, mcmc...@unca.edu wrote:

> On Jan 12, 5:16 am, Paul Abbott <p...@physics.uwa.edu.au> wrote:
>
> > ... why does Mathematica not do this
> > automatically? For example, multiplication of large integers and
> > high-precision approximate numbers is done using interleaved schoolbook,
> > Karatsuba, three-way Toom-Cook and number-theoretic transform algorithms.
>
> I'm not sure that your comparison to integer multiplication
> (perhaps as fundamental an operation on as fundamental a
> data type as one can imagine) is a fair one. Your code
> suggests that we simply test if the input expressions are
> polynomials in a variable x; this is already far more
> complicated than checking for an integer. Furthermore, we
> really need to know if the expression is a polynomial in
> *some* variable. This seems quite a bit harder,
> particularly if the expression does not look like
> a polynomial.

I went on to write that:

> The total time of these operations is still less than that of direct
> multiplication, but there are now significant time penalties -- which
> makes one wonder why there is not an explicit Polynomial type or
> representation involving the coefficients and variables (perhaps
> invoking other tools such as sparse and/or packed arrays, as required
> for sparse polynmomials) so that multiplication and other polynomial
> operations is as fast as possible automatically?

Here I was asking why there is not an explicit Polynomial datatype. If
there was such a type, then one could either enter expressions as
polynomials -- and gain speed for certain operations -- or convert to
this type if required.

> More fundamentally, Expand and Collect are primarily
> symbolic operations while ListConvolve is primarily a
> numerical operation; each command works broadly, but is
> strongest in its own domain.

In what sense is ListConvolve a numerical operation? Do you mean it is
optimised for numerical arguments? Actually, as your code below shows,
it works fine with symbolic coefficients.

The majority of the time is taken testing whether p1 and p2 are
polynomials in x. Try

Timing[PolynomialQ[p1, x] && PolynomialQ[p2, x]]

to see this. If there was an explicit Polynomial datatype this test
would not be required, or at least would be _much_ faster. Indeed,
without this test, for the examples that I've tried, PolynomialMultiply
_beats_ Expand even for symbolic coefficients!

> Another observation: the different commands may return
> answers in different forms. For example:
>
> p1 = (a+b)x+c;
> p2 = d*x+(e+f);
>
> Expand[p1*p2] //InputForm
> c*e + c*f + c*d*x + a*e*x + b*e*x + a*f*x + b*f*x + a*d*x^2 + b*d*x^2
>
> Expand[p1*p2,x] //InputForm
> c*e + c*f + c*d*x + (a + b)*e*x + (a + b)*f*x + (a + b)*d*x^2
>
> Collect[p1*p2,x] //InputForm
> c*e + c*f + (c*d + (a + b)*e + (a + b)*f)*x + (a + b)*d*x^2
>
> PolynomialMultiply[p1*p2,x] //InputForm
> c*(e + f) + (c*d + (a + b)*(e + f))*x + (a + b)*d*x^2
>
> Which is the correct answer?

Perhaps the best answer is the list of coefficients?

CoefficientList[p1 p2, x] // Expand
{c*e + c*f, c*d + a*e + b*e + a*f + b*f, a*d + b*d}

[Roman Pearce]

> Actually, there is another issue. I have tried ListConvolve for a few
> cases, and compared its performance to polynomial multiplication in
> certain other computer algebra systems.

...

> But for small degree, larger coefficients (say degree 100, coefficients
> in [0, 10^500]), Mathematica was something like 20 times slower.

I bet it is incrementally adding up large numbers. For example, if I
have 1000 numbers c[0], ..., c[999], and those numbers are big, then
the obvious algorithm (in C):

for(i=0, sum=0; i < 1000; i++)
sum += c[i];

is a disaster. Each addition will be linear time in the length of the
partial sum, which is O(n^2). Adding up n numbers with n digits this
way is O(n^3), whereas a divide-and-conquer approach would be O(n^2).
I have no idea if Mathematica in fact does this, but the mistake is
surprisingly common. People do this with polynomials as well as
integers, and then they wonder why their algorithms are slow. It
doesn't matter if it's coded in C, this happens whenever the object
being constructed is large (ie: more than a machine word).

[Roman Pearce]

Roman Pearce wrote:
> I bet it is incrementally adding up large numbers. For example, if I
> have 1000 numbers c[0], ..., c[999], and those numbers are big, then
> the obvious algorithm (in C):
>
> for(i=0, sum=0; i < 1000; i++)
> sum += c[i];
>
> is a disaster.

Sorry, this is nonsense. I was confusing addition with multiplication.