Login/Register
Is there a C++ function that returns exactly the value of the built-in CPU operation RSQRTSS for inverse square root?
Asked Answered
Solved c++x86floating-pointssesqrt
W

3

6

I am looking for a C++ function that returns the inverse sqrt of a float: rsqrt(x) = 1/sqrt(x) by using the exact method like the built-in XMM operation RSQRTSS (cf. https://www.felixcloutier.com/x86/rsqrtss). (I.e., I want the built-in approximation rather than the more precise 1/sqrtf, and I don't care about speed (a lot).)

According to this question:

Is there a fast C or C++ standard library function for double precision inverse square root?

...there is at least no "fast way with double precision" to accomplish this with a standard C++ library. But how to do it slow, non-standard and with floats?

Wretch answered 29/10, 2019 at 19:31 Comment(4)
not fast and not precise, sounds interesting :P, out of curiosity: what do you need it for? – Schwing
From "non-standard" it sounds like you don't care about portability. And you seem to know exactly the assembly instruction you would like executed. So you can just use asm and write a bit of assembly into your code. – Husein
@formerlyknownas_463035818 Speculation: Could be useful for constant propagation in a compiler. – Hollister
@formerlyknownas_463035818 I'm working within a larger project which allows mathematical formulas to be either interpreted for evaluation or pre-compiled (which should both gain the same result). When compiled, the rsqrt function uses the RSQRTSS assembly operation, but the interpreted version doesn't - which, during unit tests, leads to failure messages. I want to fix that mostly to avoid the failure messages (the exact result is not even that important), not really caring about speed on the interpreter side since it's a 1000 times slower than the compiled version anyway. – Wretch
H
6

The RSQRTSS instruction is readily accessible via the _mm_rsqrt_ss() intrinsic declared in immintrin.h. But we can also emulate the instruction in software, as is done in the my_rsqrtf() function below. By simply observing the output of RSQRTSS one easily finds that its function values are based on a (virtual) table of 211 entries, each 12 bit in size.

Note the "virtual" attribute, because it is unlikely that the hardware employs a straight 24 Kbit table. My analysis of the patterns in the table entries doesn't suggest the use of bipartite tables. A much simpler compression scheme -- like the one I used in the code below -- based on a table of base values and a table of offsets may be used. My scheme requires only a narrow adder but reduces ROM storage to 13 Kbit, i.e. almost half.

The implementation below was developed on and tested against an Intel Xeon Processor E3-1270 V2 that uses the Ivy Bridge architecture. There may be some functional differences in the implementation of RSQRTSS between various Intel architectures, and such differences are likely between architectures from different x86-84 vendors.

The framework below checks that the emulation by my_rsqrtf() delivers bit-wise identical results to RSQRTSS for all four rounding modes, two DAZ (denormals are zero) modes, and two FTZ (flush to zero) modes. We find that the function results are not affected by any modes, which matches with how Intel specified RSQRTSS in the Intel® 64 and IA-32 Architectures Software Developer’s Manual:

The RSQRTSS instruction is not affected by the rounding control bits in the MXCSR register. When a source value is a 0.0, an ∞ of the sign of the source value is returned. A denormal source value is treated as a 0.0 (of the same sign). When a source value is a negative value (other than −0.0), a floating-point indefinite is returned. When a source value is an SNaN or QNaN, the SNaN is converted to a QNaN or the source QNaN is returned.

#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;
}
Hollister answered 30/10, 2019 at 21:7 Comment(0)
C
5

There is no function in the standard library that does this, but your compiler might optimize the expression 1 / sqrt(value) such that it does emit the RSQRTSS instruction. For example, with the compiler flags -ffast-math -march=native, GCC will emit that instruction, see: https://godbolt.org/z/cL6seG

Calia answered 29/10, 2019 at 19:43 Comment(1)
I see, this appears to gain the operation I want. I'm afraid I can't change compiler flags, though. Also, why did you use -march=native, it appears to work without as well: godbolt.org/z/nG4mZC – Wretch
W
2

For what it's worth, I ended up implementing it in plain assembly within C++, as @François Andrieux suggested (more precisely I used ASMJIT).

This works well, though it comes with the drawback of losing portability (less than with plain asm, though). But this is somewhat inherent to my question since I WANT to use a very specific x86 function.

Here's my code:

   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.
Wretch answered 30/10, 2019 at 9:20 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.