Subscribe to GameDev.net Direct to receive the latest updates and exclusive content.
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.
Posted 07 June 2007 - 07:52 PM
Posted 07 June 2007 - 09:45 PM
Posted 07 June 2007 - 10:06 PM
Quote:
Original post by Zipster
Considering that a quartic is the highest degree polynomial with a closed-form solution . . .
Posted 07 June 2007 - 10:41 PM
Quote:
Original post by dragongame Quote:
Original post by Zipster
Considering that a quartic is the highest degree polynomial with a closed-form solution . . .
That is not true! You can go up to 4th degree, but as stated the solution will be very complicated. I'm sorry but I can't give you the solutions without my math book.
Posted 08 June 2007 - 02:40 AM
#include <math.h>
#include <complex>
const double PI = 3.14159265358979323846;
//----------------------------------------------------------------------------
bool solveQuadratic(double &a, double &b, double &c, double &root)
{
if(a == 0.0 || abs(a/b) < 1.0e-6)
{
if(abs(b) < 1.0e-4)
return false;
else
{
root = -c/b;
return true;
}
}
double discriminant = b * b - 4.0 * a * c;
if(discriminant >= 0.0)
{
discriminant = sqrt(discriminant);
root = (b + discriminant) * -0.5 / a;
return true;
}
return false;
}
//----------------------------------------------------------------------------
bool solveCubic(double &a, double &b, double &c, double &d, double &root)
{
if(a == 0.0 || abs(a/b) < 1.0e-6)
return solveQuadratic(b, c, d, root);
double B = b/a, C = c/a, D = d/a;
double Q = (B*B - C*3.0)/9.0, QQQ = Q*Q*Q;
double R = (2.0*B*B*B - 9.0*B*C + 27.0*D)/54.0, RR = R*R;
// 3 real roots
if(RR<QQQ)
{
/* This sqrt and division is safe, since RR >= 0, so QQQ > RR, */
/* so QQQ > 0. The acos is also safe, since RR/QQQ < 1, and */
/* thus R/sqrt(QQQ) < 1. */
double theta = acos(R/sqrt(QQQ));
/* This sqrt is safe, since QQQ >= 0, and thus Q >= 0 */
double r1, r2, r3;
r1 = r2 = r3 = -2.0*sqrt(Q);
r1 *= cos(theta/3.0);
r2 *= cos((theta+2*PI)/3.0);
r3 *= cos((theta-2*PI)/3.0);
r1 -= B/3.0;
r2 -= B/3.0;
r3 -= B/3.0;
root = 1000000.0;
if(r1 >= 0.0) root = r1;
if(r2 >= 0.0 && r2 < root) root = r2;
if(r3 >= 0.0 && r3 < root) root = r3;
return true;
}
// 1 real root
else
{
double A2 = -pow(fabs®+sqrt(RR-QQQ),1.0/3.0);
if (A2!=0.0) {
if (R<0.0) A2 = -A2;
root = A2 + Q/A2;
}
root -= B/3.0;
return true;
}
}
//----------------------------------------------------------------------------
bool solveQuartic(double &a, double &b, double &c, double &d, double &e, double &root)
{
// When a or (a and b) are magnitudes of order smaller than C,D,E
// just ignore them entirely.
if(a == 0.0 || abs(a/b) < 1.0e-6 || abs(a/c) < 1.0e-6)
return solveCubic(b, c, d, e, root);
// Uses Ferrari's Method
double aa = a*a, aaa=aa*a, bb=b*b, bbb=bb*b;
double alpha = -3.0*bb/(8.0*aa) + c/a, alpha2 = alpha * alpha;
double beta = bbb/(8.0*aaa) + b*c/(-2.0*aa) + d/a;
double gamma = -3.0*bbb*b/(256.0*aaa*a) + c*bb/(16.0*aaa) + b*d/(-4.0*aa) + e/a;
if(beta == 0.0)
{
root = b/(-4.0*a) + sqrt(0.5 * (-alpha + sqrt(alpha2 + 4.0*gamma)));
return true;
}
else
{
std::complex<double> P = -alpha2/12.0 - gamma;
std::complex<double> Q = -alpha2*alpha/108.0 + alpha*gamma/3.0 - beta*beta/8.0;
std::complex<double> R = Q*0.5 + sqrt(Q*Q*0.25 + P*P*P/27.0);
std::complex<double> U = pow(R, 1.0/3.0);
std::complex<double> y = -5.0*alpha/6.0 - U;
if(U != 0.0) y += P/(3.0*U);
std::complex<double> W = sqrt(alpha + y + y);
std::complex<double> aRoot;
bool foundRealRoot = false;
double firstPart = b/(-4.0*a);
std::complex<double> secondPart = -3.0*alpha - 2.0*y;
std::complex<double> thirdPart = 2.0*beta/W;
aRoot = firstPart + 0.5 * (-W - sqrt(secondPart + thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0)
{
root = aRoot.real();
foundRealRoot = true;
}
aRoot = firstPart + 0.5 * (-W + sqrt(secondPart + thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0 &&
(!foundRealRoot || aRoot.real() < root))
{
root = aRoot.real();
foundRealRoot = true;
}
aRoot = firstPart + 0.5 * (W - sqrt(secondPart - thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0 &&
(!foundRealRoot || aRoot.real() < root))
{
root = aRoot.real();
foundRealRoot = true;
}
aRoot = firstPart + 0.5 * (W + sqrt(secondPart - thirdPart));
if(abs(aRoot.imag()) < 1.0e-10 && aRoot.real() >= 0.0 &&
(!foundRealRoot || aRoot.real() < root))
{
root = aRoot.real();
foundRealRoot = true;
}
return foundRealRoot;
}
}
//----------------------------------------------------------------------------
Posted 08 June 2007 - 03:15 AM
Quote:
Original post by dragongame Quote:
Original post by Zipster
Considering that a quartic is the highest degree polynomial with a closed-form solution . . .
That is not true! You can go up to 4th degree, but as stated the solution will be very complicated. I'm sorry but I can't give you the solutions without my math book.
Posted 08 June 2007 - 03:20 AM
Posted 08 June 2007 - 03:40 AM
Quote:
Original post by pTymN
[EDIT] I realized that I should probably note that in any case where there are multiple roots, it tries to return the root that is closest to zero but above it. So if I had, say, root={-1, 1, 0.2}, then the algorithm returns 0.2. Assuming that ABCDE are terms that describe explicitly the movement of the two spheres over a given frame, then 0.2 is the most useful root, since it says that the spheres were able to use 20% of their motions.
Quote:
Original post by pTymN
Bonus note:
I'm not as sure how the math works out, but for my circle-circle test (2D collision test), I was able to reduce it to second order just by fixing one of the circles at the origin, and making its radius zero. The other circle had all of the motion and radius, so it ended up being a whenDoesCircleTouchOrigin test, which was much faster than 4th order.
Posted 08 June 2007 - 03:54 AM
Posted 08 June 2007 - 04:05 AM
Quote:
I'm a little confused about that, because with my previous collision detection using position and velocity (second degree), although there were only two, I always took the - root (like from the +- sqrt()). I assumed when doing to the fourth degree, it would be only one of the roots that I'd use aswell.
Posted 08 June 2007 - 10:03 AM
Posted 08 June 2007 - 12:20 PM
Posted 09 June 2007 - 03:13 AM
Quote:
Original post by Silver Phoenix
Solve equation for z (one solution is enough):
z^3 + 2Iz^2 + (I^2 - 4K)z + J^2 = 0
Posted 09 June 2007 - 12:38 PM
Quote:
Original post by johnnyBravo
I just kinda learned how to do the cubic equation, but I don't understand how there are only 3 roots, because I have the two +- from the quadratic equation,
but from that I'd assume you would be able to use the ++, --, +- and -+
Quote:
Original post by johnnyBravo
btw x as in ax^3 + bx^2 + cx + d
Quote:
Original
z^3 + 2Iz^2 + (I^2 - 4K)z + J^2 = 0
Posted 09 June 2007 - 01:42 PM
Posted 10 June 2007 - 02:33 AM
Quote:
Original post by pTymN
If you do a more thorough investigation and find a useful pattern, I'd appreciate if you'd share your learnings. :-)
Posted 13 June 2007 - 12:14 PM
bool solveQuadraticOther(double a, double b, double c, double &root)
{
if(a == 0.0 || abs(a/b) < 1.0e-4)
{
if(abs(b) < 1.0e-4)
return false;
else
{
root = -c/b;
return true;
}
}
double discriminant = b * b - 4.0 * a * c;
if(discriminant >= 0.0)
{
discriminant = sqrt(discriminant);
root = (b - discriminant) * -0.5 / a;
return true;
}
return false;
}
//----------------------------------------------------------------------------
bool solveQuartic(double a, double b, double c, double d, double e, double &root)
{
// I switched to this method, and it seems to be more numerically stable.
// http://www.gamedev.net/community/forums/topic.asp?topic_id=451048
// When a or (a and b) are magnitudes of order smaller than C,D,E
// just ignore them entirely. This seems to happen because of numerical
// inaccuracies of the line-circle algorithm. I wanted a robust solver,
// so I put the fix here instead of there.
if(a == 0.0 || abs(a/b) < 1.0e-5 || abs(a/c) < 1.0e-5 || abs(a/d) < 1.0e-5)
return solveCubic(b, c, d, e, root);
double B = b/a, C = c/a, D = d/a, E = e/a;
double BB = B*B;
double I = -3.0*BB*0.125 + C;
double J = BB*B*0.125 - B*C*0.5 + D;
double K = -3*BB*BB/256.0 + C*BB/16.0 - B*D*0.25 + E;
double z;
bool foundRoot2 = false, foundRoot3 = false, foundRoot4 = false, foundRoot5 = false;
if(solveCubic(1.0, I+I, I*I - 4*K, -(J*J), z))
{
double value = z*z*z + z*z*(I+I) + z*(I*I - 4*K) - J*J;
double p = sqrt(z);
double r = -p;
double q = (I + z - J/p)*0.5;
double s = (I + z + J/p)*0.5;
bool foundRoot = false, foundARoot;
double aRoot;
foundRoot = solveQuadratic(1.0, p, q, root);
root -= B/4.0;
foundARoot = solveQuadratic(1.0, r, s, aRoot);
aRoot -= B/4.0;
if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0)
|| root < 0.0)) || (!foundRoot && foundARoot))
{
root = aRoot;
foundRoot = true;
}
foundARoot = solveQuadraticOther(1.0, p, q, aRoot);
aRoot -= B/4.0;
if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0)
|| root < 0.0)) || (!foundRoot && foundARoot))
{
root = aRoot;
foundRoot = true;
}
foundARoot = solveQuadraticOther(1.0, r, s, aRoot);
aRoot -= B/4.0;
if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0)
|| root < 0.0)) || (!foundRoot && foundARoot))
{
root = aRoot;
foundRoot = true;
}
return foundRoot;
}
return false;
}
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.
GameDev.net™, the GameDev.net logo, and GDNet™ are trademarks of GameDev.net, LLC.