#!/usr/bin/env python3 import math import cmath import numpy def quad_discriminating_roots(a,b,c, entier = 1e-5): """For reference, the naive algorithm which shows complete loss of precision on the quadratic in question. (This function also returns a characterization of the roots.)""" discriminant = b*b - 4*a*c a,b,c,d =complex(a), complex(b), complex(c), complex(discriminant) root1 = (-b + cmath.sqrt(d))/2./a root2 = (-b - cmath.sqrt(d))/2./a if abs(discriminant) < entier: return "real and equal", abs(root1), abs(root1) if discriminant > 0: return "real", root1.real, root2.real return "complex", root1, root2 def middlebrook(a, b, c): try: q = math.sqrt(a*c)/b f = .5+ math.sqrt(1-4*q*q)/2 except ValueError: q = cmath.sqrt(a*c)/b f = .5+ cmath.sqrt(1-4*q*q)/2 return (-b/a)*f, -c/(b*f) def whatevery(a, b, c): try: d = math.sqrt(b*b-4*a*c) except ValueError: d = cmath.sqrt(b*b-4*a*c) if b > 0: return div(2*c, (-b-d)), div((-b-d), 2*a) else: return div((-b+d), 2*a), div(2*c, (-b+d)) def div(n, d): """Divide, with a useful interpretation of division by zero.""" try: return n/d except ZeroDivisionError: if n: return n*float('inf') return float('nan') testcases = [ (3, 4, 4/3), # real, equal (3, 2, -1), # real, unequal (3, 2, 1), # complex (1, -1e9, 1), # ill-conditioned "quadratic in question" required by task. (1, -1e100, 1), (1, -1e200, 1), (1, -1e300, 1), ] print('Naive:') for c in testcases: print("{} {:.5} {:.5}".format(*quad_discriminating_roots(*c))) print('\nMiddlebrook:') for c in testcases: print(("{:.5} "*2).format(*middlebrook(*c))) print('\nWhat Every...') for c in testcases: print(("{:.5} "*2).format(*whatevery(*c))) print('\nNumpy:') for c in testcases: print(("{:.5} "*2).format(*numpy.roots(c)))