Archived

This topic is now archived and is closed to further replies.

avianRR

polynomial solution needed (3rd & 4th order)

Recommended Posts

avianRR    100
I need to write functions to solve 3rd and 4th order polynomials. The function proto types are as follows... struct POLYRESULT3{ double Result[3]; bool Real[3]; }; struct POLYRESULT4{ double Result[4]; bool Real[4]; }; POLYRESULT3 SolveForT(double A,double B,double C,double D); // solve At^3 + Bt^2 + Ct + D = 0 POLYRESULT4 SolveForT(double A,double B,double C,double D,double E); // solve At^4 + Bt^3 + Ct^2 + Dt + E = 0 My problem with writing these two functions is that I don''t know how to solve these types of equations and haven''t been able to find any examples/tutorials anywhere. Any help would be appreciated.

Share this post


Link to post
Share on other sites
v71    100
you may use Ruffini''s rule ( as we call in italy ) but this may yelds to complex solutions if you look in a algebra book you''ll find how to solve this equations

Share this post


Link to post
Share on other sites
avianRR    100
quote:
Original post by CorsairK8
Well, perhaps you can write a trial''n''error based search algorithm. You can end it once it yields an answer to the correct precision.


Tried it, The problem is that it took about 10 min to get the answers for one problem. My program has to be able to solve the problem millions of times in a reasonable amount of time.

quote:
Original post by narfidoo
http://mathforum.com/dr.math/faq/faq.cubic.equations.html


Thanks, that was very interesting. I still don''t understand how to solve the problems but that page has enough info for me to study that it''ll take me a while to go through it all.

Share this post


Link to post
Share on other sites
Enigma    1410
Hi,

A little while ago I wrote a fairly fast polynomial solver. It only works with polynomials containing only positive integer powers of x, but it looks like that''s what you want. It just uses Newton-Raphson method. This means you can modify the accuracy at the expense of some speed. I used it while messing around with the collision detection code in the NeHe tutorial No. 31 (by Dimitrios Christopoulos) so I''m fairly confident that it will never fail to solve the equation. You should be able to optimise it further if you need to.

Only problem is that I never got round to wrapping it up in a simple method, it''s fairly integrated into my collision detection code at the moment. I''ll try and strip it out and post it later tonight if I can.

Share this post


Link to post
Share on other sites
Enigma    1410
Ok, i''ve wrapped my polynomial solver up into a class. It seems to work (minimal testing ''cause it took longer than I expected). I won''t post it yet because a) I haven''t finished commenting it properly and b) it''s fairly big! If you want it, post back and I''ll e-mail it to you. Speedwise, a quick test shows it should be able to solve about 750,000 quartic equations per minute on a 933MHz Pentium III. I probably still have some error checking to do, but that shouldn''t bring the speed down drastically.

Share this post


Link to post
Share on other sites
avianRR    100
quote:
Original post by Enigma
Ok, i''ve wrapped my polynomial solver up into a class. It seems to work (minimal testing ''cause it took longer than I expected). I won''t post it yet because a) I haven''t finished commenting it properly and b) it''s fairly big! If you want it, post back and I''ll e-mail it to you. Speedwise, a quick test shows it should be able to solve about 750,000 quartic equations per minute on a 933MHz Pentium III. I probably still have some error checking to do, but that shouldn''t bring the speed down drastically.


That would be fast enough for a begining. I could research a faster way while I use that. my e-mail is avian@prairie.lakes.com

Share this post


Link to post
Share on other sites
Guest Anonymous Poster   
Guest Anonymous Poster
Look at the first gem in "Graphics Gems V", "Solving Quartics and Cubics for graphics" It covers this in a lot of detail, with calculation steps that would be easy to turn into code and duscussions of the many issues of doing this.

Share this post


Link to post
Share on other sites
Guest Anonymous Poster   
Guest Anonymous Poster
If I recall my college calculas, there is no formula that will always provide the right answer to those equations. Now I havn''t had a chance to look into this in depth, and that may have just been the professors and books way of putting off till later what they were too lazy to explain right now, but I remember a few examples of attempted formulas. Most of them worked most of the time, but there would always be some trivial example (ie E = sqrt[-1] ) that would cause it to fail.


Okay, so that example isn''t trivial, but I was trying to make a point. I don''t have my calc book next to me right now anyways.

Share this post


Link to post
Share on other sites
Guest Anonymous Poster   
Guest Anonymous Poster
The last poster was incorrect, cubics always have a real root at least. Qurtics sometimes do not but then what''s wrong with comples results (e.g. sqrt(-1) = i)? Only eqautions of power 5 or above do not have a general solution (this has been proven) and would require approximation methods.

Share this post


Link to post
Share on other sites
Enigma    1410
Just a quick thought - I was looking at your structs and I see that you have an array to show which results are real and which are not. Do you need the imaginary roots? My class currently just returns an array of all real roots to the equation.

Share this post


Link to post
Share on other sites
Enigma    1410
Ok, finally got it finished. I''ll e-mail it to you in a couple of minutes. A tiny bit of optimisation pushed the quartics/second up to just over 20k. If you do use it I''d appreciate a bit of feedback (i.e any bugs, failure to find roots, obvious improvements you can see that could be made). Either post on the forum or e-mail me (address in profile).

Cheers & hope the code helps

Enigma

Share this post


Link to post
Share on other sites
FallingFrog    122
Oddly enough, I have been wrestling with the exact same problem the last couple days! I need to find whether
ax^3 + bx^2 + cx + d exceeds a certain cutoff value v anywhere within the interval x + q and x - q, (where q is an even power of 2). x is known in this case. I thought of finding the max value on the interval, and testing that, and I thought of trying to find whether there exist any roots of the function
ax^3 + bx^2 + cx + (d - v)
This would also give me the answer. But I have to know whether the roots are between x + q and x - q, or the answer is useless. This does not seem to be a trivial problem.

