高速浮動小数点の逆数による2 ** 64 /除数の効率的な計算

Question

現在、さまざまな最新のプロセッサの高速単精度浮動小数点相互補完機能を使用して、64ビットの開始近似を計算する方法を検討しています固定小数点Newton-Raphson反復に基づく符号なし整数除算。可能な限り正確に2/divisorの計算が必要で、次の固定小数点反復の要件に基づいて、初期近似が数学的結果より小さくなければならない。これは、この計算が過小評価を提供する必要があることを意味します。私は現在、広範囲なテストに基づいて、うまく機能次のコードを、持っている：

#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 

このコードが機能しているが、それはほとんどのプラットフォーム上で正確に速くありません。機械固有のコードを少し必要とする明らかな改善の1つは、ハードウェアによって提供される高速浮動小数点の逆数を利用するコードで除算r = 1.0f/tを置き換えることです。これは、数学的結果の1 ulp以内の結果を生成するために反復処理で拡張することができるため、過小評価は既存のコードのコンテキストで生成されます。 x86_64のためのサンプル実装は次のようになります。nextafterf()の

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

実装は、一般的にパフォーマンスが最適化されていません。

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

をこれらのアプローチであると仮定すると、次のように組み込み関数float_as_int()とint_as_float()を通じて迅速IEEE int32およびその逆に754 binary32を再interpreteするための手段があるプラットフォームでは、我々はnextafterf()とスケーリングの使用を組み合わせることができます可能であれば、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ビットの固定小数点ニュートン・ラフソン・ステップでそれらをフォローする方が良いようです。

私の現在のアプローチの効率を改善する方法はありますか？特定のプラットフォーム用の組み込み関数を含むポータブルおよび準ポータブルな方法が関心があります（特に、現在支配的なCPUアーキテクチャであるx86およびARMの場合）。非常に高い最適化（/O3 /QxCORE-AVX2 /Qprec-div-）でIntelコンパイラを使用して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

このソリューションは、2つのアイデアを兼ね備え：

• あなたがいる限り、単に浮動小数点としてビットを再解釈し、定数を減算して浮動小数点に変換することができます数は特定の範囲内です。したがって、定数を追加し、再解釈し、その定数を減算します。これにより切り捨てられた結果が得られます（したがって、常に望ましい値以下です）。
• 指数と仮数の両方を無効にして逆数を近似することができます。これは、ビットをintとして解釈することによって達成される。

オプション1は特定の範囲でのみ動作するので、範囲を確認し、使用する定数を調整します。これは64ビットで動作します。なぜなら、必要なフロートは23ビットの精度しか持たないからです。

このコードの結果は2倍になりますが、浮動小数点への変換は簡単で、ハードウェアによってはビットまたは直接実行できます。

この後、ニュートンラフソン反復をしたいと思うでしょう。

このコードの多くは、単にマジックナンバーに変換されます。インテルコア7でこれをコンパイル

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

は、命令の数（分岐）を付与するものではありません、しかし、もちろん、何の乗算またはまったく分割します。 intとdoubleの間のキャストが速ければ、これはかなり早く動くはずです。