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@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@jyu.fi
University of Jyväskylä
Department of Mathematical Information Technology
-------------------------------------------------