unsigned int frcp_rn_core (unsigned int z)
{
unsigned int x, y;
int sign;
int expo;
sign = z & 0x80000000;
expo = (z >> 23) & 0xff;
x = expo - 1;
if (x > 0xfd) {
/* handle special cases */
x = z << 1;
/* zero or small denormal returns infinity of like sign */
if (x <= 0x00400000) {
return sign | 0x7f800000;
}
/* infinity returns zero of like sign */
else if (x == 0xff000000) {
return sign;
}
/* convert SNaNs to QNaNs */
else if (x > 0xff000000) {
return z | 0x00400000;
}
/* large denormal, normalize it */
else {
y = x < 0x00800000;
z = x << y;
expo = expo - y;
}
}
y = z & 0x007fffff;
#if USE_TABLE
z = 256 + rcp_tab[y >> 15]; /* approx */
#else
x = y << 3;
z = 0xe39ad7a0;
z = mul_32_hi (x, z) + 0x0154bde4;
z = mul_32_hi (x, z) + 0xfff87521;
z = mul_32_hi (x, z) + 0x00001ffa;
z = z >> 4;
#endif /* USE_TABLE */
y = y | 0x00800000;
/* cubically convergent approximation to reciprocal */
x = (unsigned int)-(int)(y * z); /* x = 1 - arg * approx */
x = mul_32_hi (x, x) + x; /* x = x * x + x */
z = z << 15;
z = mul_32_hi (x, z) + z; /* approx = x * approx + approx */
/* compute result exponent */
expo = 252 - expo;
if (expo >= 0) {
/* result is a normal */
x = y * z + y;
z = (expo << 23) + z;
z = z | sign;
/* round result */
if ((int)x <= (int)(y >> 1)) {
z++;
}
return z;
}
/* result is a denormal */
expo = -expo;
z = z >> expo;
x = y * z + y;
z = z | sign;
/* round result */
if ((int)x <= (int)(y >> 1)) {
z++;
}
return z;
}
/* 32-bit int, 64-bit long long int */
int mul_32_hi (int a, int b)
{
return (int)(unsigned int)(((unsigned long long)(((long long)a)*b)) >> 32);
}
基于表的变体使用的256项互易表可以构造如下：
static unsigned char rcp_tab[256];
for (int i = 0; i < 256; i++) {
rcp_tab[i] = (unsigned char)(((1./(1.+((i+.5)/256.)))*512.+.5)-256.);
}
3
取决于你想要多么复杂 . 保持相当简单，你可以尝试通过倒数近似除法 .
而不是计算：（n / d）你可以解决：n *（1 / d） .
要做到这一点，你需要使用某种方法计算倒数，例如，Newton-Raphson在完成最后的乘法步骤之前使用牛顿's method to calculate successively more accurate estimates of the reciprocal of the divisor until it' s "adequately"准确 .
2 回答
如果您的硬件具有足够快的整数乘数（比如4-5个周期），则使用迭代方法计算recip = 1 / divisor可能是最快的方法 . 然后，您将计算商数为dividend * recip . 如果硬件提供具有双倍宽度结果的整数乘法（即完整乘积）或提供两个32位整数乘积的高32位的“mulhi”指令，则对于必要的定点计算非常有用 .
您需要动态地重新定标定点计算，以在中间结果中保留接近32位，以获得精确到舍入后最终结果所需的24位的结果 .
最快的方法可能是产生9位起始近似值“约”1 /除数，然后是倒数的立方收敛迭代：
最简单的是预先计算起始近似并将其存储在256字节的数组中，由除数的位23:16（即尾数的8个最高有效小数位）索引 . 由于所有近似值都是0x100 ... 0x1FF范围内的数字（对应于0.5到0.998046875），因此不需要存储每个值的最高有效位，因为它将为“1” . 只需将0x100添加到检索到的表格元素中即可获得初始近似值的最终值 .
如果您无法承受256字节的表存储，则生成精确到9位的起始近似的另一种方法是3次多项式，其近似为1 /（1 f），其中f是除数尾数m的小数部分 . 由于IEEE-754已知m为[1.0,2.0]，因此f为[0.0,1.0] .
正确舍入到最近 - 或 - 偶数（如果需要）可以通过除数的初等商的反乘来 Build 余数，并选择最终商以使得余数最小化 .
下面的代码演示了上面讨论的近似原理，使用更简单的倒数情况，根据IEEE-754最近或偶数模式进行适当的舍入，并适当处理特殊情况（零，非正规，无穷大，NaN） . 它假定一个32位IEEE-754单精度浮点数可以从32位无符号整数位传输到32位无符号整数 . 然后，所有操作都在32位整数上执行 .
函数mul_32_high（）是一个占位符，用于特定于机器的操作，它返回两个32位位整数的有符号乘法的高32位 . 代替机器特定版本的半便携式实现是
基于表的变体使用的256项互易表可以构造如下：
取决于你想要多么复杂 . 保持相当简单，你可以尝试通过倒数近似除法 .
而不是计算：（n / d）你可以解决：n *（1 / d） .
要做到这一点，你需要使用某种方法计算倒数，例如，Newton-Raphson在完成最后的乘法步骤之前使用牛顿's method to calculate successively more accurate estimates of the reciprocal of the divisor until it' s "adequately"准确 .
编辑
刚看到你的更新 . 毕竟这可能对你有用，也可能不对你有用！