Newton's Cubic equation solver

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);}




			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);


			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.


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; = ad, = bd; = 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.

Apparently min{1,2,3} and max{1,2,3} are designed to take one of two values, 0 or 1. The “^1” inverts the least significant bit of an integer, thus turning 0 into 1, and 1 into 0.

The expression ((unsigned char)&coeffs+7) extracts the most significant byte of coeffs, the subsequent right shift extracts the sign bit, which is the most significant bit in that byte.

IMHO, it would have been clearer to use signbit(), but the code may have been written prior to its introduction in C99. This would result in the following code:

min3=!!signbit(coeffs[3]); max3=min3^1;        

  min2=!!signbit(coeffs[2]); max2=min2^1;

  min1=!!signbit(coeffs[1]); max1=min1^1;

Since signbit() returns a non-zero result (not necessarily 1) when the sign bit of the argument is set, the “!!” is needed to ensure the output of the expression is either 0 or 1.

I understood the signal of the coefficients changes the behavior of the cubic function, so the left and right input values have different meanings depending on the case. But now I have problems to visualize this part:

// 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];

First I thought those lines compute the value of the function for the left and right of the interval. Now I think it isn’t true, since min3, min2 and min1 could have different values. In this case the function computed will even not be exactly the one desired. I’m a little confused about this… Could someone help me?

Thanks in advance.

Ok, now I understood it. To visualize the minor and major values I had to think of the cubic function as a sum of simpler four functions. The first is


which is a trivial cubic function.

The second is

coeffs[2] * v2

which is a trivial parabolic function.

The third is


which is a trivial line function.

And the last is


, a constant function.

Since the values of v* are always positive (they express time), the minor value in the interval of each simpler function will be the left one ( l ) if the respective coefficient is positive and the right one ( r ) if negative. Similarly the major value will be r if the coefficient is positive and l if it is negative. So the sum of these values will be the minor/major of the more complex cubic too.