Prev: Minimization with constraint expresses as CDF[] >=
Next: pursuit curve (differential equations)
From: Giovanni Resta on 9 Jan 2007 07:01
dmharvey(a)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.

From: Paul Abbott on 12 Jan 2007 05:16
In article <ent4s4$etq$1(a)smc.vnet.net>, Carl Woll <carlw(a)wolfram.com>
wrote:

> dmharvey(a)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.

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

From: mcmcclure on 13 Jan 2007 03:49
On Jan 12, 5:16 am, Paul Abbott <p...(a)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

From: dmharvey on 13 Jan 2007 04:06

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

From: Paul Abbott on 16 Jan 2007 02:16
In article <eoa6fh$mfs$1(a)smc.vnet.net>, mcmcclure(a)unca.edu wrote:

> On Jan 12, 5:16 am, Paul Abbott <p...(a)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.

> 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}

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}

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

First  |  Prev  |  Next  |  Last
Pages: 1 2 3
Prev: Minimization with constraint expresses as CDF[] >=
Next: pursuit curve (differential equations)