9

I have a need to check if a point in R^2 lies in a circle of relatively large radius r (up to 10^5). Obviously I would normally just compare the inner product to r^2, but this is in an embedded environment and this isn't going to work on int32_t values that are large enough since the quadratures will overflow the type (max 32 bit types).

Possible solutions:

I could manually kludge a 64 bit product out of two 32 bit ints (probably what I'll end up doing).

I could divide everything by 10 (or any value) then do the usual inner product comparison, but I loose precision.

I could try to check inside an n-gon inscribed in the circle, but that is a lot of calculation, tables, etc. and I still loose precision.

Is there an algorithm that is typically used for this sort of thing?

11
  • What's the range of r? Commented Jan 9, 2021 at 22:38
  • Do you care more about worst-case timing or average time for each point? Is code size a concern?
    – user3185968
    Commented Jan 9, 2021 at 22:39
  • 5
    If you're trying to decide whether x^2 + y^2 < r^2 is true, you can test if either x or y is greater than r fairly quickly, and whether both x and y are less than r/sqrt(2). Depending on how your x, y and r are distributed, this could be a useful fast test.
    – user3185968
    Commented Jan 9, 2021 at 22:51
  • 3
    You could do a crude initial test with a simple bounding box. If the point is outside the box, it's outside the circle. Commented Jan 9, 2021 at 22:52
  • 5
    When dealing with large values, don't compute x^2 + y^2 <= r^2, because the addition could overflow even when x<r and y<r. Instead compute x^2 <= r^2 - y^2. As long as y <= r the subtraction is guaranteed to work. Commented Jan 9, 2021 at 23:08

2 Answers 2

6

I'm afraid computing the 64-bit results is the simplest solution. Check if your compiler can generate efficient inline code for this:

int check_distance(int x, int y, int r) {
    return (long long)x * x + (long long)y * y <= (long long)r * r;
}

If the generated code seems too slow, you can add a test to check if 64-bit operation is required. Assuming x, y and r are positive, here is a solution using unsigned arithmetics and exact width types from <stdint.h>:

int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if (x <= 46340 && y <= 46340 && r <= 0xffff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x + y * y <= r * r;
    } else {
        return (uint64_t)x * x + (uint64_t)y * y <= (uint64_t)r * r;
    }
}

Notice the constant 46340 which is floor(sqrt(pow(2, 31))): If both x and y are greater than this value, x*x + y*y will exceed 232.

Here is an alternative with a quicker test, but that will fall back to 64-bit operation for slightly smaller values:

int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if ((x | y | r) <= 0x7fff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x + y * y <= r * r;
    } else {
        return (uint64_t)x * x + (uint64_t)y * y <= (uint64_t)r * r;
    }
}

Then if you really don't want to use the compiler's 64-bit arithmetics, you can write the computation explicitly. Considering the range of x, y and r specified as <= 100000, shifting the values right 2 bits keeps x and y below 46340:

int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if (x <= 46340 && y1 <= 46340 && r1 <= 0xffff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x + y * y <= r * r;
    } else {
        /* shift all values right 2 bits to keep them below 46340 */
        uint32_t x0 = x & 3;
        uint32_t y0 = y & 3;
        uint32_t r0 = r & 3;
        uint32_t x1 = x >> 2;
        uint32_t y1 = y >> 2;
        uint32_t r1 = r >> 2;
        uint32_t x2_lo = x0 * (x0 + x1 * 8);
        uint32_t y2_lo = y0 * (y0 + y1 * 8);
        uint32_t d2_lo = x2_lo + y2_lo;
        uint32_t d2_hi = x1 * x1 + y1 * y1 + (d2_lo >> 4);
        uint32_t r2_lo = r0 * (r0 + r1 * 8);
        uint32_t r2_hi = r1 * r1 + (r2_lo >> 4);
        return d2_hi < r2_hi || (d2_hi == r2_hi && (d2_lo & 15) <= (r2_lo & 15));
    }
}

Finally, shifting values by 5 bits allows for numbers up to 1000000:

