Newton's Cubic Equation Solver code

[dsilvavini]

Hello, everybody.

I'm trying to change some code from a library 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 that I couldn't understand yet. I want to understand it to verify errors later on. Any help would be very appreciated.

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;

	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;

	// 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.