通过快速浮点倒数高效计算 2**64 / 除数

Question

我目前正在研究如何使用各种现代处理器的快速单精度浮点倒数功能来计算基于定点 Newton-Raphson 迭代的 64 位无符号整数除法的起始近似值。它需要尽可能准确地计算 2⁶⁴ / 除数，其中初始近似值必须小于或等于数学结果，基于以下定点迭代的要求。这意味着该计算需要提供低估。根据广泛的测试，我目前有以下代码，效果很好：

#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() 的实现通常没有进行性能优化。在可以通过内在函数 float_as_int() 和 int_as_float() 快速将 IEEE 754 binary32 重新解释为 int32 的平台上，我们可以结合使用 nextafterf() 和缩放如下：

s = int_as_float (float_as_int (r) + 0x1fffffff);

假设这些方法在给定的平台上是可行的，这给我们留下了float 和uint64_t 之间的转换作为主要障碍。大多数平台不提供使用静态舍入模式执行从uint64_t 到float 转换的指令（这里：朝向正无穷大=向上），并且有些平台不提供在uint64_t 和浮动之间转换的任何指令-point 类型，使其成为性能瓶颈。

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

我研究了各种拆分和位旋转方法来处理转换（例如，在整数端进行舍入，然后使用正常转换为 float，它使用 IEEE 754 舍入模式舍入到最近- 或 - 甚至），但是从性能的角度来看，这产生的开销使得通过快速浮点倒数进行的计算没有吸引力。就目前而言，看起来我最好通过使用带有插值的经典 LUT 或定点多项式逼近来生成起始逼近，然后使用 32 位定点 Newton-Raphson 步长。

有没有办法提高我当前方法的效率？ 涉及特定平台内在函数的可移植和半可移植方法会很有趣（特别是对于 x86 和 ARM 作为当前占主导地位的 CPU 架构)。使用 Intel 编译器以非常高的优化 (/O3 /QxCORE-AVX2 /Qprec-div-) 为 x86_64 编译，初始近似的计算比迭代需要更多的指令，迭代需要大约 20 条指令。下面是完整的除法代码供参考，在上下文中显示了近似值。

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

Answer 1

这个解决方案结合了两个想法：

• 只要数字在特定范围内，只需将位重新解释为浮点数并减去一个常数，即可转换为浮点数。所以添加一个常数，重新解释，然后减去该常数。这将给出截断的结果（因此始终小于或等于所需值）。
• 您可以通过取反指数和尾数来近似倒数。这可以通过将位解释为 int 来实现。

这里的选项1只在一定范围内有效，所以我们检查范围并调整使用的常数。这适用于 64 位，因为所需的浮点数只有 23 位精度。

此代码中的结果将是双精度，但转换为浮点数是微不足道的，可以在位上完成，也可以直接完成，具体取决于硬件。

在此之后，您需要进行 Newton-Raphson 迭代。

大部分代码只是简单地转换为幻数。

double                                                       
u64tod_inv( uint64_t u64 ) {                                 
  __asm__( "#annot0" );                                      
  union {                                                    
    double f;                                                
    struct {                                                 
      unsigned long m:52; // careful here with endianess     
      unsigned long x:11;                                    
      unsigned long s:1;                                     
    } u64;                                                   
    uint64_t u64i;                                           
  } z,                                                       
        magic0 = { .u64 = { 0, (1<<10)-1 + 52, 0 } },        
        magic1 = { .u64 = { 0, (1<<10)-1 + (52+12), 0 } },   
        magic2 = { .u64 = { 0, 2046, 0 } };                  

  __asm__( "#annot1" );                                      
  if( u64 < (1UL << 52UL ) ) {                               
    z.u64i = u64 + magic0.u64i;                              
    z.f   -= magic0.f;                                       
  } else {                                                   
    z.u64i = ( u64 >> 12 ) + magic1.u64i;                    
    z.f   -= magic1.f;                                       
  }                                                          
  __asm__( "#annot2" );                                      

  z.u64i = magic2.u64i - z.u64i;                             

  return z.f;                                                
}                                                            

在 Intel core 7 上编译它会给出许多指令（和一个分支），但当然，根本没有乘法或除法。如果 int 和 double 之间的转换很快，这应该运行得很快。

我怀疑浮点数（只有 23 位精度）将需要超过 1 或 2 次 Newton-Raphson 迭代才能获得您想要的精度，但我还没有计算过...