通过快速浮点倒数有效计算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()通常不是性能优化的.在那里有手段迅速重新interprete一个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和浮点类型之间的任何指令,这使其成为性能瓶颈.

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架构).编译使用英特尔编译器在非常高的优化x86_64的(/O3 /QxCORE-AVX2 /Qprec-div-)的初始近似的计算时间比迭代,大约需要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 位精度。

此代码中的结果将是双精度的，但转换为浮点值很简单，并且可以按位或直接完成，具体取决于硬件。

之后，您需要进行牛顿-拉夫森迭代。

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

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 之间的转换很快，那么它应该运行得很快。

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