online - Cálculo eficiente de 2** 64/divisor mediante recíproco de punto flotante rápido

punto flotante ieee 754 ejemplos

Question

Actualmente estoy buscando formas de utilizar la capacidad recíproca de coma flotante de precisión simple de varios procesadores modernos para calcular una aproximación inicial para una división de enteros sin signo de 64 bits basada en iteraciones Newton-Raphson de punto fijo. Requiere un cálculo de 2 ⁶⁴ / divisor, con la mayor precisión posible, donde la aproximación inicial debe ser menor que, o igual a, el resultado matemático, en función de los requisitos de las siguientes iteraciones de punto fijo. Esto significa que este cálculo debe proporcionar una subestimación. Actualmente tengo el siguiente código, que funciona bien, basado en pruebas exhaustivas:

#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 

Si bien este código es funcional, no es exactamente rápido en la mayoría de las plataformas. Una mejora obvia, que requiere un poco de código específico de la máquina, es reemplazar la división r = 1.0f / t con un código que hace uso de un recíproco de punto flotante rápido proporcionado por el hardware. Esto puede ser aumentado con iteración para producir un resultado que está dentro de 1 ulp del resultado matemático, por lo que se produce una subestimación en el contexto del código existente. Una implementación de muestra para x86_64 sería:

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

Las implementaciones de nextafterf() normalmente no están optimizadas para el rendimiento. En las plataformas donde hay medios para reinterpretar rápidamente un IEEE 754 binary32 en un int32 y viceversa, a través de los elementos intrínsecos float_as_int() e int_as_float() , podemos combinar el uso de nextafterf() y la escala de la siguiente manera:

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

Suponiendo que estos enfoques son posibles en una plataforma determinada, esto nos deja con las conversiones entre float y uint64_t como los principales obstáculos. La mayoría de las plataformas no proporcionan una instrucción que realiza una conversión desde uint64_t para float con el modo de redondeo estático (aquí: hacia infinito positivo = arriba), y algunas no ofrecen instrucciones para convertir entre uint64_t y tipos de punto flotante, lo que hace un cuello de botella de rendimiento.

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 */

Una implementación portátil, pero lenta, de uint64_to_float_ru usa cambios dinámicos en el modo de redondeo 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;
}

He analizado varios enfoques de división y de mezcla de bits para tratar las conversiones (por ejemplo, redondeo en el lado entero, luego uso una conversión normal a float que usa el modo de redondeo IEEE 754 redondo-a-más cercano-o-incluso) , pero la sobrecarga que esto crea hace que este cálculo sea poco atractivo recíproco de punto flotante desde una perspectiva de rendimiento. Tal como está, parece que sería mejor generar una aproximación inicial usando un LUT clásico con interpolación, o una aproximación polinómica de punto fijo, y seguir con un paso Newton-Raphson de punto fijo de 32 bits.

¿Hay formas de mejorar la eficiencia de mi enfoque actual? Las formas portátiles y semi-portátiles que implican intrínsecos para plataformas específicas serían de interés (en particular para x86 y ARM como las arquitecturas de CPU actualmente dominantes). Compilando para x86_64 usando el compilador de Intel a una optimización muy alta ( /O3 /QxCORE-AVX2 /Qprec-div- ) el cálculo de la aproximación inicial toma más instrucciones que la iteración, que toma alrededor de 20 instrucciones. A continuación se muestra el código de división completo para referencia, que muestra la aproximación en contexto.

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() generalmente se correlacionaría con un código intrínseco específico de plataforma, o un poco de código ensamblador en línea. En x86_64 actualmente uso esta implementación:

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

Esta solución combina dos ideas:

• Puede convertir a coma flotante simplemente reinterpretando los bits como punto flotante y restando una constante, siempre que el número esté dentro de un rango particular. Así que agrega una constante, reinterpreta y luego resta esa constante. Esto dará un resultado truncado (que por lo tanto siempre es menor o igual que el valor deseado).
• Puedes aproximar el recíproco al negar tanto el exponente como la mantisa. Esto puede lograrse interpretando los bits como int.

La opción 1 aquí solo funciona en un rango determinado, por lo que verificamos el rango y ajustamos las constantes utilizadas. Esto funciona en 64 bits porque el flotador deseado solo tiene 23 bits de precisión.

El resultado en este código será el doble, pero la conversión a float es trivial, y se puede hacer en los bits o directamente, dependiendo del hardware.

Después de esto, querrás hacer la (s) iteración (es) de Newton-Raphson.

Gran parte de este código simplemente se convierte en números mágicos.

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

Compilar esto en un núcleo Intel 7 da una serie de instrucciones (y una rama), pero, por supuesto, no se multiplica ni se divide en absoluto. Si los lanzamientos entre int y double son rápidos, esto debería ejecutarse bastante rápido.

Sospecho que la flotación (con solo 23 bits de precisión) requerirá más de 1 o 2 iteraciones de Newton-Raphson para obtener la precisión que desea, pero no he hecho la matemática ...