Issue with Square Root in C Algorithm

问题:

I have a bit of code that finds a point on a unit sphere. Recall, for a unit sphere:

1 = sqrt( x^2 + y^2 + z^2 )

The algorithm picks two random points (the x and y coordinates) between zero and one. Provided their magnitude is less than one we have room to define a third coordinate by solving the above equation for z.

void pointOnSphere(double *point){
double x, y;

do {
    x = 2*randf() - 1;
    y = 2*randf() - 1;
} while (x*x + y*y > 1);

double mag = sqrt(fabs(1 - x*x - y*y));

point[0] = 2*(x*mag);
point[1] = 2*(y*mag);
point[2] = 1 - 2*(mag*mag);
}

Technically, I inherited this code. The previous owner compiled using -Ofast which "Disregards strict standards compliance". TL;DR it means your code doesn't need to follow strict IEEE standards. So when I tried to compile without optimization I ran into an error.

 undefined reference to `sqrt'

What are IEEE standards? Well, because computers can't store floating point numbers to infinite precision, rounding errors pop up during certain calculations if you're not careful.

After some googling I ran into this question which got me on the right track about using proper IEEE stuff. I even read this article about floating point numbers (which I recommend). Unfortunately it didn't answer my questions.

I'd like to use sqrt() in my function as opposed to something like Newton Iteration. I understand the issue in my algorithm probably comes from the fact I could potentially (even though not really) pass a negative number to the sqrt() function. I'm just not quite sure how to remedy the issue. Thanks for all the help!

Oh, and if it's relevant I'm using a Mersenne Twister number generator.

Just to clarify, I am linking libm with -lm! I have also confirmed it is pointing to the correct library.

回答1:

As for the undefined reference to sqrt you need to link with libm, usually with -lm or similar option.

Also note that

Provided their magnitude is less than one we have room to define a third coordinate by solving the above equation for z.

is wrong. The x and y must satisfy x * x + y * y <= 1 in order for there to be a solution for z.

回答2:

I'd use spherical coordinates

theta = randf()*M_PI;
phi = randf()*2*M_PI;
r = 1.0;
x = r*sin(theta)*cos(phi);
y = r*sin(theta)*sin(phi);
z = r*cos(theta);

回答3:

To insure the points meet a condition, test for the condition itself as part of the while loop, rather than a derivation of the condition.

// functions like `sqrt(), hypot()` benefit with declaration before use
//   and without it may generate "undefined reference to `sqrt'"
// Some functions like `sqrt()` are understood and optimized out by a smart compiler.
// Still, best to always declare them.
#include <math.h>

void pointOnSphere(double *point){
  double x, y, z;
  do {
    x = 2*randf() - 1;
    y = 2*randf() - 1;
    double zz = 1.0 - hypot(x,y); 
    if (zz < 0.) continue;  // On rare negative values due to imprecision
    z = sqrt(zz);
    if (rand()%2) z = -z;  // Flip z half the time
  } while (x*x + y*y + z*z > 1);  // Must meet this condition

  point[0] = x;
  point[1] = y;
  point[2] = z;
}