Hello, everybody.
I was trying to convert a CPU code from a Continuous Collision Detection library to CUDA code when stepped on a implementation of a cubic equation solver based on Newton’s numerical method. In the beginning of the solver function there are some strange bitwise operations (for me at least) that I couldn’t understand yet. They have to do with the function interval considered by the Newton method but I really dont understand the hack. It’s probably a newbie question, but I need to get it to test the CUDA code later on.
inline bool solveCubicWithIntervalNewton(double &l, double &r, vec3f& baryc, bool bVF,
NewtonCheckData &data, double coeffs[])
{
double v2[2]={l*l,r*r};
double v[2]={l,r};
double rBkUp;
//Strange bitwise operations
unsigned char min3, min2, min1, max3, max2, max1;
min3=*((unsigned char*)&coeffs[3]+7)>>7;max3=min3^1;
min2=*((unsigned char*)&coeffs[2]+7)>>7;max2=min2^1;
min1=*((unsigned char*)&coeffs[1]+7)>>7;max1=min1^1;
//Bitwise operations end
// bound the cubic
double minor= coeffs[3]*v2[min3]*v[min3] + coeffs[2]*v2[min2] + coeffs[1]*v[min1] + coeffs[0];
double major= coeffs[3]*v2[max3]*v[max3] + coeffs[2]*v2[max2] + coeffs[1]*v[max1] + coeffs[0];
if (major<0) return false;
if (minor>0) return false;
// starting here, the bounds have opposite values
double m=0.5*(r+l);
// bound the derivative
double dminor= 3.0*coeffs[3]*v2[min3]+ 2.0*coeffs[2]*v[min2]+coeffs[1];
double dmajor= 3.0*coeffs[3]*v2[max3]+ 2.0*coeffs[2]*v[max2]+coeffs[1];
if ((dminor>0)||(dmajor<0)) // we can use Newton
{
double m2=m*m;
double fm= coeffs[3]*m2*m + coeffs[2]*m2 + coeffs[1]*m + coeffs[0];
double nl=m;
double nu=m;
if (fm>0) {nl-=fm*(1.0/dminor);nu-=fm*(1.0/dmajor);}
else {nu-=fm*(1.0/dminor);nl-=fm*(1.0/dmajor);}
//intersect with [l,r]
if (nl>r) return false; // pas de solution
if (nu<l) return false; // pas de solution
if (nl>l)
{
if (nu<r) {l=nl;r=nu;m=0.5*(l+r);}
else {l=nl;m=0.5*(l+r);}
}
else
{
if (nu<r) {r=nu;m=0.5*(l+r);}
}
}
// sufficient temporal resolution, check root validity
if ((r-l)<ccdTimeResolution)
if (bVF)
return checkRootValidity_VF(r, baryc, data);
else
return checkRootValidity_EE(r, baryc, data);
rBkUp = r, r = m;
if (solveCubicWithIntervalNewton(l,r,baryc, bVF, data, coeffs)) return true;
l = m, r = rBkUp;
return (solveCubicWithIntervalNewton(l,r,baryc, bVF, data, coeffs));
}
And this is the function that call it. To give some backgound: The code tries to get the first time of contact between a moving triangle and a point in a time interval. ta0, tb0, tc0 and q0 are repectively the triangle vertices and point positions on the interval beginning. Similarly ta1, tb1, tc1 and q1 means the same for the interval ending. The other arguments are to store information of the solution and aren’t important as the _equateCubic_VF function wich calculates the coeficients of the cubic equation.
float
Intersect_VF(const vec3f &ta0, const vec3f &tb0, const vec3f &tc0,
const vec3f &ta1, const vec3f &tb1, const vec3f &tc1,
const vec3f &q0, const vec3f &q1,
vec3f &qi, vec3f &baryc)
{
/* Default value returned if no collision occurs */
float collisionTime = -1.0f;
vec3f qd, ad, bd, cd;
/* diff. vectors for linear interpolation */
qd = q1 - q0, ad = ta1 - ta0, bd = tb1 - tb0, cd = tc1 - tc0;
/*
* Compute scalar coefficients by evaluating dot and cross-products.
*/
float a, b, c, d; /* cubic polynomial coefficients */
_equateCubic_VF(ta0, ad, tb0, bd, tc0, cd, q0, qd, a, b, c, d);
if (IsZero(a) && IsZero(b) && IsZero(c) && IsZero(d))
return collisionTime;
NewtonCheckData data;
data.a0 = ta0, data.b0 = tb0;
data.c0 = tc0, data.p0 = q0;
data.ad = ad, data.bd = bd;
data.cd = cd, data.pd = qd;
/*
* iteratively solve the cubic (scalar) equation and test for validity of the solution.
*/
double l = 0;
double r = 1;
double coeffs[4];
coeffs[3] = a, coeffs[2] = b, coeffs[1] = c, coeffs[0] = d;
if (solveCubicWithIntervalNewton(l, r, baryc, true, data, coeffs)) {
collisionTime = (l+r)*0.5f;
qi = qd*collisionTime + q0; //pont in the time of collision
}
return collisionTime;
}
Thanks in advance,
VinÃcius da Silva.