# [Numpy-discussion] [ANN]: Taylorpoly, an implementation of vectorized Taylor polynomial operations and request for opinions

Friedrich Romstedt friedrichromstedt@gmail....
Sat Feb 27 15:02:03 CST 2010

To the core developers (of numpy.polynomial e.g.): Skip the mess and

The other things I will post back to the list, where they belong to.
I just didn't want to have off-topic discussion there.

> I wanted to stress that one can do arithmetic on Taylor polynomials in
> a very similar was as with complex numbers.

I do not understand completely.  What is the analogy with complex
numbers which one cannot draw to real numbers?  Or, more precise: What
/is/ actually the analogy besides that there are operations?  With i,
i * i = -1, thus one has not to discard terms, contrary to the
polynomial product as you defined it, no?

> I guess there are also situations when you have polynomials   z =
> \sum_{d=0}^{D-1} z_d t^d, where z_d are complex numbers.

>> I like it more to implement operators as overloads of the __mul__
>
> I thought about something like
> def __mul__(self, other):
>    return mul(self,other)

Yes, I know!  And in fact, it may symmetrise the thing a bit.

>> etc., but this is a matter of taste.
>>  In fact, you /have/ to provide
>> external binary operators, because I guess you also want to have
>> numpy.ndarrays as left operand.  In that case, the undarray will have
>> higher precedence, and will treat your data structure as a scalar,
>> applying it to all the undarray's elements.
>
> well, actually it should treat it as a scalar since the Taylor
> polynomial is something like a real or complex number.

Maybe I misunderstood you completely, but didn't you want to implement
arrays of polynomials using numpy?  So I guess you want to apply a
vector from numpy pairwise to the polynomials in the P-object?

> [z]_{D+1} &=_{D+1}& [x]_{D+1} y_{D+1} \\
> \sum_{d=0}^D z_d t^d & =_{D+1} & \left( \sum_{d=0}^D x_d t^d \right)
> \left( \sum_{c=0}^D y_c t^c \right)

Did your forget the [] around the y in line 1, or is this intentional?
Actually, I had to compile the things before I could imagine what you
could mean.

Why don't you use multidimensional arrays?  Has it reasons in the C
accessibility?  Now, as I see it, you implement your strides manually.
With a multidimensional array, you could even create arrays of shape
(10, 12) of D-polynomials by storing an ndarray of shape (10, 12, D)
or (D, 10, 12).

Just because of curiosity:  Why do you set X_{0, 1} = ... X_{0, P} ?

Furthermore, I believe there is some elegant way formulating the
product for a single polynomial.  Let me think a bit ...

For one entry of [z]_E, you want to sum up all pairs:
x_{0} y{E} + ... + x{D} y{E - D} ,     (1)
right?  In a matrix, containing all mutual products x_{i} y{j}, this
are diagonals.  Rotating the matrix by 45° counterclockwise, they are
sums in columns.  Hmmm... how to rotate a matrix by 45°?

Another fresh look: (1) looks pretty much like the discrete
convolution of x_{i} and y_{i} at argument E.  Going into Fourier
space, this becomes a product.  x_i and y_i have finite support,
because one sets x_{i} and y_{i} = 0 outside 0 <= i <= D.  The support
of (1) in the running index is at most [0, D].  The support in E is at
most [0, 2 D].  Thus you don't make mistakes by using DFT, when you
pad x and y by D zeros on the right.  My strategy is: Push the padded
versions of x and y into Fourier space, calculate their product,
transform back, cut away the D last entries, and you're done!

I guess you mainly want to parallelise the calculation instead of
improving the singleton calculation itself?  So then, even non-FFT
would incorporate 2 D explicit python sums for Fourier transform, and
again 2 D sums for transforming back, is this an improvement over the
method you are using currently?  And you could code it in pure Python
:-)

I will investigate this tonight.  I'm curious.  Irrespective of that I
have also other fish to fry ... ;-)

iirc, in fact DFT is used for efficient polynomial multiplication.
Maybe Chuck or another core developer of numpy can tell whether
numpy.polynomial does it actually that way?

Friedrich