Plot roots

import numpy as np
import itertools

def squareRoots( x ):
    return np.roots( [1,0,-x] )
def cubeRoots( x ):
    return np.roots( [1,0,0,-x] )

xs = squareRoots(50)
ys = [ y for x in xs for y in cubeRoots(7+x) ]
solutions = [ a[0]+a[1] for a in itertools.product( ys, ys ) ]
uniqueSolutions = np.unique( np.array( solutions ) )
print "Found",len( solutions ),"out of which",len( uniqueSolutions ),"are unique:"
for sol in uniqueSolutions:
    print "    "+str( sol )

import matplotlib.pyplot as plt
fig = plt.figure( figsize = (12,12) )
ax = fig.add_subplot( 111, xlim=[-6,6], ylim=[-6,6] )
ax.set_aspect( 'equal', 'datalim' )

ax.scatter( [ x.real for x in uniqueSolutions ], [ x.imag for x in uniqueSolutions ], s = 100, color = 'r', zorder = 5, label = r"$(7 \pm \sqrt{50})^{\frac{1}{3}} + (7 \pm \sqrt{50})^{\frac{1}{3}}$" )
ax.scatter( [ x.real for x in np.unique( np.array(ys) ) ], [ x.imag for x in np.unique( np.array(ys) ) ], s = 50, color = '0.5', zorder = 4, label = r"$(7 \pm \sqrt{50})^{\frac{1}{3}}$" )

# round to make np.unique work with rounding errors
rs  = np.unique( np.array( [ abs(x) for x in uniqueSolutions ] ).round( decimals = 12 ) )
rys = np.unique( np.array( [ abs(x) for x in ys              ] ).round( decimals = 12 ) )
print "rs  =",rs
print "rys =",rys
for r in rs:
    ax.add_artist( plt.Circle( (0,0), r, edgecolor = 'r', alpha = 0.5, linestyle = '--', fill = False, linewidth = 2 ) )
for y in ys:
  for ry in rys:
    ax.add_artist( plt.Circle( ( y.real, y.imag ), ry, edgecolor = '0.5', fill = False, linewidth = 2 ) )
ax.legend()
fig.tight_layout()
fig.savefig( "cuberoot-solutions.png" )

plt.show()