Algorithm to compute the n-th derivative of a polynomial in Python

Question

As a personal exercise, I'm trying to write an algorithm to compute the n-th derivative of an ordered, simplified polynomial (i.e., all like terms have been combined). The polynomial is passed as an ordered list where the i-th index corresponds (though is not equivalent) to the coefficient of x to the n-th power.

Example: Take the derivative of: \3ドルx^3 + 5x^2 + 2x + 2\$ -> [3,5,2,2]

• 1st derivative: \9ドルx^2 + 10x + 2\$
• 2nd derivative: \18ドルx + 10\$
• 3rd derivative: \18ドル\$
• 4th...n-th derivative: \0ドル\$

Implementation in Python:

def poly_dev(poly, dev): 
 """
 :param poly: a polynomial
 :param dev: derivative desired, e.g., 2nd
 """
 c = 0 # correction
 r = dev # to remove
 while dev > 0:
 for i in range(1, len(poly)-c):
 poly[i-1] = (len(poly)-(i+c))*poly[i-1]
 dev -= 1 # I suspect this
 c += 1 # this can be simplified
 return poly[:-r]

E.g., print(poly_dev(poly = [3,5,2,2], dev = 2))

I have a math background, but I'm only starting to learn about computer science concepts like Complexity Theory.

I've intentionally tried to avoid reversing a list, as I know that can be expensive. What other steps can I change to decrease this procedure's run time?

• Gareth Rees: Clearly the polynomial should be stored the other way round.

Answer 1

You are going one by one derivative, when you could easily do them all at once. With a little multiply-and-divide trick, you could save on complexity.

Also, your function will return wrong results for negative dev (it should raise an exception instead).

Further, your function will mess up the original list:

>>> poly = [3, 5, 2, 2]
>>> print(poly_dev(poly, 2))
[18, 10]
>>> print(poly)
[18, 10, 2, 2]

Here is my take on it:

def poly_dev(poly, dev):
 if dev == 0:
 return poly[:]
 if dev > len(poly):
 return list()
 if dev < 0:
 raise ValueError("negative derivative")
 p = 1
 for k in range(2, dev+1):
 p *= k
 poly = poly[:-dev]
 n = len(poly)-1
 for i in range(len(poly)):
 poly[n-i] *= p
 p = p * (i+dev+1) // (i+1)
 return poly

So, the first thing I do after handling trivial cases is computing the multiplication factor for the first coefficient and take only a copy of the part of the list that I need (to avoid messing up the original).

Then I multiply each coefficient with the multiplication factor p, followed by computing the next one, meaning "divide by the smallest member of the product that defines p and multiply by the one bigger then the biggest one".

Notice that the division is //, which is integer division. That way you don't lose precision (and I think it's a bit quicker than the floating point one, but I'm not sure right now).

It may look redundant to multiply and later divide by the same number, but it's less work than multiplying dev numbers in each go.

Also, it would be more Pythonic to have enumerate(reversed(poly)) instead of range(len(poly)) in the for loop, but then we'd still need n-i somewhere, so this seemed cleaner to me for this particular case. Maybe I'm missing something obvious here.

• lnNoam: 1. I intentionally left out error handling, but thanks for including it in a final answer. 2. As this exercise is purely pedagogical, I'm rather curious as to why creating a copy of poly is not an expensive operation. I intentionally tried to avoid such a step, but that effort appears to have been in vain. Do you know why? Thanks.
• Vedran Šego: If you're OK with messing up the original list, you can replace the creation of the new list poly = poly[:-dev] with removing the list's elements: poly[-dev:] = list().
• Paul Wintz: Why use poly[n-i] *= p but not p *= (i+dev+1) in the second and third to last lines?
• Vedran Šego: @PaulrBear There is a subtle difference that is important here: p *= (i+dev+1) // (i+1) would've been equivalent to p = p * ((i+dev+1) // (i+1)), which would get you the wrong result in case i+dev+1 is not divisble by i+1. One could, ov course, do p *= i + dev + 1 and then p //= i + 1, but I prefer having a single expression.
• Paul Wintz: Ohhh, I see. Coming from a Java background, I interrupted // (i+1) as a comment. Thanks for the clarification!

Answer 2

I would use this function signature:

def poly_deriv(poly, nth=1):
 ...

"dev" is too cryptic for my taste. The second parameter should default to 1, as I would expect it to be the common case.

Your code works on poly in place, rather than on a copy, and thus trashes the contents of the array. That's unacceptable, especially considering that the result is shorter than the input, so a copy is unavoidable.

It is not obvious to me what c and r represent.

Counting loops are usually done using for ... in range(...). I think that it would be easier to work term by term, repeatedly differentiating it.

def poly_deriv(poly, nth=1):
 degree = len(poly) - 1
 if degree - nth < 0:
 return [0] # or an empty list? Your decision.
 poly = poly[: -nth]
 for i in range(len(poly)):
 for exponent in range(degree - i, degree - i - nth, -1):
 poly[i] = exponent * poly[i]
 return poly

As @GarethRees has pointed out in a comment, the little-endian convention for representing polynomials is clearly technically superior, since the index of each element corresponds to the degree. The big-endian representation just happens to match the human writing convention.

• jonny: Agreed wholeheartedly. I'm often confounded by the overuse of the contraction of words in programming - even poly_deriv is too contracted IMO. Writing polynomial_derivative is a difference of a few characters but it is immediately plain to anyone looking at the signature what it does.

Answer 3

Some notes:

• Following @GarethRees' sound advice, I'd definitely use reversed indexes in the polynomials.

• In programming, like in math, you take advantage of abstractions. Instead of a monolithic algorithm, you split it into meaningful functions. Here a product(start, end) function that multiplies a range of numbers seems useful.

• Use list-comprehensions instead of imperative loops.

For example:

import functools
import operator
def product(start, end):
 return functools.reduce(operator.mul, range(start, end + 1), 1) 
def poly_derivative(poly, n):
 if n < 0:
 raise ValueError("The nth order derivative must be non-negative")
 return [x * product(idx + 1, idx + n) for (idx, x) in enumerate(poly[n:])]

• oliverpool: It looks nice. Maybe product is a better name than mul (and instead of using range as an argument, you could provide def product(start, end):). You should also check that n is strictly greater than 0. I didn't know about the second argument of enumerate: maybe naming it is a good idea! enumerate(poly[n:], start=n + 1)
• tokland: All good points! updated