Существует ли функция C ++, которая возвращает точно значение встроенной операции ЦП RSQRTSS для обратного квадратного корня? - PullRequest
Новая
       33

Существует ли функция C ++, которая возвращает точно значение встроенной операции ЦП RSQRTSS для обратного квадратного корня?

1 голос
/ 29 октября 2019

Я ищу функцию C ++, которая возвращает обратный sqrt числа с плавающей запятой: rsqrt(x) = 1/sqrt(x) с использованием точного метода, такого как встроенная операция XMM RSQRTSS (ср. https://www.felixcloutier.com/x86/rsqrtss). (т. Е. Яхочу встроенное приближение, а не более точное 1/sqrtf, и мне наплевать на скорость (много).)

По этому вопросу:

есть быстрая стандартная библиотечная функция C или C ++ для обратного квадратного корня двойной точности?

... по крайней мере, не существует "быстрого способа с двойной точностью", который бы выполнял это с помощью стандартной библиотеки C ++. как это сделать медленно, нестандартно и с float с?

Ответы [ 3 ]

3 голосов
/ 29 октября 2019

В стандартной библиотеке нет функции, которая бы делала это, но ваш компилятор мог бы оптимизировать выражение 1 / sqrt(value) таким образом, чтобы он выдавал инструкцию RSQRTSS. Например, с флагами компилятора -ffast-math -march=native, GCC будет выдавать эту инструкцию, см .: https://godbolt.org/z/cL6seG

1 голос
/ 31 октября 2019

Инструкция RSQRTSS легко доступна через _mm_rsqrt_ss(), объявленный в immintrin.h. Но мы также можем эмулировать инструкцию в программном обеспечении, как это делается в функции my_rsqrtf() ниже. Просто наблюдая за выводом RSQRTSS, легко обнаружить, что значения его функций основаны на (виртуальной) таблице из 2 11 записей, каждая размером 12 бит.

Обратите внимание на «виртуальный» атрибут, поскольку маловероятно, что в оборудовании используется прямая таблица 24 Кбит. Мой анализ шаблонов в записях таблицы не предполагает использования двудольных таблиц . Можно использовать гораздо более простую схему сжатия - например, ту, которую я использовал в приведенном ниже коде - на основе таблицы базовых значений и таблицы смещений . Моя схема требует только узкого сумматора, но уменьшает объем ПЗУ до 13 Кбит, то есть почти вдвое.

Приведенная ниже реализация была разработана и протестирована на процессоре Intel Xeon E3-1270 V2, использующем архитектуру Ivy Bridge. Могут быть некоторые функциональные различия в реализации RSQRTSS между различными архитектурами Intel, и такие различия, вероятно, существуют между архитектурами разных поставщиков x86-84.

Приведенная ниже структура проверяет, обеспечивает ли эмуляция my_rsqrtf()побитовые идентичные результаты равны RSQRTSS для всех четырех режимов округления, двух режимов DAZ (денормали равны нулю) и двух режимов FTZ (сброса в ноль). Мы обнаруживаем, что на результаты функции не влияют никакие режимы, что соответствует тому, как Intel указала RSQRTSS в Руководстве разработчика программного обеспечения Intel® 64 и IA-32 Architectures :

На команду RSQRTSS не влияют управляющие биты округления в регистре MXCSR. Когда исходное значение равно 0.0, возвращается ∞ знака исходного значения. Денормальное значение источника рассматривается как 0.0 (с тем же знаком). Когда исходное значение является отрицательным значением (кроме -0.0), возвращается неопределенное число с плавающей точкой. Когда исходным значением является SNaN или QNaN, SNaN преобразуется в QNaN или возвращается исходный QNaN. 

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include "immintrin.h"

/* SSE reference for RSQRTSS instruction */
float sse_rsqrtf (float a, uint32_t daz, uint32_t ftz, uint32_t rnd)
{
    __m128 b, t;
    float res;
    uint32_t old_mxcsr;
    old_mxcsr = _mm_getcsr();
    _MM_SET_DENORMALS_ZERO_MODE (daz);
    _MM_SET_FLUSH_ZERO_MODE (ftz);
    _MM_SET_ROUNDING_MODE (rnd);
    b = _mm_set_ss (a);
    t = _mm_rsqrt_ss (b);
    _mm_store_ss (&res, t);
    _mm_setcsr (old_mxcsr);
    return res;
}

inline uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

inline float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

#define LOG2_NBR_TAB_ENTRIES (11)
#define NBR_TAB_ENTRIES      (1 << LOG2_NBR_TAB_ENTRIES)
#define TAB_ENTRY_BITS       (12)
#define BASE_TAB_ENTRY_BITS  (9)

/* 128 9-bit entries = 1152 bits */
const uint16_t base_tab[128] = {
    0x0ce, 0x0c9, 0x0c3, 0x0be, 0x0b9, 0x0b4, 0x0af, 0x0aa,
    0x0a6, 0x0a1, 0x09d, 0x098, 0x094, 0x090, 0x08b, 0x087,
    0x083, 0x07f, 0x07c, 0x078, 0x074, 0x070, 0x06d, 0x069,
    0x066, 0x062, 0x05f, 0x05c, 0x058, 0x055, 0x052, 0x04f,
    0x04c, 0x049, 0x046, 0x043, 0x040, 0x03d, 0x03a, 0x038,
    0x035, 0x032, 0x030, 0x02d, 0x02a, 0x028, 0x025, 0x023,
    0x021, 0x01e, 0x01c, 0x019, 0x017, 0x015, 0x013, 0x010,
    0x00e, 0x00c, 0x00a, 0x008, 0x006, 0x004, 0x002, 0x000,
    0x1f8, 0x1f0, 0x1e9, 0x1e1, 0x1da, 0x1d3, 0x1cc, 0x1c5,
    0x1bf, 0x1b8, 0x1b2, 0x1ab, 0x1a5, 0x19f, 0x199, 0x194,
    0x18e, 0x188, 0x183, 0x17e, 0x178, 0x173, 0x16e, 0x169,
    0x164, 0x15f, 0x15a, 0x156, 0x151, 0x14d, 0x148, 0x144,
    0x13f, 0x13b, 0x137, 0x133, 0x12f, 0x12b, 0x127, 0x123,
    0x11f, 0x11b, 0x118, 0x114, 0x110, 0x10d, 0x109, 0x106,
    0x102, 0x0ff, 0x0fc, 0x0f8, 0x0f5, 0x0f2, 0x0ef, 0x0eb,
    0x0e8, 0x0e5, 0x0e2, 0x0df, 0x0dc, 0x0d9, 0x0d7, 0x0d4,
};

/* 2048 6-bit entries = 12288 bits */
const uint8_t ofs_tab[2048] = {
    0x2f, 0x2c, 0x2a, 0x27, 0x24, 0x21, 0x1e, 0x1c, 0x19, 0x16, 0x13, 0x10, 0x0e, 0x0b, 0x08, 0x05,
    0x2b, 0x28, 0x25, 0x22, 0x1f, 0x1d, 0x1a, 0x17, 0x15, 0x12, 0x0f, 0x0c, 0x0a, 0x07, 0x04, 0x02,
    0x2f, 0x2c, 0x29, 0x27, 0x24, 0x21, 0x1f, 0x1c, 0x19, 0x17, 0x14, 0x11, 0x0f, 0x0c, 0x09, 0x07,
    0x2c, 0x29, 0x27, 0x24, 0x22, 0x1f, 0x1c, 0x1a, 0x17, 0x15, 0x12, 0x0f, 0x0d, 0x0a, 0x08, 0x05,
    0x2a, 0x28, 0x25, 0x23, 0x20, 0x1e, 0x1b, 0x18, 0x16, 0x13, 0x11, 0x0e, 0x0c, 0x09, 0x07, 0x04,
    0x2a, 0x27, 0x24, 0x22, 0x1f, 0x1d, 0x1a, 0x18, 0x15, 0x13, 0x10, 0x0e, 0x0b, 0x09, 0x07, 0x04,
    0x2a, 0x27, 0x25, 0x22, 0x20, 0x1d, 0x1b, 0x18, 0x16, 0x13, 0x11, 0x0f, 0x0c, 0x0a, 0x07, 0x05,
    0x2a, 0x28, 0x26, 0x23, 0x21, 0x1e, 0x1c, 0x1a, 0x17, 0x15, 0x12, 0x10, 0x0e, 0x0b, 0x09, 0x07,
    0x24, 0x22, 0x1f, 0x1d, 0x1b, 0x18, 0x16, 0x14, 0x11, 0x0f, 0x0d, 0x0a, 0x08, 0x06, 0x03, 0x01,
    0x27, 0x24, 0x22, 0x20, 0x1d, 0x1b, 0x19, 0x16, 0x14, 0x12, 0x10, 0x0d, 0x0b, 0x09, 0x06, 0x04,
    0x22, 0x20, 0x1d, 0x1b, 0x19, 0x17, 0x14, 0x12, 0x10, 0x0e, 0x0b, 0x09, 0x07, 0x05, 0x02, 0x00,
    0x26, 0x24, 0x21, 0x1f, 0x1d, 0x1b, 0x19, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0b, 0x09, 0x07, 0x05,
    0x23, 0x20, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x09, 0x06, 0x04, 0x02,
    0x20, 0x1e, 0x1c, 0x1a, 0x17, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x09, 0x06, 0x04, 0x02, 0x00,
    0x26, 0x24, 0x22, 0x20, 0x1e, 0x1c, 0x19, 0x17, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x09, 0x07,
    0x25, 0x23, 0x21, 0x1f, 0x1d, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06,
    0x24, 0x22, 0x20, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06,
    0x24, 0x22, 0x20, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06,
    0x1d, 0x1b, 0x19, 0x17, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x09, 0x07, 0x05, 0x03, 0x01, 0x00,
    0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0b, 0x09, 0x07, 0x05, 0x03, 0x01,
    0x1f, 0x1d, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0d, 0x0b, 0x09, 0x07, 0x05, 0x03,
    0x21, 0x20, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x0a, 0x08, 0x06,
    0x1c, 0x1a, 0x19, 0x17, 0x15, 0x13, 0x11, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x07, 0x05, 0x03, 0x01,
    0x1f, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x15, 0x13, 0x11, 0x0f, 0x0e, 0x0c, 0x0a, 0x08, 0x07, 0x05,
    0x1b, 0x19, 0x18, 0x16, 0x14, 0x12, 0x11, 0x0f, 0x0d, 0x0b, 0x0a, 0x08, 0x06, 0x04, 0x03, 0x01,
    0x1f, 0x1e, 0x1c, 0x1a, 0x18, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0e, 0x0c, 0x0b, 0x09, 0x07, 0x06,
    0x1c, 0x1a, 0x19, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0e, 0x0d, 0x0b, 0x09, 0x08, 0x06, 0x04, 0x03,
    0x19, 0x17, 0x16, 0x14, 0x12, 0x11, 0x0f, 0x0d, 0x0c, 0x0a, 0x08, 0x07, 0x05, 0x03, 0x02, 0x00,
    0x1f, 0x1d, 0x1b, 0x1a, 0x18, 0x16, 0x15, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0b, 0x09, 0x08, 0x06,
    0x1d, 0x1b, 0x19, 0x18, 0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0c, 0x0b, 0x09, 0x08, 0x06, 0x04,
    0x1b, 0x19, 0x18, 0x16, 0x15, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0b, 0x0a, 0x08, 0x06, 0x05, 0x03,
    0x1a, 0x18, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0f, 0x0d, 0x0c, 0x0a, 0x09, 0x07, 0x06, 0x04, 0x02,
    0x19, 0x17, 0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0b, 0x0a, 0x08, 0x07, 0x05, 0x03, 0x02,
    0x18, 0x17, 0x15, 0x14, 0x12, 0x11, 0x0f, 0x0e, 0x0c, 0x0b, 0x09, 0x08, 0x06, 0x05, 0x03, 0x02,
    0x18, 0x17, 0x15, 0x14, 0x12, 0x11, 0x0f, 0x0e, 0x0d, 0x0b, 0x0a, 0x08, 0x07, 0x05, 0x04, 0x02,
    0x19, 0x17, 0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0c, 0x0a, 0x09, 0x07, 0x06, 0x04, 0x03,
    0x19, 0x18, 0x16, 0x15, 0x14, 0x12, 0x11, 0x0f, 0x0e, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x05, 0x04,
    0x1a, 0x19, 0x18, 0x16, 0x15, 0x13, 0x12, 0x10, 0x0f, 0x0e, 0x0c, 0x0b, 0x09, 0x08, 0x07, 0x05,
    0x1c, 0x1a, 0x19, 0x18, 0x16, 0x15, 0x13, 0x12, 0x11, 0x0f, 0x0e, 0x0c, 0x0b, 0x0a, 0x08, 0x07,
    0x15, 0x14, 0x13, 0x11, 0x10, 0x0f, 0x0d, 0x0c, 0x0a, 0x09, 0x08, 0x06, 0x05, 0x04, 0x02, 0x01,
    0x17, 0x16, 0x15, 0x13, 0x12, 0x11, 0x0f, 0x0e, 0x0d, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x04, 0x03,
    0x1a, 0x18, 0x17, 0x16, 0x14, 0x13, 0x12, 0x10, 0x0f, 0x0e, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x06,
    0x14, 0x13, 0x12, 0x10, 0x0f, 0x0e, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x05, 0x03, 0x02, 0x01,
    0x17, 0x16, 0x15, 0x13, 0x12, 0x11, 0x0f, 0x0e, 0x0d, 0x0c, 0x0a, 0x09, 0x08, 0x06, 0x05, 0x04,
    0x1b, 0x19, 0x18, 0x17, 0x15, 0x14, 0x13, 0x12, 0x10, 0x0f, 0x0e, 0x0c, 0x0b, 0x0a, 0x09, 0x07,
    0x16, 0x15, 0x13, 0x12, 0x11, 0x10, 0x0e, 0x0d, 0x0c, 0x0b, 0x09, 0x08, 0x07, 0x06, 0x04, 0x03,
    0x1a, 0x19, 0x17, 0x16, 0x15, 0x14, 0x12, 0x11, 0x10, 0x0f, 0x0d, 0x0c, 0x0b, 0x0a, 0x08, 0x07,
    0x16, 0x15, 0x13, 0x12, 0x11, 0x10, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x05, 0x03,
    0x12, 0x11, 0x10, 0x0f, 0x0d, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x05, 0x04, 0x02, 0x01, 0x00,
    0x17, 0x16, 0x14, 0x13, 0x12, 0x11, 0x10, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x05,
    0x14, 0x12, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0b, 0x0a, 0x09, 0x08, 0x07, 0x05, 0x04, 0x03, 0x02,
    0x19, 0x18, 0x16, 0x15, 0x14, 0x13, 0x12, 0x11, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x08, 0x07,
    0x16, 0x15, 0x14, 0x13, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x09, 0x08, 0x07, 0x06, 0x05,
    0x14, 0x13, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03,
    0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02, 0x01,
    0x18, 0x16, 0x15, 0x14, 0x13, 0x12, 0x11, 0x10, 0x0f, 0x0e, 0x0c, 0x0b, 0x0a, 0x09, 0x08, 0x07,
    0x16, 0x15, 0x14, 0x12, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x09, 0x08, 0x06, 0x05,
    0x14, 0x13, 0x12, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04,
    0x13, 0x12, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0b, 0x0a, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03,
    0x12, 0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0a, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02,
    0x11, 0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x09, 0x08, 0x07, 0x06, 0x04, 0x03, 0x02, 0x01,
    0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02, 0x01,
    0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02, 0x01,
    0x10, 0x0f, 0x0e, 0x0d, 0x0c, 0x0b, 0x0a, 0x09, 0x08, 0x07, 0x06, 0x05, 0x04, 0x03, 0x02, 0x01,
    0x3e, 0x3a, 0x36, 0x32, 0x2e, 0x2a, 0x26, 0x22, 0x1e, 0x1a, 0x16, 0x12, 0x0e, 0x0b, 0x07, 0x03,
    0x3f, 0x3b, 0x37, 0x33, 0x2f, 0x2b, 0x27, 0x24, 0x20, 0x1c, 0x18, 0x14, 0x10, 0x0c, 0x09, 0x05,
    0x39, 0x35, 0x31, 0x2e, 0x2a, 0x26, 0x22, 0x1e, 0x1b, 0x17, 0x13, 0x0f, 0x0c, 0x08, 0x04, 0x00,
    0x3d, 0x39, 0x35, 0x31, 0x2e, 0x2a, 0x26, 0x23, 0x1f, 0x1b, 0x18, 0x14, 0x10, 0x0d, 0x09, 0x05,
    0x3a, 0x36, 0x32, 0x2f, 0x2b, 0x27, 0x24, 0x20, 0x1d, 0x19, 0x15, 0x12, 0x0e, 0x0b, 0x07, 0x03,
    0x38, 0x34, 0x31, 0x2d, 0x2a, 0x26, 0x22, 0x1f, 0x1b, 0x18, 0x14, 0x11, 0x0d, 0x0a, 0x06, 0x03,
    0x37, 0x34, 0x30, 0x2d, 0x29, 0x26, 0x22, 0x1f, 0x1b, 0x18, 0x15, 0x11, 0x0e, 0x0a, 0x07, 0x03,
    0x38, 0x35, 0x31, 0x2e, 0x2a, 0x27, 0x24, 0x20, 0x1d, 0x19, 0x16, 0x13, 0x0f, 0x0c, 0x09, 0x05,
    0x32, 0x2e, 0x2b, 0x28, 0x24, 0x21, 0x1e, 0x1a, 0x17, 0x14, 0x11, 0x0d, 0x0a, 0x07, 0x03, 0x00,
    0x35, 0x31, 0x2e, 0x2b, 0x28, 0x24, 0x21, 0x1e, 0x1b, 0x17, 0x14, 0x11, 0x0e, 0x0a, 0x07, 0x04,
    0x31, 0x2e, 0x2a, 0x27, 0x24, 0x21, 0x1e, 0x1a, 0x17, 0x14, 0x11, 0x0e, 0x0b, 0x07, 0x04, 0x01,
    0x36, 0x33, 0x30, 0x2c, 0x29, 0x26, 0x23, 0x20, 0x1d, 0x1a, 0x17, 0x13, 0x10, 0x0d, 0x0a, 0x07,
    0x34, 0x31, 0x2e, 0x2b, 0x28, 0x25, 0x21, 0x1e, 0x1b, 0x18, 0x15, 0x12, 0x0f, 0x0c, 0x09, 0x06,
    0x33, 0x30, 0x2d, 0x2a, 0x27, 0x24, 0x21, 0x1e, 0x1b, 0x18, 0x15, 0x12, 0x0f, 0x0c, 0x09, 0x06,
    0x33, 0x30, 0x2d, 0x2a, 0x27, 0x24, 0x21, 0x1e, 0x1b, 0x18, 0x15, 0x13, 0x10, 0x0d, 0x0a, 0x07,
    0x2c, 0x29, 0x26, 0x23, 0x20, 0x1d, 0x1a, 0x18, 0x15, 0x12, 0x0f, 0x0c, 0x09, 0x06, 0x03, 0x01,
    0x2e, 0x2b, 0x28, 0x25, 0x22, 0x1f, 0x1d, 0x1a, 0x17, 0x14, 0x11, 0x0e, 0x0c, 0x09, 0x06, 0x03,
    0x30, 0x2e, 0x2b, 0x28, 0x25, 0x22, 0x20, 0x1d, 0x1a, 0x17, 0x14, 0x12, 0x0f, 0x0c, 0x09, 0x07,
    0x2c, 0x29, 0x26, 0x24, 0x21, 0x1e, 0x1b, 0x19, 0x16, 0x13, 0x10, 0x0e, 0x0b, 0x08, 0x06, 0x03,
    0x28, 0x25, 0x23, 0x20, 0x1d, 0x1b, 0x18, 0x15, 0x13, 0x10, 0x0d, 0x0b, 0x08, 0x05, 0x03, 0x00,
    0x2d, 0x2b, 0x28, 0x25, 0x23, 0x20, 0x1d, 0x1b, 0x18, 0x15, 0x13, 0x10, 0x0e, 0x0b, 0x08, 0x06,
    0x2b, 0x28, 0x26, 0x23, 0x21, 0x1e, 0x1b, 0x19, 0x16, 0x14, 0x11, 0x0f, 0x0c, 0x09, 0x07, 0x04,
    0x2a, 0x27, 0x25, 0x22, 0x1f, 0x1d, 0x1a, 0x18, 0x15, 0x13, 0x10, 0x0e, 0x0b, 0x09, 0x06, 0x03,
    0x29, 0x26, 0x24, 0x21, 0x1f, 0x1c, 0x1a, 0x17, 0x15, 0x12, 0x10, 0x0d, 0x0b, 0x08, 0x06, 0x03,
    0x29, 0x26, 0x24, 0x21, 0x1f, 0x1d, 0x1a, 0x18, 0x15, 0x13, 0x10, 0x0e, 0x0b, 0x09, 0x06, 0x04,
    0x2a, 0x27, 0x25, 0x22, 0x20, 0x1d, 0x1b, 0x19, 0x16, 0x14, 0x11, 0x0f, 0x0d, 0x0a, 0x08, 0x05,
    0x2b, 0x29, 0x26, 0x24, 0x21, 0x1f, 0x1d, 0x1a, 0x18, 0x15, 0x13, 0x11, 0x0e, 0x0c, 0x0a, 0x07,
    0x25, 0x23, 0x20, 0x1e, 0x1b, 0x19, 0x17, 0x14, 0x12, 0x10, 0x0d, 0x0b, 0x09, 0x06, 0x04, 0x02,
    0x27, 0x25, 0x23, 0x20, 0x1e, 0x1c, 0x1a, 0x17, 0x15, 0x13, 0x10, 0x0e, 0x0c, 0x09, 0x07, 0x05,
    0x23, 0x20, 0x1e, 0x1c, 0x1a, 0x17, 0x15, 0x13, 0x10, 0x0e, 0x0c, 0x0a, 0x07, 0x05, 0x03, 0x01,
    0x26, 0x24, 0x22, 0x20, 0x1d, 0x1b, 0x19, 0x17, 0x15, 0x12, 0x10, 0x0e, 0x0c, 0x09, 0x07, 0x05,
    0x23, 0x21, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x08, 0x06, 0x04, 0x02,
    0x28, 0x25, 0x23, 0x21, 0x1f, 0x1d, 0x1b, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x09, 0x07,
    0x25, 0x23, 0x21, 0x1f, 0x1d, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x09, 0x07, 0x05,
    0x23, 0x21, 0x1f, 0x1d, 0x1b, 0x19, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06, 0x04,
    0x22, 0x1f, 0x1d, 0x1b, 0x19, 0x17, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x09, 0x07, 0x05, 0x03,
    0x21, 0x1f, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06, 0x04, 0x02,
    0x20, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06, 0x04, 0x02,
    0x20, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06, 0x04, 0x02,
    0x20, 0x1f, 0x1d, 0x1b, 0x19, 0x17, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0b, 0x09, 0x07, 0x05, 0x03,
    0x21, 0x1f, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x07, 0x05,
    0x23, 0x21, 0x1f, 0x1d, 0x1b, 0x19, 0x17, 0x15, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x08, 0x06,
    0x1c, 0x1b, 0x19, 0x17, 0x15, 0x13, 0x11, 0x0f, 0x0e, 0x0c, 0x0a, 0x08, 0x06, 0x04, 0x02, 0x01,
    0x1f, 0x1d, 0x1b, 0x19, 0x17, 0x15, 0x14, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x09, 0x07, 0x05, 0x03,
    0x21, 0x1f, 0x1e, 0x1c, 0x1a, 0x18, 0x16, 0x15, 0x13, 0x11, 0x0f, 0x0d, 0x0c, 0x0a, 0x08, 0x06,
    0x1c, 0x1b, 0x19, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0e, 0x0c, 0x0a, 0x09, 0x07, 0x05, 0x03, 0x02,
    0x20, 0x1e, 0x1c, 0x1a, 0x19, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0e, 0x0c, 0x0b, 0x09, 0x07, 0x05,
    0x1c, 0x1a, 0x18, 0x16, 0x15, 0x13, 0x11, 0x0f, 0x0e, 0x0c, 0x0a, 0x08, 0x07, 0x05, 0x03, 0x01,
    0x20, 0x1e, 0x1c, 0x1b, 0x19, 0x17, 0x15, 0x14, 0x12, 0x10, 0x0f, 0x0d, 0x0b, 0x09, 0x08, 0x06,
    0x1c, 0x1b, 0x19, 0x17, 0x15, 0x14, 0x12, 0x10, 0x0f, 0x0d, 0x0b, 0x0a, 0x08, 0x06, 0x05, 0x03,
    0x19, 0x17, 0x16, 0x14, 0x12, 0x11, 0x0f, 0x0d, 0x0c, 0x0a, 0x08, 0x07, 0x05, 0x03, 0x02, 0x00,
    0x1e, 0x1d, 0x1b, 0x19, 0x18, 0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0c, 0x0b, 0x09, 0x07, 0x06,
    0x1c, 0x1a, 0x19, 0x17, 0x16, 0x14, 0x12, 0x11, 0x0f, 0x0d, 0x0c, 0x0a, 0x08, 0x07, 0x05, 0x04,
    0x1a, 0x18, 0x17, 0x15, 0x14, 0x12, 0x10, 0x0f, 0x0d, 0x0b, 0x0a, 0x08, 0x07, 0x05, 0x03, 0x02,
    0x18, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0f, 0x0d, 0x0c, 0x0a, 0x08, 0x07, 0x05, 0x04, 0x02, 0x00,
    0x1f, 0x1d, 0x1c, 0x1a, 0x19, 0x17, 0x15, 0x14, 0x12, 0x11, 0x0f, 0x0e, 0x0c, 0x0a, 0x09, 0x07,
    0x1e, 0x1c, 0x1b, 0x19, 0x18, 0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0b, 0x0a, 0x08, 0x07,
    0x1d, 0x1c, 0x1a, 0x18, 0x17, 0x15, 0x14, 0x12, 0x11, 0x0f, 0x0e, 0x0c, 0x0b, 0x09, 0x08, 0x06,
    0x1d, 0x1b, 0x1a, 0x18, 0x17, 0x15, 0x13, 0x12, 0x10, 0x0f, 0x0d, 0x0c, 0x0a, 0x09, 0x07, 0x06,
    0x1c, 0x1b, 0x19, 0x18, 0x16, 0x15, 0x13, 0x12, 0x10, 0x0f, 0x0d, 0x0c, 0x0b, 0x09, 0x08, 0x06,
    0x1d, 0x1b, 0x1a, 0x18, 0x17, 0x15, 0x14, 0x12, 0x11, 0x0f, 0x0e, 0x0c, 0x0b, 0x09, 0x08, 0x06,
    0x1d, 0x1c, 0x1a, 0x19, 0x17, 0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0b, 0x0a, 0x09, 0x07,
    0x16, 0x14, 0x13, 0x11, 0x10, 0x0e, 0x0d, 0x0c, 0x0a, 0x09, 0x07, 0x06, 0x04, 0x03, 0x02, 0x00,
    0x17, 0x15, 0x14, 0x12, 0x11, 0x10, 0x0e, 0x0d, 0x0b, 0x0a, 0x08, 0x07, 0x06, 0x04, 0x03, 0x01,
};

#define IEEE_BINARY32_EXPO_BIAS (127)
#define IEEE_BINARY32_MANT_BITS (23)
#define IEEE_BINARY32_EXPO_BITS (8)
#define IEEE_BINARY32_EXPO_MASK (0x7f800000)
#define IEEE_BINARY32_NAN_INDEF (0xffc00000)
#define IEEE_BINARY32_POS_INF   (0x7f800000)
#define IEEE_BINARY32_POS_ZERO  (0x00000000)
#define IEEE_BINARY32_MIN_NORM  (0x00800000)
#define IEEE_BINARY32_SIGN_BIT  (0x80000000)

/* Emulate the RSQRTSS instruction in software */
float my_rsqrtf (float x)
{
    float r;
    uint32_t arg, res, idx, expo, mant;

    arg = float_as_uint32 (x);
    /* zeros and subnormals */
    if ((arg & ~IEEE_BINARY32_SIGN_BIT) < IEEE_BINARY32_MIN_NORM) {
        res = IEEE_BINARY32_POS_INF | (arg & IEEE_BINARY32_SIGN_BIT);
        r = uint32_as_float (res);
    } 
    /* NaNs */
    else if ((arg & ~IEEE_BINARY32_SIGN_BIT) > IEEE_BINARY32_POS_INF) {
        r = x + x; // convert SNaN to QNaN
    }
    /* negative arguments */
    else if (arg & IEEE_BINARY32_SIGN_BIT) {
        res = IEEE_BINARY32_NAN_INDEF;
        r = uint32_as_float (res);
    } 
    /* positive infinity */
    else if (arg == IEEE_BINARY32_POS_INF) {
        res = IEEE_BINARY32_POS_ZERO;
        r = uint32_as_float (res);
    } 
    /* positive normals */
    else {
        /* extract exponent lsb and leading mantissa bits for table index */
        expo = (arg & IEEE_BINARY32_EXPO_MASK) >> IEEE_BINARY32_MANT_BITS;
        idx = (arg >> (IEEE_BINARY32_MANT_BITS - LOG2_NBR_TAB_ENTRIES + 1))
               & (NBR_TAB_ENTRIES - 1);
        /* compute exponent and mantissa of reciprocal square root */
        expo = (3 * IEEE_BINARY32_EXPO_BIAS + ~expo) >> 1;
        mant = (((base_tab [idx >> 4] << (TAB_ENTRY_BITS - BASE_TAB_ENTRY_BITS)) + ofs_tab [idx])
                << (IEEE_BINARY32_MANT_BITS - TAB_ENTRY_BITS));
        /* combine exponent and mantissa bits to compute final result */
        res = (expo << IEEE_BINARY32_MANT_BITS) | mant;
        r = uint32_as_float (res);
    }
    return r;
}

#define NBR_RND_MODES (4)
#define NBR_DAZ_MODES (2)
#define NBR_FTZ_MODES (2)

int main (void)
{
    const uint32_t rnd_mode [NBR_RND_MODES] = 
    {
        _MM_ROUND_NEAREST,
        _MM_ROUND_TOWARD_ZERO, 
        _MM_ROUND_DOWN, 
        _MM_ROUND_UP
    };
    const uint32_t ftz_mode [NBR_FTZ_MODES] = 
    {
        _MM_FLUSH_ZERO_OFF,
        _MM_FLUSH_ZERO_ON
    };
    const uint32_t daz_mode [NBR_DAZ_MODES] = 
    {
        _MM_DENORMALS_ZERO_OFF,
        _MM_DENORMALS_ZERO_ON
    };
    uint32_t iarg, ires, iref;
    float arg, res, ref;
    double relerr, maxrelerr;

    for (int rnd = 0; rnd < NBR_RND_MODES; rnd++) {
        printf ("rnd=%d\n", rnd);
        for (int ftz = 0; ftz < NBR_FTZ_MODES; ftz++) {
            printf (" ftz=%d\n", ftz);
            for (int daz = 0; daz < NBR_DAZ_MODES; daz++) {
                printf ("  daz=%d\n", daz); fflush(stdout);

                maxrelerr = 0;

                iarg = 0;
                do {
                    arg = uint32_as_float (iarg);
                    ref = sse_rsqrtf (arg, daz_mode[daz], ftz_mode[ftz], rnd_mode[rnd]);
                    res = my_rsqrtf (arg);
                    if ((arg >= 1.17549435e-38f) && (arg < 3.40282347e+38f)) { /* normals only */
                        relerr = fabs ((ref - sqrt(1.0/(double)arg)) / sqrt(1.0/(double)arg));
                        if (relerr > maxrelerr) maxrelerr = relerr;
                    }

                    iref = float_as_uint32 (ref);
                    ires = float_as_uint32 (res);
                    if (ires != iref) {
                        printf ("!!!! rnd=%d ftz=%d daz=%d  arg=%08x  res=%08x  ref=%08x\n", 
                                rnd, ftz, daz, iarg, ires, iref);
                        return EXIT_FAILURE;
                    }
                    iarg++;
                } while (iarg);

                printf ("    maxrelerr = %15.8e\n", maxrelerr);
            }
        }
    }
    printf ("RSQRTSS emulation test passed\n");
    return EXIT_SUCCESS;
}
1 голос
/ 30 октября 2019

Для чего бы это ни стоило, я в итоге реализовал его в виде простой сборки в C ++, как предложил @ Francois Andrieux (точнее я использовал ASMJIT ).

Это работает хорошо, хотяимеет недостаток, заключающийся в потере переносимости (хотя и меньше, чем при простом ассемблере). Но это несколько присуще моему вопросу, поскольку я ХОЧУ использовать очень специфическую функцию x86.

Вот мой код: 

   typedef float(*JITFunc)();

   JITFunc func;
   asmjit::JitRuntime jit_runtime;
   asmjit::CodeHolder code;
   code.init(jit_runtime.getCodeInfo());

   asmjit::X86Compiler cc(&code);
   cc.addFunc(asmjit::FuncSignature0<float>());

   float value = 2.71; // Some example value.
   asmjit::X86Xmm x = cc.newXmm();
   setXmmVar(cc, x, value);
   cc.rsqrtss(x, x);   // THE asm function.

   cc.ret(x);

   cc.endFunc();
   cc.finalize();

   jit_runtime.add(&func, &code);

   return func(); // Or something to that effect. func() is the result, anyway.
...