FLT_EPSILON用于具有SSE / AVX的第n个根查找器

我正在尝试转换一个函数,在C中找到第n个根,从以下链接http://rosettacode.org/wiki/Nth_root#C获取一个double值,使用AVX一次找到8个浮点数的第n个根。

部分代码使用DBL_EPSILON * 10.但是,当我将其转换为使用浮点数和AVX时,我必须使用FLT_EPSILON * 1000或代码挂起并且不会收敛。 当我打印出FLT_EPSILON时,我看到它是订单1E-7。 但是这个链接http://www.cplusplus.com/reference/cfloat/表示它应该是1E-5。 当我打印出DBL_EPSILON时它是1E-16但链接说它应该只是1E-9。 这是怎么回事?

这是迄今为止的代码(未完全优化)。

#include  #include  #include  // AVX inline double abs_(double x) { return x >= 0 ? x : -x; } double pow_(double x, int e) { double ret = 1; for (ret = 1; e; x *= x, e >>= 1) { if ((e & 1)) ret *= x; } return ret; } double root(double a, int n) { double d, x = 1; x = a/n; if (!a) return 0; //if (n < 1 || (a = abs_(x) * (DBL_EPSILON * 10)); printf("%d\n", cnt); return x; } __m256 pow_avx(__m256 x, int e) { __m256 ret = _mm256_set1_ps(1.0f); for (; e; x = _mm256_mul_ps(x,x), e >>= 1) { if ((e & 1)) ret = _mm256_mul_ps(x,ret); } return ret; } inline __m256 abs_avx (__m256 x) { return _mm256_max_ps(_mm256_sub_ps(_mm256_setzero_ps(), x), x); //return x >= 0 ? x : -x; } int get_mask(const __m256 d, const __m256 x) { __m256 ad = abs_avx(d); __m256 ax = abs_avx(x); __m256i mask = _mm256_castps_si256(_mm256_cmp_ps(ad, ax, _CMP_GT_OQ)); return _mm_movemask_epi8(_mm256_castsi256_si128(mask)) + _mm_movemask_epi8(_mm256_extractf128_si256(mask,1)); } __m256 root_avx(__m256 a, int n) { printf("%e\n", FLT_EPSILON); printf("%e\n", DBL_EPSILON); printf("%e\n", FLT_EPSILON*1000.0f); __m256 d; __m256 x = _mm256_set1_ps(1.0f); //if (!a) return 0; //if (n < 1 || (a = abs_(x) * (DBL_EPSILON * 10)); } while (get_mask(d, xtmp)); return x; } int main() { __m256 d = _mm256_set1_ps(16.0f); __m256 out = root_avx(d, 4); float result[8]; int i; _mm256_storeu_ps(result, out); for(i=0; i<8; i++) { printf("%f\n", result[i]); } printf("\n"); //double x = 16; //printf("root(%g, 15) = %g\n", x, root(x, 4)); //double x = pow_(-3.14159, 15); //printf("root(%g, 15) = %g\n", x, root(x, 15)); return 0; } 

_mm256_rcp_ps映射到rcpps指令,只执行近似倒数。 英特尔64和IA-32架构软件开发人员手册称其相对误差可能高达1.5•2 -12 。 这不足以使根探测器以精确度100*FLT_EPSILON收敛。

您可以使用精确的除法,例如:

 d = pow_avx(x, n-1); d = _mm256_sub_ps(_mm256_div_ps(a, d), x); 

或者为倒数估计添加一些细化步骤。

顺便提一下,如果您的编译器支持将常规C运算符与SIMD对象一起使用,请考虑使用常规C运算符:

 d = pow_avx(x, n-1); d = a/d - x; 

1e-5只是C标准允许实现用于FLT_EPSILON 。 在实践中,您将使用IEEE-754单精度,其ε为2 -23 ,大约为1e-7