2

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.

user3556814
  • 21
  • 1
  • 4
  • `-Ofast` probably inlined the `sqrt` call. Compile (specifically, link) with `-lm` and you should be good – Nemo Nov 01 '17 at 21:14
  • Does [this answer](https://stackoverflow.com/a/25439802/5522303) from the question you linked not help? – Kevin Nov 01 '17 at 21:14
  • Insure code first uses `#include ` – chux - Reinstate Monica Nov 01 '17 at 21:14
  • "algorithm picks two random points (the x and y coordinates) between zero and one." and `x = 2*randf() - 1;` is not consistent. `x` can have negative values too between `[-1 0]`. – chux - Reinstate Monica Nov 01 '17 at 21:19
  • @Nemo, I already link libm. That is not the issue. – user3556814 Nov 01 '17 at 21:22
  • @chux , I can assure you that randf() only returns a number between zero and one. – user3556814 Nov 01 '17 at 21:24
  • 3
    You might think you're linking libm, but it's telling you quite clearly you're not. Listen to your compiler. One thing that's an issue sometimes is that the `-lm` has to come last on the `gcc` command. Also, maybe the compiler and/or libraries aren't installed right? All the other stuff here is distraction. Find the library! – Lee Daniel Crocker Nov 01 '17 at 21:26
  • what if you use spherical coordinates? – MFisherKDX Nov 01 '17 at 21:27
  • @user3556814 "randf() only returns a number between zero and one" is fuzzy on the edges. Does it create 0.0? 1.0? Yet an important concern is that with the possibilities of FP optimizations, `FLT_EVAL_METHOD`, does the generation of `mag, point[0], point[1], point[2]` given any set of results from `randf()` generate an unacceptable value set? Near the edges, FP math differs importantly from true math. – chux - Reinstate Monica Nov 01 '17 at 21:37
  • Something about this code stinks. If x^2 + y^2 <= 1 then the fabs is basically useless unless to handle edge cases around x^2 +y^2 ~ 1 only. Also, mag is guaranteed to be between 0 and 1 and equal to the value z such that x^2 + y^2 + z^2 ~ 1 ... So (x, y, mag) is on the +z side of the unit sphere. But (2.x.mag, 2.y.mag and 1-2.mag^2) has basically no relation to the unit sphere. To see this, let x ~ y ~ sqrt(1/3). Then mag ~ sqrt(1/3) and the point (2/3, 2/3, 1/3) is not on the unit sphere. Or x ~ y ~ sqrt(2) gives point (0, 0, 0). What gives? – Patrick87 Nov 01 '17 at 23:37
  • Ignore sqrt(1/3) example and use sqrt (1/2) only. Sqrt(1/3) actually turns out OK – Patrick87 Nov 01 '17 at 23:56
  • @Patrick87: There are edge cases where `x*x+y*y` is exactly 1 due to rounding but `1-x*x-y*y` is negative. The `fabs` deals with that. However, a better solution may have been to use the same expression in both places, either `1-x*x-y*y < 0` in the loop test or `1-(x*x+y*y)` in the `sqrt`. – Eric Postpischil Nov 02 '17 at 01:01
  • I don't want to argue with you about whether or not you linked libm... but in any case your code does not produce uniformly distributed points. It is a miraculous property of spheres that you can pick a uniformly distributed point by choosing a uniform random z, and then a uniform random angle around the z axis to get x,y. You don't need a rejection loop. – Matt Timmermans Nov 02 '17 at 01:14
  • "I already link libm. That is not the issue" No you don't and yes it is – Nemo Nov 02 '17 at 02:38

3 Answers3

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.

Jens
  • 69,818
  • 15
  • 125
  • 179
  • "The x and y must satisfy x * x + y * y <= 1 in order for there to be a solution for z." --> is not that insured with `while (x*x + y*y > 1);`? – chux - Reinstate Monica Nov 01 '17 at 21:17
  • @chux Well, maybe. The code disagrees with the text, I'd say. Or the text is extremely ambiguous. – Jens Nov 01 '17 at 21:18
  • I checked this box of very early in my trouble shooting process. I already link libm. That is not the issue. And you are correct about the math stuff. – user3556814 Nov 01 '17 at 21:21
  • Then you checked off that box prematurely. Everything else is distraction--get the compiler command pointed at the correct library. – Lee Daniel Crocker Nov 01 '17 at 21:28
  • Is it just me or is this the only answer that actually attempts to answer the question? – SirGuy Nov 01 '17 at 21:45
  • @LeeDanielCrocker, I know it's linked to the correct library because I wrote a test function using just sqrt() and it worked fine. – user3556814 Nov 01 '17 at 22:05
  • It's still something about the build process, not the code. You're no longer using `-Ofast` or `--fast-math`, right? – Lee Daniel Crocker Nov 01 '17 at 22:24
  • I still want to see your Makefile/other build script. – Lee Daniel Crocker Nov 01 '17 at 22:41
  • To a mathematician, the magnitude of x and y, considered as a pair, is the magnitude of the vector [x, y], hence sqrt(x\*x+y\*y). If they wanted to talk about their magnitudes separately they would say “their magnitudes” rather than “their magnitude.” Although the sentence is wrong because it says “less than” but should say “less than or equal to.” – Eric Postpischil Nov 01 '17 at 22:42
  • @user3556814: A test function using `sqrt` might compile and link even without the correct libm linked for various reasons, such as it took the square root of a constant, so the compiler evaluated it at compile time and did not generate a run-time reference to `sqrt`, or because its result was not used, so the compiler discarded it. – Eric Postpischil Nov 01 '17 at 22:45
0

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;
}
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
0

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);
MFisherKDX
  • 2,840
  • 3
  • 14
  • 25
  • Note: Using spherical coordinates does generate a different distribution of points on a sphere than OPs original approach: perhaps better, perhaps worse. Yet the final _computed_ `|x,y,z|` still may exceed 1.0 by just a tad. If that is an issue for OP, a loop to re-generate on excessive magnitude would easily solve that. – chux - Reinstate Monica Nov 01 '17 at 21:46