[Numpy-discussion] Cython-based OpenMP-accelerated quartic polynomial solver

[Juha Jeronen]

Hi,

On 30.09.2015 04:37, Charles R Harris wrote:
>
>
> On Tue, Sep 29, 2015 at 6:48 PM, Chris Barker - NOAA Federal 
> <chris.barker at noaa.gov <mailto:chris.barker at noaa.gov>> wrote:
>
>     This sounds pretty cool -- and I've had a use case. So it would be
>     nice to get into Numpy.
>
>     But: I doubt we'd want OpenMP dependence in Numpy itself.
>
>     But maybe a pure Cython non-MP version?
>
>     Are we putting Cuthon in Numpy now? I lost track.
>
>
> Yes, but we need to be careful of Numeric Recipes.

True.

As I said, only the algorithms come from NR, and the code was written 
from scratch. I haven't even looked at any of their code, precisely due 
to the restrictive license.


In this particular case, the referred pages (pp. 227-229, 3rd ed.) 
contain only a description of the solution process in English, with 
equations - there is no example code printed into the book in section 
5.6, Quadratic and Cubic Equations, which in its entirety is contained 
on pp. 227-229.

Furthermore, on p. 228, equation (5.6.12), which contains the solutions 
of the cubic equation, is attributed to Francois Viete. The reference 
given is the treatise "De emendatione", published in 1615.

The same solution, also with attribution to Francois Viete, is given in 
Wikipedia, but without especially defining temporary variables to 
eliminate some common subexpressions:

https://en.wikipedia.org/wiki/Cubic_function#Trigonometric_.28and_hyperbolic.29_method


The solution to the quadratic equation is the standard one, but with 
special care taken to avoid catastrophic cancellation in computing the 
other root.

Wikipedia mentions that in the standard formula, this issue exists:

https://en.wikipedia.org/wiki/Quadratic_equation#Avoiding_loss_of_significance

but does not give an explicit remedy. References given are

Kahan, Willian (November 20, 2004), On the Cost of Floating-Point 
Computation Without Extra-Precise Arithmetic. ( 
http://www.cs.berkeley.edu/~wkahan/Qdrtcs.pdf , see esp. p. 8, which 
provides the same remedy as NR )

Higham, Nicholas (2002), Accuracy and Stability of Numerical Algorithms 
(2nd ed.), SIAM, p. 10, ISBN 978-0-89871-521-7.


(Browsing through Kahan's text, I see now that I could improve the 
discriminant computation to handle marginal cases better. Might do this 
in a future version.)


For both of these cases, I've used NR as a reference, because the 
authors have taken special care to choose only numerically stable 
approaches - which is absolutely critical, yet something that sadly not 
all mathematics texts care about.


The quartic equation is not treated in NR. The reduction trick is 
reported in Wikipedia:

https://en.wikipedia.org/wiki/Quartic_function#Solving_by_factoring_into_quadratics

where the original reference given is

Brookfield, G. (2007). "Factoring quartic polynomials: A lost art". 
Mathematics Magazine 80 (1): 67–70.

(This seemed an elegant approach, given stable solvers for cubics and 
quadratics.)


  -J

-------------------------------------------------
Juha Jeronen, Ph.D.
juha.jeronen at jyu.fi
University of Jyväskylä
Department of Mathematical Information Technology
-------------------------------------------------

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20150930/202f29df/attachment.html>