73 lines
2 KiB
Python
73 lines
2 KiB
Python
#!/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)))
|