我需要检查R ^ 2中的点是否位于半径r较大(最大为10 ^ 5)的圆上。显然,我通常只将内部乘积与r ^ 2进行比较,但这是在嵌入式环境中进行的,并且不会对int32_t值起作用,因为整数会溢出类型(最大32位类型),因此int32_t值足够大。
可能的解决方案:
我可以手动从两个32位整数中提取一个64位乘积(可能最终会做)。
我可以将所有内容除以10(或任何值),然后进行通常的内部乘积比较,但我失去了精度。
我可以尝试在圆内刻一个n形内进行检查,但这需要大量的计算,表格等操作,但我仍然无法满足精度要求。
是否有通常用于此类事情的算法?
最佳答案
恐怕计算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;
}
关于计算范数/内积的C算法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/65648300/