sqrt(1.0 - pow(1.0,2)) returns -nan [duplicate]

This question already has answers here: Is floating point math broken? (31 answers) Why pow(10,5) = 9,999 in C++ (8 answers) 开发者_如何学Python Closed 4 years ago.

I've found an interesting floating point problem. I have to calculate several square roots in my code, and the expression is like this:

sqrt(1.0 - pow(pos,2))

where pos goes from -1.0 to 1.0 in a loop. The -1.0 is fine for pow, but when pos=1.0, I get an -nan. Doing some tests, using gcc 4.4.5 and icc 12.0, the output of

1.0 - pow(pos,2) = -1.33226763e-15

and

1.0 - pow(1.0,2) = 0

or

poss = 1.0
1.0 - pow(poss,2) = 0

Where clearly the first one is going to give problems, being negative. Anyone knows why pow is returning a number smaller than 0? The full offending code is below:

int main() {
  double n_max = 10;
  double a = -1.0;
  double b = 1.0;
  int divisions = int(5 * n_max);
  assert (!(b == a));

  double interval = b - a;
  double delta_theta = interval / divisions;
  double delta_thetaover2 = delta_theta / 2.0;
  double pos = a;
  //for (int i = 0; i < divisions - 1; i++) {
   for (int i = 0; i < divisions+1; i++) {

    cout<<sqrt(1.0 - pow(pos, 2)) <<setw(20)<<pos<<endl;

     if(isnan(sqrt(1.0 - pow(pos, 2)))){
      cout<<"Danger Will Robinson!"<<endl;
      cout<< sqrt(1.0 - pow(pos,2))<<endl;
      cout<<"pos "<<setprecision(9)<<pos<<endl;
      cout<<"pow(pos,2) "<<setprecision(9)<<pow(pos, 2)<<endl;
      cout<<"delta_theta "<<delta_theta<<endl;
      cout<<"1 - pow "<< 1.0 - pow(pos,2)<<endl;
      double poss = 1.0;
      cout<<"1- poss "<<1.0 - pow(poss,2)<<endl;


  }

  pos += delta_theta;

}

 return 0;
 }

---

When you keep incrementing pos in a loop, rounding errors accumulate and in your case the final value > 1.0. Instead of that, calculate pos by multiplication on each round to only get minimal amount of rounding error.

---

The problem is that floating point calculations are not exact, and that 1 - 1^2 may be giving small negative results, yielding an invalid sqrt computation.

Consider capping your result:

double x = 1. - pow(pos, 2.);
result = sqrt(x < 0 ? 0 : x);

or

result = sqrt(abs(x) < 1e-12 ? 0 : x);

---

setprecision(9) is going to cause rounding. Use a debugger to see what the value really is. Short of that, at least set the precision beyond the possible size of the type you're using.

---

You will almost always have rounding errors when calculating with doubles, because the double type has only 15 significant decimal digits (52 bits) and a lot of decimal numbers are not convertible to binary floating point numbers without rounding. The IEEE standard contains a lot of effort to keep those errors low, but by principle it cannot always succeed. For a thorough introduction see this document

In your case, you should calculate pos on each loop and round to 14 or less digits. That should give you a clean 0 for the sqrt.

You can calc pos inside the loop as

pos = round(a + interval * i / divisions, 14);

with round defined as

double round(double r, int digits) 
{
    double multiplier = pow(digits,10);
    return floor(r*multiplier + 0.5)/multiplier;
}