woelen
Super Administrator
Posts: 8043
Registered: 20-8-2005
Location: Netherlands
Member Is Offline
Mood: interested
|
|
Accurate solver of quartic equations.
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
Harmless
Posts: 43
Registered: 20-11-2023
Member Is Offline
|
|
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
Hazard to Others
Posts: 391
Registered: 17-9-2011
Member Is Offline
Mood: No Mood
|
|
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.
Any other SF Bay chemists?
|
|
annaandherdad
Hazard to Others
Posts: 391
Registered: 17-9-2011
Member Is Offline
Mood: No Mood
|
|
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?
Any other SF Bay chemists?
|
|
woelen
Super Administrator
Posts: 8043
Registered: 20-8-2005
Location: Netherlands
Member Is Offline
Mood: interested
|
|
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
Hazard to Others
Posts: 391
Registered: 17-9-2011
Member Is Offline
Mood: No Mood
|
|
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.
Any other SF Bay chemists?
|
|
woelen
Super Administrator
Posts: 8043
Registered: 20-8-2005
Location: Netherlands
Member Is Offline
Mood: interested
|
|
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
Hazard to Others
Posts: 391
Registered: 17-9-2011
Member Is Offline
Mood: No Mood
|
|
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]
Any other SF Bay chemists?
|
|
woelen
Super Administrator
Posts: 8043
Registered: 20-8-2005
Location: Netherlands
Member Is Offline
Mood: interested
|
|
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
|
|