Actually I lied, I actually want to find whether the sum of
ax^3 + bx^2 + cx + d +
ey^3 + fy^2 + gy +
hz^3 + iz^2 + jz
is greater than the cutoff anywhere. But it''s a similar problem. I can afford to be approximate, though- for instance if there was an easy way to detect whether it was impossible, or certain that it was greater than the cutoff. It has to be fast- a few adds, a couple shifts, but no multiplies. It has to run about 5 million times per second.

FallingFrog

Share this post


Link to post
Share on other sites
Enigma    1410
FallingFrog - when you say you can afford to be approximate, how approximate do you mean? I''ve quickly put together a short program that uses floats for your coefficients and tries to solve five million problems and it runs in seven seconds on my machine (problem is there are a lot of approximations, some of them may make the answer invalid). I haven''t had time to test it to check that the output is valid, but I think it should be.

Share this post


Link to post
Share on other sites
FallingFrog    122
The program is automatically going to subdivide q into 8 subvolumes if it thinks that there exist any solutions. So the checker needs to overestimate: A false positive is ok, but a false negative is bad.

FallingFrog

Share this post


Link to post
Share on other sites
FallingFrog    122
By the way, currently I''m using fixed point integers as the coefficients, because then I can shift by log q instead of multiplying by q to find the value of the function at a point.

FallingFrog

Share this post


Link to post
Share on other sites
FallingFrog    122
I think I may have come up with a way of doing it. If I take the first corner, and find where the tangent plane at that corner intersects the nearest 2 midpoints of the sides of the square, and the center of the square, that subsquare might have a hit. So, recurse. Then do the same for the other 3 subsquares. Only do it in 3 dimensions. Is that roughly how your program works?

FallingFrog

Share this post


Link to post
Share on other sites
Enigma    1410
No, that''s not how my program works at all. (Actually, I can''t really figure out what you''re doing).

My program does what you were thinking about doing to start with - find the maxima and minima of the curve and test these and the start and end points against v. I''ve now got it down to 2.1 seconds for 5 million repeats.

Just out of curiosity, why do you need to be able to do this so many times? You seemed to state earlier that you were recursing and testing the interval again. Would it not be faster just to find the values where the curve crosses the line v, or would this not work for what you''re trying to do?

Share this post


Link to post
Share on other sites
FallingFrog    122
I''m going to render the volume in 3d. I have an input ray to a bounded volume (an octree node) and I know the formula for the cubic function inside the volume of the node, and I want to find where the line intersects the surface that you get wherever the function equals the cutoff point. So, I figure, for each subtree of the node, if there are any points inside where the value of the function exceeds the cutoff, trace the ray inside that node. Otherwise, skip it, and do the same for the next subnode that the ray hits.

FallingFrog

Share this post


Link to post
Share on other sites
Guest Anonymous Poster   
Guest Anonymous Poster

y^3 + p*y^2 + q*y + r = 0
Set y = (x - p/3) and we have

x^3 + a*x + b = 0 were
a = (3*q - P^2)/3
b = (2*p^3 - 9*p*q + 27*r)/27

Let
A = cubicroot(sqrt(b^2/4 + ^3/27) - b/2)
B = - cubicroot(sqrt(b^2/4 + a^3/27) + b/2)

And the roots are

A+B
sqrt(-3)*(A - B)/2 - (A + B)/2
sqrt(-3)*(B - A)/2 - (A + B)/2

If p,q,r are real then

If (b^2/4 + a^3/27).GT.0 there one real root and a pair
of conjugate complex ones

If (b^2/4 + a^3/27).EQ.0 there three real roots and
at least two equal

If (b^2/4 + a^3/27).LT.0 there three real roots and all
unequal

All from CRC Math Table... (any errors is for sure mine)

But you got some to do to make use of it ...

micca

Share this post


Link to post
Share on other sites
ragonastick    134
Newton''s Method is pretty cool (and it works for pretty much anything). It uses an iterative approach, so you can be as accurate or inaccurate as you want.

NewX = OldX - (f(OldX) / f''(oldX))

Provided that f''(OldX) isn''t 0.

What it actually does, is starting from a point, the tangent to the curve at that point is found, then at the point where that tangent hits the x axis, is the new x coordinate. This will eventually reach (or converge on) a solution most of the time. Problems come when you can''t get the solution you are looking for, in which case, your best bet is to find the turning points and other useful points which are often found near roots (points of inflexion) and use them for initial values of x.

The other possible problem can come when f''(oldX) is a really small number for a few values of x which totally puts off the search because the tangent is almost horizontal and intercepts the x-axis a long way away.

And, if you don''t know calculus, then this will be a bit harder. But here is the crash course:

Polynomial f(x) = ax^n + bx^m + cx^o + dx...

Derivative f''(x) = anx^(n - 1) + bmx^(m - 1) + cox^(o - 1) + d

(so you take the power, multiply it by the coefficient and then subtract 1 from the power)

Trying is the first step towards failure.

Share this post


Link to post
Share on other sites
Enigma    1410
I''m a bit busy at the moment, but I''ll try and post my source some time tonight.

FallingFrog - my cutoff code will now run 50M times in 10.6 seconds (1.06s for 5M times). I know this is a bit slower than what you wanted, but I can''t see any further way to speed it up (at least not without knowing what kind of values it will be given and trying to approximate based on these).

Share this post


Link to post
Share on other sites