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 exact

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.