int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if (x <= 46340 && y1 <= 46340 && r1 <= 0xffff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x + y * y <= r * r;
    } else {
        /* shift all values right 5 bits to keep them below 46340 */
        uint32_t x0 = x & 31;
        uint32_t y0 = y & 31;
        uint32_t r0 = r & 31;
        uint32_t x1 = x >> 5;
        uint32_t y1 = y >> 5;
        uint32_t r1 = r >> 5;
        uint32_t x2_lo = x0 * (x0 + x1 * 64);
        uint32_t y2_lo = y0 * (y0 + y1 * 64);
        uint32_t d2_lo = x2_lo + y2_lo;
        uint32_t d2_hi = x1 * x1 + y1 * y1 + (d2_lo >> 10);
        uint32_t r2_lo = r0 * (r0 + r1 * 64);
        uint32_t r2_hi = r1 * r1 + (r2_lo >> 10);
        return d2_hi < r2_hi || (d2_hi == r2_hi && (d2_lo & 1023) <= (r2_lo & 1023));
    }
}

All of the above versions produce exact results for the specified ranges. If you do not require exact result, you can just shift the values to bring them within the proper range and perform the 32-bit computation:

int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    while (x > 46340 || y > 46340 || r > 0xffff) {
        x >>= 1;
        y >>= 1;
        r >>= 1;
    }
    /* 32-bit unsigned expression no longer overflows */
    return x * x + y * y <= r * r;
}
7
  • Thanks! Yea, that's actually what I'm doing now. It's not bad - it just feels like there should be something easier. But then math is a cruel mistress sometimes. Commented Jan 9, 2021 at 22:46
  • @TrivialCase: Did you look at the assembly code generated for the above function? The processor should have a 32x32 -> 64 multiply opcode
    – chqrlie
    Commented Jan 9, 2021 at 22:51
  • Yes unfortunately it just overflows a 32, so I will just stick with cobbling together a 64 but value from two 32s split into 16 bit words. Commented Jan 9, 2021 at 23:19
  • chqrlie's answer is as good as it gets for checking "within a circle." A potential optimization is to check if -r <= x <= r and -r <= y <= r, before checkingx*x + y*y <= r*r. You would have to benchmark this to see if it makes a difference. With full retail optimizations turned on, it might not make any difference.
    – selbie
    Commented Jan 10, 2021 at 13:05
  • It's going to be very difficult to beat the compiler in generating assembly with the above code snippet.
    – selbie
    Commented Jan 10, 2021 at 13:06
2

The maximum radius you can have, using the premises you have described is 100000, which leads to a squared modulus of 10 000 000 000, which requires (as all numbers are positive) 35 bit integer (of unsigned square radius) to be represented.

Based on these premises, and the fact that you have no easy way to get 64bit integers, and having quite low extra bits, we can scale the results four bits in norm (two bits in the source coordinates) to achieve full capacity to handle upto 100000 coordinates in a 32bit unsigned integer.

In my first edition of this answer, I assumed using only one shift in the coordinates was enough to handle the full set of values (two bits in the calculated norm), and losing 1 bit of precision was considered, but I was wrong and one extra bit was needed. It is needed to shift the results at least three bits to hold the full set of inputs, so I decided to scale the coordinates two bits, and so the results would be scaled by four. As I decided also to always scale, and return the fraction of a square unit as an integer ranging from 0 to 15 (in sixteenths of a square unit). So you will achieve exact results by comparing first the integer parts of the two points and use the fractional parts, in case the integer parts match. This makes the computation and the meaning of the results returned back more coherent than earlier, and gives you complete exactitude with integer coordinates.

You requested a working implementation, so I have posted one for you below:

#include <stdio.h>
#include <stdint.h>

/* calculate the square of a divided by four number and
 * accumulate the fraction (in sixteenths of a square unit)
 * into the reference pointed by frac_p.  */
uint32_t
square_of_div16(uint32_t x, int *frac_p)
{
    /* we use (IP + FP)^2 = IP^2 + 2*IP*FP + FP^2 */

    uint32_t int_part    = x >>  2;                /* divide by four */
    uint32_t frac_part   = x & 0x3;                /* mod 4 */
    uint32_t int_result  = int_part  * int_part;   /* square of IP */
    int      frac_result = frac_part * frac_part;  /* square of FP */
    uint32_t mixed_prod  = int_part  * frac_part;  /* IP*FP */

    int_result  +=  mixed_prod >> 1;
    frac_result += (mixed_prod & 1) << 3;
    if (frac_result >= 0x10) { /* carry process */
        int_result  += frac_result >> 4;
        frac_result &= 0x0f;
    }
    if (frac_p) *frac_p += frac_result; /* accumulate */
    return int_result;
}

