Sciencemadness Discussion Board

Accurate solver of quartic equations.

woelen - 21-12-2023 at 10:45

I developed a solver for quartic equations (general 4th degree equations with 5 coefficients). I think it is one of the best solvers available worldwide. I did not find any quartic, which made it fail, and it is very fast (the Java implementation takes appr. 400 ns for solving a quartic equation at nearly the maximum accuracy, which is available theoretically).

A description can be found here: https://woelen.homescience.net/science/math/exps/accurate_so...

The website links to software and a paper. I intend to publish the paper in an open source journal. Still looking for a suitable place.

akmetal - 3-1-2024 at 00:57

Its not clear how Q(phi) is arrived at with the random free variable being introduced in the anti-diagonal. This feels very much like the first term of the Bessel function lol. The gamma function is really specific yet completely arbitrary ....
Quote: Originally posted by woelen  
I developed a solver for quartic equations (general 4th degree equations with 5 coefficients). I think it is one of the best solvers available worldwide. I did not find any quartic, which made it fail, and it is very fast (the Java implementation takes appr. 400 ns for solving a quartic equation at nearly the maximum accuracy, which is available theoretically).

A description can be found here: https://woelen.homescience.net/science/math/exps/accurate_so...

The website links to software and a paper. I intend to publish the paper in an open source journal. Still looking for a suitable place.


[Edited on 3-1-2024 by akmetal]

annaandherdad - 26-12-2024 at 19:09

woelen, did you ever publish your paper on finding the roots of quartics?

I had a research project once that required me to find the roots of quartics. I followed the advice of Numerical Recipes, which was to map the problem into one of linear algebra, in which the quartic equation was the characteristic polynomial of a 4x4 matrix, so that the roots were the eigenvalues. The latter were determined by standard matrix operations. This worked well in my application.

annaandherdad - 28-12-2024 at 16:56

Hi woelen, I was looking at your home science web page on roots of polynomials. It is very nice that you have given such a careful and pedagogical explanation of things! It's also very cool that you have identified a hole in the literature (on how to avoid loss of precision in the use of the closed solutions for cubics and quartics), and filled it.

I have some comments about the part on quadratic equations. First, perhaps you should note that you are thinking of the case in which a,b,c are all real. It's clear from the context that that's what you mean, but it might help to say it explicitly.

Next, near the end of the section on quadratics, you say that there is a "little issue", namely, that the code gives a zero divide if b=c=0. I don't see that, looking over your code. It looks to me like your code should give the right answer (x1=x2=0) in that case. Am I missing something?

woelen - 29-12-2024 at 06:03

Thanks for your comments. I'll make it more explicit that I am considering the case of real roots.

The issue with b and c both equal to 0 is in the line of code

xre[1] = c/(a*xre[0]);

The variable c equals zero in this case, and xre[0] also equals zero. So, you get 0/0 in that code, which leads to a so-called NaN. So, the special case that both b and c are equal to zero must be handled separately, by simply setting both solutions to 0 and then returning from the solve-method..

annaandherdad - 2-1-2025 at 08:08

Quote: Originally posted by woelen  


The issue with b and c both equal to 0 is in the line of code

xre[1] = c/(a*xre[0]);

The variable c equals zero in this case, and xre[0] also equals zero. So, you get 0/0 in that code, which leads to a so-called NaN. So, the special case that both b and c are equal to zero must be handled separately, by simply setting both solutions to 0 and then returning from the solve-method..


Thanks for the clarification, yes, I see the issue. It's amazing how complicated even the usual quadratic formula is.

I spent quite some time reading and rereading the papers by you and by Kahan. It took me a while to realize that Wilco Oelen was actually you. Have you succeeded in getting your paper published? It should be publishable, you have corrected some errors in the existing literature and provided new analysis and new software. In the meantime, why don't you put it on the arXiv? It's not refereed, but it does allow you to establish a kind of precedence for the ideas.

I know Kahan personally, although not that well. I talked to him for several hours once about 10 years ago, concerning a project of mine. You know he played a major role in the establishment of the IEEE standards for floating point, including the use of denormalized numbers that smooth out the jump at overflow or underflow. He also played a major role in the design of HP calculators. At the time he complained to me about the designers of GPU's, who were evidently not following the IEEE standards because they took too much hardware to implement in detail. The logic was, if it's only for gaming, who cares if the answers come out wrong once in a while? The problem, of course, is that people started to use GPU's for designing bridges and other applications. Now they play a big role in AI, which is why Nvidida's stock has exploded. Kahan also complained about Bill Gates, who left no opportunity for the users of Windows to set the floating point rounding option.

In the case of the quadratic formula, there are 3 sources of potential error. One, the loss of precision in computing -b +- sqrt(b^2 - 4ac) (one of the two signs has a loss of precisiion); two, the loss of precisiion in computing b^2 - 4ac when roots are close together; three, overflow/underflow. Your notes address the first of these, which is the easiest to fix. Kagan has a long discussion of the second of these, and I spent a long time thinking about it. He also talks about the third.

I've used the cubic solution of Cardano and Ferrari, but not much. It was fun learning about the quartic solution.

woelen - 4-1-2025 at 06:46

Putting it on arXiv may indeed be an option. Is there some procedure for doing that?

Interesting to have met Kahan personally. I read quite a few things from him and I know he has contributed to the IEEE standards for floating point numbers.

In my software library, I address all three kinds of potential issues you mention.
The computation of b*b - 4ac I do with the help of so-called fused multiply add instructions, which can be used to compute a*b + c without intermediate rounding. It is a single instruction, accessible from Java by means of the Math.fma intrinsic code. It performs exceptionally well. Using two fma-calls and some bookkeeping, you can also compute a*b + c*d without intermediate rounding and this is what I do in computing the value of b*b - 4ac.
Underflow/overflow also is covered in my software library: I do that by scaling a, b, and c in a particular way, using Java's Math.scalb method, which allows exact scaling of numbers, without adding numerical noise, which may make arithmetic less accurate. I perform scaling, then use the abc-formula (together with the trick, mentioned in my web page), and the result is scaled back. Using this strategy, I can solve equations with coefficients which are very close to the representable range of an IEEE double.

In all my code, also for higher degree equations, and in the generic polynomial equation solvers, I do my best to avoid overflow/underflow, but I try to keep things reasonable. If e.g. my code works fine within a range of 1e-250 to 1e250, instead of 1e-308 to 1e308 then it is good enough for me. The quadratic solver is very close to the range of IEEE doubles, but the other solvers work in a somewhat smaller range, but levels of 1e250 still are handled fine. In practice that kind of input hardly occurs.

annaandherdad - 4-1-2025 at 19:25

Quote: Originally posted by woelen  
Putting it on arXiv may indeed be an option. Is there some procedure for doing that?


You should be able just to go to the arXiv web site and set up an account. It is free, and posting manuscripts is free. There is some basic checking for scientific integrity but you should have no trouble with that. https://arxiv.org/

I did not know about the fused multiply-add instruction. It shows how out of date I am. Thanks for explaining it. I read the wikipedia article on it and found a link to another interesting article by Kahan, "How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?" Like his other article (the one you posted), this one is unpublished.

[Edited on 5-1-2025 by annaandherdad]

woelen - 6-1-2025 at 04:27

The fma-instructions are really interesting. I use them for computing a*b+c*d without any intermediate rounding, but they can be used for computing general dot-products and Horner-evaluations of polynomials as well.

Performance also is quite good, well-crafted code is only about 1.5 times as slow as standard precision code: https://www-pequan.lip6.fr/~graillat/papers/posterRNC7.pdf