C algorithm for calculating norm/inner product

Question

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?

question from:https://stackoverflow.com/questions/65648300/c-algorithm-for-calculating-norm-inner-product

Answer 1

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 to 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 2³².

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;
}