/* this calculates the squared norm scaled to one sixteenth
 * of the original coordinates (scaled by one fourth).
 * The ref_fraction pointer is a reference of a variable to
 * accumulate the fraction sixteenths of a square unit.  If
 * you are not interested in the fraction value, you can just
 * pass NULL as parameter. */
uint32_t
norm_scaled(uint32_t x, uint32_t y, int *ref_fraction)
{
    int fraction = 0;
    uint32_t result = 0;

    result += square_of_div16(x, &fraction);
    result += square_of_div16(y, &fraction);

    if (ref_fraction)
        *ref_fraction += fraction; /* the excess */

    return result;
}

/* TEST MAIN PROGRAM.  Just input pairs of coordinates in the
 * same line (separated by spaces) and calculate the squared
 * norm of the vector, scaled by 1/16 (accumulating the
 * fraction of the value in 1/16s of a square unit in the
 * location referenced.  This is done using double floating
 * point numbers and uint32_t integers. */
int main()
{
    char line[256];
    while (fgets(line, sizeof line, stdin) != NULL) {
        int x = 0, y = 0, fraction = 0;

        sscanf(line, "%u%u", &x, &y);

        uint32_t norm_16th = norm_scaled(x, y, &fraction);

        printf("Trying (%u, %u) => %u (fraction = %d/16)\n",
                x, y, norm_16th, fraction);

        double norm_sq_16th
            = (double) x/4.0 * (double)x/4.0
            + (double) y/4.0 * (double)y/4.0;

        printf("squared norm scaled: %.8f\n", norm_sq_16th);
    }
    printf("Program finished\n");
}

The function square_of_div16 calculates a scaled modulus divided by 16 of a number, so we can reuse it to calculate the squares of x and y coordinates. The function takes a pointer frac_p to an integer variable to store the fraction part (in sixteenths of a square unit)

The function norm_scaled then calculates the norm, by using the square_of_div16 function and adding both results. The fractional part is accumulated for both calls and the result is accumulated to the referred variable by pointer ref_fraction. A carry processing is done in here, to give correct results.

Finally a main() routine is in charge of querying the user to input pairs of coordinates and calculate the scaled norm of the resulting vector by calling the function and using the squares pithagorean formula applied to double values. The results should be the same in all cases.

6
  • 1
    This is really the answer I think. Any sane person would know that the reduction in precision is likely well within the margin for error of, say, the ADC that I am getting the coordinates from. In my case, I am on the hook to explain to the client the reduction in precision, so I might not actually implement it this way - but this is still a very good answer. Thanks! Commented Jan 10, 2021 at 18:07
  • Your approach is fine, but you should "normalize* the fraction in square_of_div16 so the resulting fraction is always <= 15: use int frac_result = frac_part * frac_part + *frac_p;. The test if (frac_result >= 0x10) is probably more costly than always performing the carry propagation. Also you did not post the function that checks if the point is within the circle.
    – chqrlie
    Commented Jan 11, 2021 at 23:24
  • @chqrlie as you see, i do normalization (i call it carry processing) after I add the numbers. Before is nonsense, as I'm accumulating the value of several sums. I do all at once in the last part, after all the sums have been done. The cost of making a bit shift against a full integer multiplication doesn't need to be demonstrated (what you do in your comment) The right shift is a division by 16 while the and mask is a remainder calculation. Do you really think your recommendation is going to be more efficient? don't believe so. Commented Jan 12, 2021 at 19:23
  • @chqrlie, you can see that the last operation I do with the fraction is to mask it with fraction &= 0x0f;, so it is going always to be a number in the range 0 to 15. The adding done with the external reference is not normalized, as any normalization should be done in the external procedure (it is possible that the calling process is doing several sums and want to normalize, after all sums have been done) This is so for two reasons: 1) to improve efficiency (only one carry proc is done at the end) and 2) overflowing an int with sums of numbers in the range 0 to 15 is not probable. Commented Jan 12, 2021 at 19:55
  • The last operation in square_of_div16() is not masking fraction with 0x0f, it is if (frac_p) *frac_p += frac_result; /* accumulate */ which may produce a value in *frac_p larger than 15. Adding the current of *frac_p to frac_result is harmless and cheap. and the final operation becomes *frac_p = frac_result;. Obviously you would need to remove the possibility to pass a null pointer for frac_p. With this approach the result is always normalized, which makes the comparison much simpler and avoids the pitfall of forgetting to perform a final normalisation.
    – chqrlie
    Commented Jan 12, 2021 at 20:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.