CFD Online Discussion Forums - Bug in eigenValue(const tensor& t)

I have found a subtle bug in the eigenValue(const tensor& t) function. It has the following code:

...
else
{
scalar Q = (a*a - 3*b)/9;
scalar R = (2*a*a*a - 9*a*b + 27*c)/54;

scalar R2 = sqr(R);
scalar Q3 = pow3(Q);

// Three different real roots
if (R2 < Q3)
{
scalar sqrtQ = sqrt(Q);
scalar theta = acos(R/(Q*sqrtQ));
...

It is assumed here that if R2 < Q3 then |R|<=Q*sqrtQ, so that the argument of acos never exceeds 1.0 by absolute value. But it is not always the case. Consider the following values:

R = 00111011001011011110100111010100100101111010111101 01101010000111
= (8419873611799175)/(680564733841876926926749214863536422912) ~ 1.23719e-23,

Q = 00111100110000110100010111001001101110000010100100 01101110001100 = (1356189309486819)/(2535301200456458802993406410752) ~= 5.34922e-16.

Here the first value is 64 bits representing the corresponding double value according to the IEEE 754 binary64 format and the second is the exact value of these 64 bits in the form of a rational number. The last value is the approximation of this rational number in scientific notation. We have

R2 = 00110110011010111111011001111101111110000101001100 01101001001101 = (7870845268728397)/(5142201741628768881734278695491720328071049580104 9370729644032),

Q3 = 00110110011010111111011001111101111110000101001100 01101001001110 = (3935422634364199)/(2571100870814384440867139347745860164035524790052 4685364822016),

so that

R2-Q3 = -(1)/(5142201741628768881734278695491720328071049580104 9370729644032) < 0,

and R2 < Q3. On the other hand, we have

Q*sqrtQ = 00111011001011011110100111010100100101111010111101 01101010000110 = (4209936805899587)/(340282366920938463463374607431768211456),

and

R - Q*sqrtQ = (1)/(680564733841876926926749214863536422912) > 0,

and thus R/(Q*sqrtQ) > 1.0. In fact, we have

R/(Q*sqrtQ) = 00111111111100000000000000000000000000000000000000 00000000000001 = (4503599627370497)/(4503599627370496) > 1.

Most probably, the simplest possible solution of this problem is to replace the call

acos(R/(Q*sqrtQ))

by the call

acos(sqrt(R2/Q3))

It is not clear that there are no new problems with this approach, but at least it avoids arguments outside the range [-1, 1].