|
// Public Domain under http://unlicense.org, see link for details. |
|
|
|
// compile with (or equiv): |
|
// clang -O3 -march=native -Wall -Wextra -Wconversion -Wpedantic -Wno-unused-function -fno-math-errno rsqrt.c -o rsqrt -lm |
|
|
|
// VALIDATE_CR options additionally requires MPFR: -lmpfr |
|
|
|
// This is intended to be a hackable testbed. Only a subset of |
|
// included routines are ran as is (those that seem potentially |
|
// interesting to me). Also note that a number of implementations |
|
// differ from those in the original papers. |
|
|
|
|
|
//******************************************************** |
|
// CONFIG STUFF HERE |
|
|
|
// not really exhaustive (add notes). actually I removed this |
|
// because I'm too lazy to check the domains of the inputs. |
|
// (general slightly smaller than all normal inputs). yeah |
|
// this is useful information. |
|
//#define EXHAUSTIVE_TEST |
|
|
|
// following options requires: MPFR installed (https://www.mpfr.org/) |
|
// if defined. |
|
|
|
// validate correct rounding of double promotion method (for normal |
|
// and denormal input) |
|
//#define VALIDATE_CR |
|
|
|
// if defined (with VALIDATE_CR) then tests all (about 10 minutes) |
|
//#define VALIDATE_CR_ALL |
|
|
|
// END OF CONFIG |
|
//******************************************************** |
|
|
|
|
|
#if defined(VALIDATE_CR) |
|
#define NEEDS_MPFR |
|
#endif |
|
|
|
#include <stdint.h> |
|
#include <stdlib.h> |
|
#include <stdio.h> |
|
#include <math.h> |
|
#include <string.h> |
|
#include <assert.h> |
|
|
|
|
|
//******************************************************** |
|
|
|
const float f32_min_normal = 0x1.0p-126f; |
|
const float f32_max_normal = 0x1.fffffep+127f; |
|
|
|
static inline uint32_t f32_to_bits(float x) |
|
{ |
|
uint32_t u; memcpy(&u,&x,4); return u; |
|
} |
|
|
|
inline float f32_from_bits(uint32_t x) |
|
{ |
|
float u; memcpy(&u,&x,4); return u; |
|
} |
|
|
|
static inline uint64_t f64_to_bits(double x) |
|
{ |
|
uint64_t u; memcpy(&u, &x, 8); return u; |
|
} |
|
|
|
static inline double f64_from_bits(uint64_t x) |
|
{ |
|
double f; memcpy(&f, &x, 8); return f; |
|
} |
|
|
|
// JUST SAY NO TO MATH ERRNO |
|
|
|
// to cut some of the pain of math errno not being disabled |
|
// (-fno-math-errno). But you really should do that. Unless |
|
// your using visual-c: Then I feel your pain and I can't |
|
// think of anything better to do than this. (and nag |
|
// MS to make a work around) |
|
inline float f32_sqrt(float a) |
|
{ |
|
#if defined(__NO_MATH_ERRNO__) || defined(_M_FP_FAST) |
|
return sqrtf(a); |
|
#elif defined(F32_INTEL) |
|
return F32_SSE_WRAP(_mm_sqrt_ss, a); |
|
#elif defined(__ARM_ARCH) |
|
return __sqrtf(a); |
|
#else |
|
// do nothing for warning ATM |
|
#endif |
|
} |
|
|
|
//******************************************************** |
|
|
|
// faithfully rounded. |
|
static inline float f32_rsqrt(float x) { return f32_sqrt(1.f/x); } |
|
|
|
// promote/demote double is correctly rounded |
|
static inline float f32_rsqrt_cr(float x) |
|
{ |
|
return (float)sqrt(1.0/(double)x); |
|
} |
|
|
|
// DEAD |
|
// dep-chain len: 2 |
|
// correctly: 74.013561% (12417415) |
|
// faithfully: 25.986439% (4359801) |
|
// uiCA predicted throughput: 23.00 |
|
float rsqrt_n1(float x) { return f32_sqrt(x)/x; } |
|
|
|
// DEAD |
|
// dep-chain len: 2 |
|
// correctly: 73.995733% (12414424) |
|
// faithfully: 26.004267% (4362792) |
|
// uiCA predicted throughput: 23.00 |
|
float rsqrt_n2(float x) { return 1.f/f32_sqrt(x); } |
|
|
|
// MOSTLY DEAD |
|
// NOTE: if sqrt(x) and 1/x "can" be computed in parallel |
|
// then this becomes potentially interesting. This is |
|
// not the case in current mainstream CPUs. |
|
|
|
// dep-chain len: 2 |
|
// correctly: 64.843583% (10878948) |
|
// faithfully: 34.837395% (5844745) |
|
// 2 ulp: 0.319022% (53523) |
|
float rsqrt_n3(float x) { return f32_sqrt(x)*(1.f/x); } |
|
|
|
//******************************************************** |
|
|
|
// Newton-Raphson step for 1/sqrt(x) using FMA when |
|
// we haven't computed (or approximated) 1/x |
|
// y = estimated in (y_n) |
|
// h = -x/2 |
|
// { don't use see: XX} |
|
inline float f32_rsqrt_nr_step_common(float y, float h) |
|
{ |
|
return y*fmaf(h,y*y, 1.5f); |
|
} |
|
|
|
// rename |
|
static inline float f32_rsqrt_nr_fma_step(float y, float h) |
|
{ |
|
float r = y*y; |
|
float v = fmaf(h,r, 0.5f); |
|
return fmaf(y,v,y); |
|
} |
|
|
|
// Newton-Raphson step for 1/sqrt(x) using FMA when |
|
// y = estimate in (y_n) |
|
// h = -x/2 |
|
// r = ~1/x |
|
// RENAME AND NO |
|
inline float f32_rsqrt_nr_step_cb(float y, float h, float r) |
|
{ |
|
float s = fmaf(h,r,0.5f); |
|
float t = fmaf(y,y, -r); |
|
float v = fmaf(h,t, s); |
|
return fmaf(y,v,y); |
|
} |
|
|
|
// Halley method step for 1/sqrt(x) |
|
// y = estimate in (y_n) |
|
// h = -x/2 |
|
// r = ~1/x |
|
inline float f32_rsqrt_hm_step(float y, float h, float r) |
|
{ |
|
float s = fmaf(r,h,0.5f); |
|
float t = fmaf(y,y,-r); |
|
float v = fmaf(h,t, s); |
|
float w = fmaf(1.5f*v,v,v); |
|
return fmaf(y,w,y); |
|
} |
|
|
|
//******************************************************** |
|
// expansion helpers |
|
|
|
static inline float rsqrt_r1_apply(float (*f)(float), float x) |
|
{ |
|
float y = f(x); |
|
return f32_rsqrt_nr_fma_step(y, -0.5f*x); |
|
} |
|
|
|
static inline float rsqrt_r2_apply(float (*f)(float), float x) |
|
{ |
|
float y = f(x); |
|
return f32_rsqrt_nr_step_cb(y, -0.5f*x, y*y); |
|
} |
|
|
|
static inline float rsqrt_hm_apply(float (*f)(float), float x) |
|
{ |
|
float y = f(x); |
|
return f32_rsqrt_hm_step(y, -0.5f*x, y*y); |
|
} |
|
|
|
//******************************************************** |
|
|
|
// Portable versions by Carlos F. Borges |
|
|
|
// "Fast Compensated Algorithms for the Reciprocal Square Root, the Reciprocal Hypotenuse, and Givens Rotations", Carlos F. Borges, 2021 |
|
// https://arxiv.org/abs/2103.08694 |
|
// -- |
|
// additional analysis: |
|
// "High-level algorithms for correctly-rounded reciprocal square roots", Borges, Jeannerod & Muller, 2022 |
|
// https://hal.inria.fr/hal-03728088 |
|
|
|
// dep-chain len: 5 |
|
// correctly: 99.999994% (16777215) |
|
// faithfully: 0.000006% (1) |
|
// Borges Algorithm 2 (rsqrtNewton in Borges/Jeannerod/Muller) |
|
// uiCA: https://bit.ly/3OkNAk5 |
|
// predicted throughput: 35.00 |
|
float rsqrt_hq(float x) |
|
{ |
|
float r = 1.f/x; |
|
float h = -0.5f*x; |
|
|
|
x = f32_sqrt(r); |
|
x = f32_rsqrt_nr_step_cb(x,h,r); |
|
|
|
return x; |
|
} |
|
|
|
// correctly rounded (in round to even/ties to nearest) |
|
// on [0x1p-125, 0x1.fffffcp+127] (almost all normal numbers) |
|
// Borges Algorithm 3 (rsqrtHalley in Borges/Jeannerod/Muller) |
|
// uiCA: https://bit.ly/3AzVQXJ |
|
// predicted throughput: 40.00 |
|
float rsqrt_cb(float x) |
|
{ |
|
float r = 1.f/x; |
|
float h = -0.5f*x; |
|
|
|
x = f32_sqrt(r); |
|
x = f32_rsqrt_hm_step(x,h,r); |
|
|
|
return x; |
|
} |
|
|
|
|
|
//******************************************************** |
|
// so-call "bithack" initial steps |
|
|
|
static inline float f32_rsqrt_initial(float x, uint32_t K) |
|
{ |
|
return f32_from_bits(K - (f32_to_bits(x) >> 1)); |
|
} |
|
|
|
// Day observes there's a second choice which gives some |
|
// extra wiggle room in optimizing the "magic" K |
|
static inline float f32_rsqrt_initial_md(float x, uint32_t K) |
|
{ |
|
return f32_from_bits((K-f32_to_bits(x)) >> 1); |
|
} |
|
|
|
//******************************************************** |
|
|
|
// Mike Day core routines as in paper. Without star is paper |
|
// and contractions are disabled. Starred versions are explict FMAs |
|
// where version choosen uses the greedy measure to minimize |
|
// the counts of >2 ulp results. |
|
// | func |max ULP| CR| FR| 2 ULP| > 2 ULP| |
|
// | ---: | ---:| ---:| ---:| ---:| ---:| |
|
// | FRSR_Mon0 | 564177| 12| 28| 23| 25165761| |
|
// | FRSR_Deg0 | 403258| 11| 26| 24| 25165763| |
|
// | FRSR_Mon1 | 14751| 557| 1070| 1063| 25163134| |
|
// | FRSR_Mon1*| 14751| 545| 1057| 1099| 25163123| |
|
// | FRSR_Deg1 | 10215| 941| 1802| 1939| 25161142| |
|
// | FRSR_Deg1*| 10215| 1051| 1848| 1963| 25160962| |
|
// | FRSR_Deg1Alt | 10219| 874| 1858| 2033| 25161059| |
|
// | FRSR_Deg1Alt*| 10219| 972| 1920| 1901| 25161031| |
|
// | FRSR_Mon2 | 317|105658| 214279| 222937| 24622950| |
|
// | FRSR_Mon2*| 317|105658| 214279| 222937| 24622950| |
|
|
|
// listing 2 (FRSR_Mon0) |
|
// max-ulp: 564177 |
|
float FRSR_Mon0(float x) |
|
{ |
|
return f32_rsqrt_initial(x, 0x5F37642F); |
|
} |
|
|
|
// listing 3 (FRSR_Deg0) |
|
// max-ulp: 403258 |
|
float FRSR_Deg0(float x) |
|
{ |
|
return 0.79247999f * f32_rsqrt_initial_md(x, 0xBEBFFDAA); |
|
} |
|
|
|
// listing 4 (FRSR_Mon1) |
|
// max-ulp: 14751 |
|
// | CR| FR | 2 ULP| > 2 ULP | |
|
// | 557| 1070| 1063| 25163134| 1 |
|
// | 497| 1094| 1097| 25163136| 2 |
|
// | 545| 1057| 1099| 25163123| 3 |
|
// 3 fma/mul |
|
float FRSR_Mon1(float x) |
|
{ |
|
float y = f32_rsqrt_initial_md(x, 0xBE167122); |
|
|
|
//return y * (1.8909901f - x*y*y); // 1 |
|
//return y * fmaf(-x, y*y, 1.8909901f); // 2 |
|
return y * fmaf(-y, x*y, 1.8909901f); // 3 |
|
} |
|
|
|
// listing 5 (FRSR_Deg1) |
|
// max-ulp: 10215 |
|
// | CR| FR | 2 ULP| > 2 ULP | |
|
// | 941| 1802| 1939| 25161142| 1 |
|
// | 992| 1907| 1882| 25161043| 2 |
|
// | 1051| 1848| 1963| 25160962| 3 |
|
// | 976| 1907| 1894| 25161047| 4 |
|
// 4 fma/mul |
|
float FRSR_Deg1(float x) |
|
{ |
|
static const float c0 = 1.1893165f; |
|
static const float c1 = -0.24889956f; |
|
|
|
float y = f32_rsqrt_initial(x, 0x5F5FFF00); |
|
|
|
//return y * (c0 + x*y*y*c1); // 1 |
|
//return y * fmaf(x, y*y*c1, c0); // 2 |
|
return y * fmaf(x*y, y*c1, c0); // 3 |
|
//return y * fmaf(x*c1,y*y, c0); // 4 |
|
} |
|
|
|
// listing 6 (FRSR_Deg1Alt) |
|
// max-ulp: 10219 |
|
// | CR| FR | 2 ULP| > 2 ULP | |
|
// | 874| 1858| 2033| 25161059| 1 |
|
// | 972| 1920| 1901| 25161031| 2 |
|
// | 950| 1991| 1908| 25160975| 3 |
|
// | 972| 1920| 1901| 25161031| 4 |
|
// 4 fma/mul |
|
float FRSR_Deg1Alt(float x) |
|
{ |
|
static const float c0 = 1.1891762f; |
|
static const float c1 = -0.24881148f; |
|
|
|
float y = f32_rsqrt_initial(x, 0x5F6004CC); |
|
|
|
//return y * (c0 + y*y*x*c1); // 1 |
|
//return y * fmaf(y*y, x*c1, c0); // 2 |
|
//return y * fmaf(x*y, y*c1, c0); // 3 |
|
return y * fmaf(x*c1,y*y, c0); // 4 |
|
} |
|
|
|
// listing 7 (FRSR_Mon2) |
|
// max-ulp: 317 |
|
// | CR| FR | 2 ULP| > 2 ULP | |
|
// | 105658| 214279| 222937| 24622950| 1 |
|
// | 105658| 214279| 222937| 24622950| 2 |
|
// | 106447| 213978| 221946| 24623453| 3 |
|
// 4 fma/mul, 1 add |
|
float rsqrt_day(float x) |
|
{ |
|
static const float c0 = 2.2825186f; |
|
static const float c1 = -2.253305f; |
|
|
|
float y = f32_rsqrt_initial(x, 0x5F11107D); |
|
float z = x*y*y; |
|
|
|
//return y * (c0 + z *(z+c1)); // 1 |
|
return y * fmaf(z, z+c1, c0); // 2 |
|
//return y * fmaf(z,z,fmaf(z,c1,c0)); // 3 |
|
} |
|
|
|
// max-ulp: 2 |
|
// correctly: 70.604507% (17768206) |
|
// faithfully: 29.345020% (7384916) |
|
// 2 ulp: 0.050473% (12702) |
|
// 8 fma/mul, 1 add |
|
float rsqrt_day_r1(float x) { return rsqrt_r1_apply(&rsqrt_day, x); } |
|
|
|
// max-ulp: 1 |
|
// correctly: 99.667398% (25082122) |
|
// faithfully: 0.332602% (83702) |
|
// 10 fma/mul, 1 add |
|
float rsqrt_day_r2(float x) { return rsqrt_r2_apply(&rsqrt_day, x); } |
|
|
|
// max-ulp: 1 |
|
// correctly: 99.999559% (25165713) |
|
// faithfully: 0.000441% (111) |
|
// 12 fma/mul, 1 add |
|
float rsqrt_day_hm(float x) { return rsqrt_hm_apply(&rsqrt_day, x); } |
|
|
|
float rsqrt_day_cr(float x) |
|
{ |
|
float y = rsqrt_day(x); |
|
float h = -0.5f*x; |
|
y = f32_rsqrt_hm_step(y, h, y*y); |
|
//y = f32_rsqrt_nr_fma_step(y, -0.5f*x); // no good with this step |
|
y = f32_rsqrt_nr_step_cb(y, h, y*y); // need this |
|
return y; |
|
} |
|
|
|
//********************************************************************* |
|
|
|
// Walczyk, Moroz, Cieslinski, Cezary: |
|
// "InvSqrt2" - modified to use FMA |
|
// no constant refinement |
|
// max-ulp: 11 |
|
// correctly: 5.120015% (858996) |
|
// faithfully: 11.176844% (1814391) |
|
// 2 ulp: 12.521724% (3151195) |
|
// >2 ulp: 71.061440% (17883197) |
|
// 7 fma/mul |
|
float rsqrt_wmc(float x) |
|
{ |
|
float h = -0.5f*x; |
|
|
|
x = f32_rsqrt_initial(x, 0x5F376908); |
|
x = x*fmaf(h,x*x, 1.5008789f); |
|
x = x*fmaf(h,x*x, 1.5000006f); |
|
return x; |
|
} |
|
|
|
float rsqrt_wmc_r1(float x) { return rsqrt_r1_apply(&rsqrt_wmc, x); } |
|
float rsqrt_wmc_r2(float x) { return rsqrt_r2_apply(&rsqrt_wmc, x); } |
|
float rsqrt_wmc_hm(float x) { return rsqrt_hm_apply(&rsqrt_wmc, x); } |
|
|
|
//********************************************************************* |
|
// "Modified Fast Inverse Square Root and Square Root Approximation |
|
// Algorithms: The Method of Switching Magic Constants" |
|
// Moroz, Volodymyr, Samotyy and Horyachyy, 2021 |
|
|
|
float rsqrt_msh(float x) |
|
{ |
|
static const uint32_t R[2] = {0x5f19e8b7, 0x5ed9e8b7}; |
|
static const float A[2] = {2.1499474f, 1.0749737f}; |
|
static const float M[2] = {0.824218631f, 2.33124256f}; |
|
|
|
uint32_t id = (f32_to_bits(x) & 0x00800000) == 0 ? 0 : 1; |
|
|
|
float y; |
|
|
|
y = f32_rsqrt_initial(x, R[id]); |
|
y = M[id]*y*fmaf(-x, y*y, A[id]); |
|
|
|
return y; |
|
} |
|
|
|
float rsqrt_msh_r1(float x) { return rsqrt_r1_apply(&rsqrt_msh, x); } |
|
float rsqrt_msh_r2(float x) { return rsqrt_r2_apply(&rsqrt_msh, x); } |
|
float rsqrt_msh_hm(float x) { return rsqrt_hm_apply(&rsqrt_msh, x); } |
|
|
|
//********************************************************************* |
|
// original SSE approximation |
|
|
|
#include <x86intrin.h> |
|
|
|
// return hardware 1/sqrt(x) approx |
|
static inline float rsqrt_sse(float x) |
|
{ |
|
return _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(x))); |
|
} |
|
|
|
// 4 fma/mul |
|
float rsqrt_sse_r1(float x) |
|
{ |
|
float y = rsqrt_sse(x); |
|
return f32_rsqrt_nr_fma_step(y, -0.5f*x); |
|
} |
|
|
|
// 8 fma/mul |
|
float rsqrt_sse_r2(float x) |
|
{ |
|
float y = rsqrt_sse(x); |
|
return f32_rsqrt_nr_step_cb(y, -0.5f*x, y*y); |
|
} |
|
|
|
// 8 fma/mul |
|
float rsqrt_sse_hm(float x) |
|
{ |
|
float y = rsqrt_sse(x); |
|
return f32_rsqrt_hm_step(y, -0.5f*x, y*y); |
|
} |
|
|
|
// 13 fma/mul |
|
float rsqrt_sse_cr(float x) |
|
{ |
|
float y = rsqrt_sse(x); |
|
float h = -0.5f*x; |
|
y = f32_rsqrt_hm_step(y, h, y*y); |
|
y = f32_rsqrt_nr_step_cb(y, h, y*y); |
|
return y; |
|
} |
|
|
|
//********************************************************************* |
|
// Use hardware approximation: neon |
|
|
|
// emulate neon approximation ops (std & not increased precision ones) |
|
|
|
// ignore all special cases (including 0): don't |
|
// need them for local testing. written to be eyeball-able to spec |
|
float vrecpe_f32(float x) |
|
{ |
|
// extract the parts |
|
uint32_t ix = f32_to_bits(x); |
|
uint32_t sx = ((ix >> 31) & 0x01); |
|
uint32_t ex = ((ix >> 23) & 0xff); |
|
uint32_t mx = ((ix ) & 0x7fffff); |
|
uint32_t t; |
|
|
|
mx = 0x100|(mx >> 15); // add implied bit and top 8 explict |
|
ex = 253-ex; // adjust the exp |
|
|
|
// SEE: ARM Spec - FPRecipEstimate |
|
t = (mx<<1)+1; |
|
t = (((1<<19)/t)+1)>>1; |
|
mx = (t & 0xff) << 15; |
|
|
|
// concat the parts back together |
|
ix = (sx & 0x01) << 31; |
|
ix |= (ex & 0xff) << 23; |
|
ix |= (mx & 0x7fffff); |
|
|
|
return f32_from_bits(ix); |
|
} |
|
|
|
// max-ulp: 4947968 |
|
// correctly: 0.001293% (217) |
|
// faithfully: 0.002837% (476) |
|
// written to be eyeball-able to spec |
|
float vrsqrte_f32(float x) |
|
{ |
|
// yeah: lots of reductions are possible & don't care |
|
uint32_t ix = f32_to_bits(x); |
|
uint32_t ex = ((ix >> 23) & 0xff); |
|
uint32_t mx = ((ix ) & 0x7fffff); |
|
uint32_t s = (mx & 1); |
|
|
|
mx >>= (15+s); |
|
mx |= 1 << (8-s); |
|
|
|
if (mx < 256) |
|
mx = mx * 2 + 1; |
|
else |
|
mx = (mx|1) << 1; |
|
|
|
uint32_t t = 512; |
|
|
|
while (mx * (t+1)*(t+1) < (1<<28)) { t++; } |
|
|
|
mx = (t+1)/2; |
|
|
|
ex = (380-ex)/2; |
|
ix = (ex & 0xff) << 23; |
|
ix |= (mx & 0xff) << 15; |
|
|
|
return f32_from_bits(ix); |
|
} |
|
|
|
float vrsqrts_f32(float a, float b) |
|
{ |
|
return -0.5f*fmaf(a,b,-3.f); |
|
} |
|
|
|
float vrecps_f32(float a, float b) |
|
{ |
|
return -fmaf(a,b,-2.f); |
|
} |
|
|
|
|
|
//******************************************************** |
|
// start of testing related stuff |
|
//******************************************************** |
|
|
|
static inline uint32_t u32_abs(uint32_t x) |
|
{ |
|
return (int32_t)x >= 0 ? x : -x; |
|
} |
|
|
|
// ulp distance provided a & b are finite |
|
// and have the same sign |
|
static inline uint32_t f32_ulp_dist_ss(float a, float b) |
|
{ |
|
uint32_t ua = f32_to_bits(a); |
|
uint32_t ub = f32_to_bits(b); |
|
return u32_abs(ua-ub); |
|
} |
|
|
|
static inline uint64_t u64_abs(uint64_t x) |
|
{ |
|
return (int64_t)x >= 0 ? x : -x; |
|
} |
|
|
|
// ulp distance provided a & b are finite |
|
// and have the same sign |
|
static inline uint64_t f64_ulp_dist_ss(double a, double b) |
|
{ |
|
uint64_t ua = f64_to_bits(a); |
|
uint64_t ub = f64_to_bits(b); |
|
return u64_abs(ua-ub); |
|
} |
|
|
|
|
|
typedef struct { |
|
float (*f)(float); |
|
char* name; |
|
char* notes; |
|
uint32_t flags; |
|
uint32_t type; |
|
} func_entry_t; |
|
|
|
#define LENGTHOF(X) (sizeof(X)/sizeof(X[0])) |
|
#define STRINGIFY(S) STRINGIFY_(S) |
|
#define STRINGIFY_(S) #S |
|
#define ENTRY(X) { .f=&X, .name=STRINGIFY(X) } |
|
|
|
func_entry_t func_table[] = |
|
{ |
|
#if 1 |
|
ENTRY(vrsqrte_f32), |
|
ENTRY(FRSR_Mon0), |
|
ENTRY(FRSR_Deg0), |
|
ENTRY(FRSR_Mon1), |
|
ENTRY(FRSR_Deg1), |
|
ENTRY(FRSR_Deg1Alt), |
|
ENTRY(rsqrt_sse), |
|
ENTRY(rsqrt_msh), |
|
ENTRY(rsqrt_day), |
|
ENTRY(rsqrt_wmc), |
|
ENTRY(rsqrt_sse_r2), |
|
ENTRY(rsqrt_sse_r1), |
|
ENTRY(rsqrt_day_r1), |
|
ENTRY(f32_rsqrt), |
|
ENTRY(rsqrt_day_r2), |
|
ENTRY(rsqrt_wmc_r1), |
|
ENTRY(rsqrt_sse_hm), |
|
ENTRY(rsqrt_day_hm), |
|
ENTRY(rsqrt_wmc_r2), |
|
ENTRY(rsqrt_wmc_hm), |
|
//ENTRY(f32_rsqrt_bar), |
|
// find ranges of these |
|
ENTRY(rsqrt_cb), |
|
ENTRY(rsqrt_day_cr), |
|
ENTRY(rsqrt_sse_cr), |
|
#else |
|
// hack for testing one |
|
#endif |
|
}; |
|
|
|
|
|
typedef struct { |
|
double abs; |
|
uint32_t max; |
|
uint32_t ulp[4]; |
|
} func_error_t; |
|
|
|
// tracks total error data across multiple intervals |
|
func_error_t func_error[LENGTHOF(func_table)] = {{0}}; |
|
|
|
|
|
// add error data from current interval to the totals |
|
void error_to_totals(func_error_t* e) |
|
{ |
|
for(uint32_t fi=0; fi < LENGTHOF(func_table); fi++, e++) { |
|
func_error[fi].ulp[0] += e->ulp[0]; |
|
func_error[fi].ulp[1] += e->ulp[1]; |
|
func_error[fi].ulp[2] += e->ulp[2]; |
|
func_error[fi].ulp[3] += e->ulp[3]; |
|
|
|
if (e->max > func_error[fi].max) func_error[fi].max = e->max; |
|
if (e->abs > func_error[fi].abs) func_error[fi].abs = e->abs; |
|
} |
|
} |
|
|
|
void error_dump_i(func_error_t* e) |
|
{ |
|
printf("|%15s|%12s|%10s|%10s|%10s|%10s|%10s|%10s|%10s|%10s|%10s|\n", |
|
"func", "e", "max ULP", "CR", "FR", "2 ULP", "> 2 ULP", |
|
"CR%", "FR%", "2 ULP%","> 2 ULP%"); |
|
|
|
printf("|%15s|%12s|%10s|%10s|%10s|%10s|%10s|%10s|%10s|%10s|%10s|\n", |
|
"---:", ":---:", "---:", "---:", "---:", "---:", |
|
"---:", "---:", "---:","---:", "---:"); |
|
|
|
for(uint32_t fi=0; fi < LENGTHOF(func_table); fi++, e++) { |
|
uint32_t u0 = e->ulp[0]; |
|
uint32_t u1 = e->ulp[1]; |
|
uint32_t u2 = e->ulp[2]; |
|
uint32_t u3 = e->ulp[3]; |
|
uint32_t t = (u0+u1+u2+u3); |
|
double s = 100.0/(double)t; |
|
|
|
if (e->abs < 9.0) |
|
printf("|%15s|%12.10f|%10u|%10u|%10u|%10u|%10u|%10f|%10f|%10f|%10f|\n", |
|
func_table[fi].name, e->abs, e->max, |
|
u0, u1, u2, u3, |
|
s*u0, s*u1, s*u2, s*u3 |
|
); |
|
else |
|
printf("|%15s| -- |%10u|%10u|%10u|%10u|%10u|%10f|%10f|%10f|%10f|\n", |
|
func_table[fi].name, e->max, |
|
u0, u1, u2, u3, |
|
s*u0, s*u1, s*u2, s*u3 |
|
); |
|
|
|
} |
|
} |
|
|
|
void error_dump(void) |
|
{ |
|
printf("\nTOTAL:\n"); |
|
error_dump_i(func_error); |
|
} |
|
|
|
void test_range(uint32_t x0, uint32_t x1) |
|
{ |
|
float f0 = f32_from_bits(x0); |
|
float f1 = f32_from_bits(x1); |
|
|
|
printf("\nchecking on [%08x,%08x] [%e,%e]\n", x0,x1,f0,f1); |
|
|
|
func_error_t error[LENGTHOF(func_table)] = {{0}}; |
|
|
|
for(uint32_t ix=x0; ix<=x1; ix++) { |
|
float x = f32_from_bits(ix); |
|
double d = sqrt(1.0/(double)x); |
|
float cr = (float)d; |
|
|
|
for(uint32_t fi=0; fi < LENGTHOF(func_table); fi++) { |
|
float r = func_table[fi].f(x); |
|
double e = fabs(d-(double)r) * 0x1.0p24; |
|
|
|
//if (e == 0.f) { error[fi].ulp[0]++; continue; }; |
|
|
|
uint32_t ulp = f32_ulp_dist_ss(cr,r); |
|
|
|
if (ulp > error[fi].max) { error[fi].max = ulp; } |
|
if (e > error[fi].abs) { error[fi].abs = e; } |
|
|
|
if (ulp > 3) ulp = 3; |
|
|
|
error[fi].ulp[ulp]++; |
|
} |
|
} |
|
|
|
error_to_totals(error); |
|
error_dump_i(error); |
|
} |
|
|
|
//******************************************************** |
|
// validating double promotion method is correctly rounded |
|
|
|
#if defined(NEEDS_MPFR) |
|
#include <mpfr.h> |
|
|
|
#if defined(VALIDATE_CR) |
|
// worker for vaidate_double_promote |
|
void validate_cr_i(float s, float e) |
|
{ |
|
mpfr_t mx,mr; |
|
|
|
uint32_t xs = f32_to_bits(s); |
|
uint32_t xe = f32_to_bits(e); |
|
|
|
mpfr_init2(mx, 24); |
|
mpfr_init2(mr, 24); |
|
|
|
for(uint32_t xi=xs; xi<=xe; xi++) { |
|
float x = f32_from_bits(xi); |
|
|
|
// compute correctly rounded |
|
mpfr_set_flt(mx, x, MPFR_RNDN); |
|
mpfr_rec_sqrt(mr, mx, MPFR_RNDN); |
|
|
|
float cr = mpfr_get_flt(mr, MPFR_RNDN); |
|
float r = f32_rsqrt_cr(x); |
|
|
|
if (r == cr) continue; |
|
|
|
printf(" ERROR: not correctly rounded at: x=%a %a got %a\n", x, cr, r); |
|
|
|
xe = xi-1; |
|
break; |
|
} |
|
|
|
mpfr_clear(mr); |
|
mpfr_clear(mx); |
|
} |
|
|
|
// double promotion is correctly rounded for all normal |
|
// inputs. |
|
void vaidate_double_promote(void) |
|
{ |
|
|
|
#if defined(VALIDATE_CR_ALL) |
|
printf("starting: validate f32_sqrt_cr: should take around 10 mins and I'm not going to give you any visual feedback of progress.\n"); |
|
validate_cr_i(0, f32_max_normal); |
|
#else |
|
printf("starting: validate f32_sqrt_cr: (fast/non-exhastive)\n"); |
|
printf(" denormals\n"); |
|
validate_cr_i(0, f32_min_normal); |
|
printf(" bottom of normal range\n"); |
|
validate_cr_i(f32_min_normal, f32_min_normal*4.f); |
|
printf(" top of normal range\n"); |
|
validate_cr_i(f32_max_normal/4.f, f32_max_normal); |
|
printf(" around 1.0\n"); |
|
validate_cr_i(1.f/4.f, 4.f); |
|
#endif |
|
} |
|
#endif |
|
|
|
//******************************************************** |
|
|
|
|
|
#endif |
|
|
|
|
|
|
|
int main(void) |
|
{ |
|
#if defined(VALIDATE_CR) |
|
vaidate_double_promote(); |
|
#endif |
|
|
|
printf("<details markdown=\"1\"><summary>click for range breakdown</summary>\n\n"); |
|
|
|
#if defined(EXHAUSTIVE_TEST) |
|
printf("fill-in again\n"); |
|
#else |
|
test_range(f32_to_bits(1.f),f32_to_bits(2.f)); |
|
test_range(f32_to_bits(2.f),f32_to_bits(4.f)); |
|
#endif |
|
|
|
printf("\n</details>\n\n"); |
|
error_dump(); |
|
|
|
return 0; |
|
} |
