Plotting polynomials roots

Question

NOTE: See follow up to this question here

I created a simple python script to plot quadratic, cubic and quartic polynomials with integer coefficients between -4 and 4. It uses numpy to find the roots for the polynomials and matplotlib for the actual plotting of the points. Because I want to plot all possible polynomial with integer coefficients between -4 and 4 I figured I needed to do a recursive function with a loop. However if I'm not mistaken this will give a time complexity of \$O(n!)\$.

import matplotlib
# Use Qt4Agg because otherwise matplotlib crashes on my computer running Arch Linux
matplotlib.use('Qt4Agg')
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
min_degree = 2
max_degree = 3
colours = []
def polynomial_scatter(ax, min_x, max_x, min_degree, max_degree, cur_degree, coeff):
 global colours
 if cur_degree <= max_degree:
 for c in range(min_x, max_x):
 new_coeff = coeff + [c]
 polynomial_scatter(ax, min_x, max_x, min_degree, max_degree, cur_degree + 1, new_coeff)
 if cur_degree >= min_degree:
 roots = np.roots(new_coeff)
 # check if no possible solution exists
 if len(roots) == 0:
 continue
 # the real part of the roots will be our x values, the imaginary parts will be our y values
 x, y = zip(*[(t.real, t.imag) for t in roots])
 # scatter the roots for this polynomial with a colour corresponding to the degree of the polynomial
 ax.scatter(x, y, c=colours[cur_degree - max_degree], alpha=.5, s=1)
if __name__ == "__main__":
 f, ax = plt.subplots()
 colours = cm.rainbow(np.linspace(0, 1, max_degree - min_degree + 1))
 polynomial_scatter(ax, -4, 4, min_degree, max_degree, 0, [])
 plt.savefig("plot.svg", format="svg")
 plt.show()

Here are some pictures showing the result: enter image description here quartic polynomial cubic polynomial

if I set max_degree to 5 it takes between 7-15 hours on my computer which is very long. Perhaps there is a way to implement multithreading.

How can I improve the performance of the polynomial_scatter algorithm?

• mvds: I think there are some symmetries in the coefficient space you're exploring. If all coefficients are multiplied by a constant, the roots stay the same, right? This would mean that any solution (except for the trivial solution) is plotted at least twice; e.g. after negating coefficients. So you could choose to loop over range(0,max_x) for the last coefficient. Also I think that coeff[0]==0 implies a lower degree, correct? If so, the function should either call itself or calculate roots.
• mvds: Expanding on the previous comment: if the last coefficient is 0, you're actually solving a lower degree polynomial (ignoring the trivial solution). That implies that you could also loop over range(1,max_x). This should be 55% faster already.
• Linus: @mvds that is an interesting optimization. Not sure how to implement that into a recursive function however.
• mvds: in pseudo code: range(cur_degree==max_degree?1:min_x, max_x)
• Linus: @mvds I might implement that, but I switched to C++ and have another performance problem now. Here is the new post: codereview.stackexchange.com/questions/164790/…

Answer 1

I know that you subsequently posted c++ code for this and that post has answers, one of which is accepted.

I did notice that the code above references the global variable colours. This may be seen as bad programming style because it can cause side effects that may be challenging to detect. Alternatively, function parameters could be used, or else a class with a class/data/instance member.