4 ms·

Nicely visualized. One thing however that I realized is that conceptually you actually need a second step if a point falls exactly on the line with distance 1

Nicely visualized. One thing however that I realized is that conceptually you actually need a second step if a point falls exactly on the line with distance 1 to the origin. Half of these points should count to the circle, not all of them, or none of them. With 32 bit floating point values there won't be many of these of course...

afafsd 12y ago
Would you care to estimate the magnitude of this error, and how many steps you would need before it actually became relevant?
- darkmighty 12y ago
  I believe the error should be (roughly) ~ size of f.p. interval at 1 * pi ~ pi* 2^(-32/64), i.e. barely representable. Number of steps would be ~ 1/(error^2) ~ 2^(2* 32/64). Even giving 50% probability at the line the result would be minutely biased because of many other things involving the discreteness of the floating point arithmetic. A less biased version would have to maybe increase the resolution around the boundary as number of steps progresses.
  - MrQuincle 12y ago
    Thanks! Practically, it might be fun to have a small library that artificially reduces the precision of floating point numbers. I've never heard about such a thing though.
darkmighty 12y ago
The line is a 0-measure space: it shouldn't matter in theory but if your boundary is "squiggly" enough it indeed does, in practice.
stephencanon 12y ago
Every primitive Pythagorean triple has c odd, so strictly speaking (+/-1,0) and (0,+/-1) are the only points with binary FP coordinates that have distance exactly 1 from the origin (decimal FP is more interesting in this regard). Of course, there are more points with binary FP coordinates whose magnitude computed naively as sqrt(xx + yy) produces a rounded result of 1.0 in binary floating-point; (3/5,4/5) is one such example. I would address this by using careful computation of the residual of xx + yy - 1.0 to determine which side of the circle the point lies on; something like the following sketch: float residual = x*x + y*y - 1 if (residual == 0) { /* x*x + y*y rounded to 1; we need to do a more careful check to determine whether we're inside or outside. WLOG, assume that |x| >= |y|. */ if (fabsf(x) < fabsf(y)) { float tmp = x; x = y; y = tmp; } /* We know that 0.5 <= x*x < 1.0, from which it follows that fmaf(x,x,-1) is computed without any rounding; we can then do another fma to isolate the residual (this is also exact because we know that y*y cannot be *too* tiny unless x is exactly 1.0, in which case the first fma is exactly zero). */ residual = fmaf(y,y,fmaf(x,x,-1)); } if (residual < 0) { /* inside circle */ } else if (residual == 0) { /* on circle */ } else { /* outside circle */ } It's worth noting that fma isn't actually necessary here; it suffices to evaluate xx + yy in double precision, except in the |x| = 1 case, which is easy to detect. fma is conceptually cleaner (and faster on modern architectures that support it), however.

afafsd 12y ago

Would you care to estimate the magnitude of this error, and how many steps you would need before it actually became relevant?

darkmighty 12y ago
I believe the error should be (roughly) ~ size of f.p. interval at 1 * pi ~ pi* 2^(-32/64), i.e. barely representable. Number of steps would be ~ 1/(error^2) ~ 2^(2* 32/64). Even giving 50% probability at the line the result would be minutely biased because of many other things involving the discreteness of the floating point arithmetic. A less biased version would have to maybe increase the resolution around the boundary as number of steps progresses.
- MrQuincle 12y ago
  Thanks! Practically, it might be fun to have a small library that artificially reduces the precision of floating point numbers. I've never heard about such a thing though.

darkmighty 12y ago

The line is a 0-measure space: it shouldn't matter in theory but if your boundary is "squiggly" enough it indeed does, in practice.

stephencanon 12y ago

Every primitive Pythagorean triple has c odd, so strictly speaking (+/-1,0) and (0,+/-1) are the only points with binary FP coordinates that have distance exactly 1 from the origin (decimal FP is more interesting in this regard). Of course, there are more points with binary FP coordinates whose magnitude computed naively as sqrt(xx + yy) produces a rounded result of 1.0 in binary floating-point; (3/5,4/5) is one such example. I would address this by using careful computation of the residual of xx + yy - 1.0 to determine which side of the circle the point lies on; something like the following sketch: float residual = x*x + y*y - 1 if (residual == 0) { /* x*x + y*y rounded to 1; we need to do a more careful check to determine whether we're inside or outside. WLOG, assume that |x| >= |y|. */ if (fabsf(x) < fabsf(y)) { float tmp = x; x = y; y = tmp; } /* We know that 0.5 <= x*x < 1.0, from which it follows that fmaf(x,x,-1) is computed without any rounding; we can then do another fma to isolate the residual (this is also exact because we know that y*y cannot be *too* tiny unless x is exactly 1.0, in which case the first fma is exactly zero). */ residual = fmaf(y,y,fmaf(x,x,-1)); } if (residual < 0) { /* inside circle */ } else if (residual == 0) { /* on circle */ } else { /* outside circle */ } It's worth noting that fma isn't actually necessary here; it suffices to evaluate xx + yy in double precision, except in the |x| = 1 case, which is easy to detect. fma is conceptually cleaner (and faster on modern architectures that support it), however.