通过快速浮动有效计算2 ** 64 /除数

2018-06-30 16:30:32

我目前正在研究如何使用各种现代处理器的快速单精度浮点交互能力来计算基于定点牛顿 - 拉夫逊迭代的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()的实现通常不是性能优化的。在通过内在函数float_as_int()和int_as_float()快速将IEEE 754 binary32重新解释为int32 （反之亦然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位定点牛顿 - 拉夫逊步骤进行跟踪。

有没有方法提高我目前的方法的效率？ 涉及用于特定平台的内部函数的便携式和半便携式方法将是有趣的（特别是对于x86和ARM作为当前占主导的CPU架构）。使用英特尔编译器以非常高的优化（ /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;nt"  // rax = a
        "mulq  %2;nt"         // rdx:rax = a * b
        "movq  %%rdx, %0;nt"  // res = (a * b)<63:32>
        : "=rm" (res)
        : "rm"(a), "rm"(b)
        : "%rax", "%rdx");
    return res;
}

该解决方案结合了两个想法

只要将数字重新解释为浮点数并减去一个常数，只要该数字在特定范围内，就可以转换为浮点数。因此，添加一个常数，重新解释，然后减去该常数。这会给出一个截断结果（因此总是小于或等于期望值）。

您可以通过否定指数和尾数来近似倒数。这可以通过将位解释为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;                                                
}

在英特尔核心7上编译时会给出一些指令（和一个分支），但是，当然，根本不会有乘法或除法。如果int和double之间的转换速度很快，这应该运行得非常快。

我怀疑float（只有23位精度）需要超过1或2个Newton-Raphson迭代才能达到你想要的精度，但是我没有做数学......

链接地址: http://www.djcxy.com/p/85607.html

上一篇: Efficient computation of 2**64 / divisor via fast floating

下一篇: Vectorizable implementation of complementary error function erfcf()