c - 通过快速浮点倒数高效计算 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() 和 将 IEEE 754 binary32 快速重新解释为 int32 的平台上，反之亦然code>int_as_float()，我们可以结合使用 nextafterf() 和缩放如下:

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

假设这些方法在给定平台上可行，这给我们留下了 float 和 uint64_t 之间的转换作为主要障碍。大多数平台不提供执行从 uint64_t 到 float 静态舍入模式(此处:朝向正无穷 = 向上)的转换的指令，有些平台不提供在 uint64_t 和浮点类型之间转换的任何指令，使这成为性能瓶颈。

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

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

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

此处的选项 1 仅在一定范围内有效，因此我们检查范围并调整使用的常量。这适用于 64 位，因为所需的 float 只有 23 位精度。

此代码中的结果将是 double 的，但转换为 float 是微不足道的，可以在位上完成或直接完成，具体取决于硬件。

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

大部分代码只是简单地转换为魔数(Magic Number)。

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

在英特尔酷睿 7 上编译它会给出许多指令(和一个分支)，但是，当然，根本没有乘法或除法。如果 int 和 double 之间的转换速度很快，这应该会很快运行。

我怀疑 float (只有 23 位精度)需要超过 1 或 2 次牛顿-拉夫森迭代才能获得您想要的精度，但我还没有计算...