# Torbert, 9.2.2005

from pylab import *
from math import sqrt

def qsolve(a,b,c):
	rad=sqrt(b**2-4*a*c)
	if 4*abs(a*c) < b**2:
		if b > 0:
			return -(b+rad)/(2*a), -2*c/(b+rad)
		else:
			return (-b+rad)/(2*a), 2*c/(-b+rad)
			
	else:
		return (-b+rad)/(2*a),(-b-rad)/(2*a)

def f(x):
	r1,r2 = qsolve(1.0,9**12,-3.0)
	return (x-r1)*(x-r2)

r1,r2 = qsolve(1.0,9**12,-3.0)
print (r1, f(r1))
print (r2, f(r2))
x=arange(-4.0,4.0,0.1)
y=map(f,x)
plot(x,y)
plot([r1,r2],map(f,[r1,r2]),"o")
show()