时间戳

首页 > 解决方案 > 计算范数/内积的C算法

c - 计算范数/内积的C算法

问题描述

我需要检查 R^2 中的一个点是否位于半径 r 相对较大（最多 10^5）的圆中。显然，我通常只会将内积与 r^2 进行比较，但这是在嵌入式环境中，这不适用于足够大的 int32_t 值，因为求积会溢出类型（最多 32 位类型）。

可能的解决方案：

我可以从两个 32 位整数中手动组合出一个 64 位产品（可能是我最终会做的）。

我可以将所有内容除以 10（或任何值），然后进行通常的内积比较，但我失去了精度。

我可以尝试检查圆圈中内接的 n 边形，但这是很多计算、表格等，我仍然失去精度。

是否有通常用于此类事情的算法？

标签： calgorithmmath

解决方案

恐怕计算 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则将超过 2 32。

这是一个更快测试的替代方法，但是对于稍小的值，这将回退到 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，将值右移 2 位保持并低于 46340：yr<= 100000xy

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

推荐阅读