我目前正在探索使用现代处理器的快速单精度浮点数倒数能力来计算基于定点牛顿-拉夫逊迭代的64位无符号整数除法的起始近似值。它需要尽可能准确地计算264 / 除数，其中初始近似值必须小于或等于数学结果，基于以下定点迭代的要求。这意味着此计算需要提供低估。根据广泛测试获得以下代码可正常工作：
#include <stdint.h> // import uint64_t
#include <math.h> // import nextafterf()
uint64_t divisor, recip;
float r, s, t;
t = uint64_to_float_ru (divisor); // ensure t >= divisor
r = 1.0f / t;
s = 0x1.0p64f * nextafterf (r, 0.0f);
recip = (uint64_t)s; // underestimate of 2**64 / divisor
虽然这段代码是可行的，但在大多数平台上并不是很快。一种明显的改进方法，需要一些机器特定的代码，就是用硬件提供的快速浮点倒数替换除法r = 1.0f / t。可以通过迭代来增强，以产生一个与数学结果相差不超过1 ulp的结果，在现有代码的情况下产生一个低估值。对于x86_64的示例实现如下：
#include <xmmintrin.h>
/* Compute 1.0f/a almost correctly rounded. Halley iteration with cubic convergence */
inline float fast_recip_f32 (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (r, -a, 1.0f);
e = fmaf (e, e, e);
r = fmaf (e, r, r);
return r;
}
nextafterf() 的实现通常不是针对性能进行优化的。在存在一种快速重新解释 IEEE 754 binary32 为 int32，反之亦然的平台上，通过使用内置函数 float_as_int() 和 int_as_float()，我们可以结合使用 nextafterf() 和缩放，具体如下：
s = int_as_float (float_as_int (r) + 0x1fffffff);
假设在特定平台上这些方法是可行的，那么在
float和uint64_t之间的转换将成为主要障碍。大多数平台不提供执行从uint64_t到float的转换的指令，并且一些平台不提供任何指令来转换uint64_t和浮点类型之间的转换，这成为了性能瓶颈。请注意保留HTML标签。t = uint64_to_float_ru (divisor);
r = fast_recip_f32 (t);
s = int_as_float (float_as_int (r) + 0x1fffffff);
recip = (uint64_t)s; /* underestimate of 2**64 / divisor */
一种可移植但较慢的 uint64_to_float_ru 实现会使用动态更改 FPU 舍入模式：
#include <fenv.h>
#pragma STDC FENV_ACCESS ON
float uint64_to_float_ru (uint64_t a)
{
float res;
int curr_mode = fegetround ();
fesetround (FE_UPWARD);
res = (float)a;
fesetround (curr_mode);
return res;
}
我研究了各种分割和位操作的方法来处理转换（例如，在整数侧进行四舍五入，然后使用使用IEEE 754舍入模式"Round-to-Nearest-or-Even"的普通float转换），但是这会带来额外的开销，从性能角度考虑，通过快速浮点数倒数计算并不具有吸引力。目前看来，最好的方法是使用经典的LUT和插值生成起始近似值，或者使用定点多项式逼近方法，并配合32位定点牛顿迭代步骤。
是否有提高当前方法效率的方法？ 可移植和半可移植的方式包括特定平台的intrinsics将会是感兴趣的（特别是对于x86和ARM作为当前主导的CPU架构）。在使用Intel编译器以非常高的优化级别（/O3 /QxCORE-AVX2 /Qprec-div-）编译为x86_64时，初始逼近值的计算所需指令比迭代所需指令多。以下是完整的除法代码，显示了上下文中的逼近值。
uint64_t udiv64 (uint64_t dividend, uint64_t divisor)
{
uint64_t temp, quot, rem, recip, neg_divisor = 0ULL - divisor;
float r, s, t;
/* compute initial approximation for reciprocal; must be underestimate! */
t = uint64_to_float_ru (divisor);
r = 1.0f / t;
s = 0x1.0p64f * nextafterf (r, 0.0f);
recip = (uint64_t)s; /* underestimate of 2**64 / divisor */
/* perform Halley iteration with cubic convergence to refine reciprocal */
temp = neg_divisor * recip;
temp = umul64hi (temp, temp) + temp;
recip = umul64hi (recip, temp) + recip;
/* compute preliminary quotient and remainder */
quot = umul64hi (dividend, recip);
rem = dividend - divisor * quot;
/* adjust quotient if too small; quotient off by 2 at most */
if (rem >= divisor) quot += ((rem - divisor) >= divisor) ? 2 : 1;
/* handle division by zero */
if (divisor == 0ULL) quot = ~0ULL;
return quot;
}
umul64hi() 通常会映射到平台特定的内置函数或一些内联汇编代码。在 x86_64 上，我目前使用以下实现：
inline uint64_t umul64hi (uint64_t a, uint64_t b)
{
uint64_t res;
__asm__ (
"movq %1, %%rax;\n\t" // rax = a
"mulq %2;\n\t" // rdx:rax = a * b
"movq %%rdx, %0;\n\t" // res = (a * b)<63:32>
: "=rm" (res)
: "rm"(a), "rm"(b)
: "%rax", "%rdx");
return res;
}