计算范数/内积的C算法

Question

我需要检查R ^ 2中的点是否位于半径r较大(最大为10 ^ 5)的圆上。显然，我通常只将内部乘积与r ^ 2进行比较，但这是在嵌入式环境中进行的，并且不会对int32_t值起作用，因为整数会溢出类型(最大32位类型)，因此int32_t值足够大。
可能的解决方案:
我可以手动从两个32位整数中提取一个64位乘积(可能最终会做)。
我可以将所有内容除以10(或任何值)，然后进行通常的内部乘积比较，但我失去了精度。
我可以尝试在圆内刻一个n形内进行检查，但这需要大量的计算，表格等操作，但我仍然无法满足精度要求。
是否有通常用于此类事情的算法？

Answer 1

恐怕计算64位结果是最简单的解决方案。检查您的编译器是否可以为此生成高效的内联代码:

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

如果生成的代码似乎变慢，则可以添加测试以检查是否需要64位操作。假设 x， y和 r是正数，这是一种使用无符号算术和 <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;
    }
}

注意常数46340，它是 floor(sqrt(pow(2, 31))):如果 x和 y都大于此值，则 x*x + y*y将超过232。
这是一种测试速度更快的替代方法，但是对于较小的值，它将退回到64位操作:

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

然后，如果 真的是不想使用编译器的64位算术，则可以显式编写计算。考虑到指定为 x的 y， r和 <= 100000的范围，将值右移2位可将 x和 y保持在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));
    }
}

最后，将值移位5位可允许最大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));
    }
}

上述所有版本都会在指定范围内产生准确的结果。如果您不需要精确的结果，则可以移动这些值以使它们处于适当的范围内并执行32位计算:

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