14
\$\begingroup\$

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?

\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – mvds Jun 2 '17 at 0:11
  • \$\begingroup\$ 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. \$\endgroup\$ – mvds Jun 2 '17 at 0:14
  • \$\begingroup\$ @mvds that is an interesting optimization. Not sure how to implement that into a recursive function however. \$\endgroup\$ – Linus Jun 2 '17 at 8:06
  • \$\begingroup\$ in pseudo code: range(cur_degree==max_degree?1:min_x, max_x) \$\endgroup\$ – mvds Jun 2 '17 at 15:06
  • \$\begingroup\$ @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/… \$\endgroup\$ – Linus Jun 2 '17 at 15:21
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.