c/c++语言开发共享Payne Hanek算法在C中的实现

我正在努力理解如何实现Payne和Hanek发布的范围缩减算法(三角函数的范围缩减)

我见过这个库: http : //www.netlib.org/fdlibm/

但它看起来如此扭曲,我所创立的所有理论解释都太简单了,无法提供实现。

有一些好的……好的……好的解释吗?

    通过Payne-Hanek算法对三角函数进行参数约简实际上非常简单。 与其他参数约简方案一样,计算n = round_nearest (x / (π/2)) ,然后通过x - n * π/2计算余数。 通过计算n = round_nearest (x * (2/π))来实现更高的效率。

    在Payne-Hanek中的关键观察是,当使用完全未连接的乘积计算x - n * π/2的余数时,前导位在减法期间取消,因此我们不需要计算这些。 我们留下的问题是根据x的大小找到正确的起点(非零位)。 如果x接近π/2的倍数,则可能存在额外的消除,这是有限的。 人们可以查阅文献,了解在这种情况下取​​消的附加位数的上限。 由于相对较高的计算成本,Payne-Hanek通常仅用于幅度较大的参数,这具有额外的好处,即在减法期间,原始参数x的比特在相关比特位置中为零。

    下面我展示了我最近编写的单精度sinf()详尽测试C99代码,其中包含减少的慢速路径中的Payne-Hanek减少,请参阅trig_red_slowpath_f() 。 注意,为了实现忠实舍入的sinf() ,必须增加参数减少以将减少的参数作为头/尾方式的两个float操作数返回。

    可以进行各种设计选择,下面我选择了基于整数的计算,以便最小化2/π所需位的存储。 使用浮点计算和重叠对或浮点数的三元组来存储2/π的位的实现也是常见的。

     /* 190 bits of 2/pi for Payne-Hanek style argument reduction. */ static const unsigned int two_over_pi_f [] = { 0x00000000, 0x28be60db, 0x9391054a, 0x7f09d5f4, 0x7d4d3770, 0x36d8a566, 0x4f10e410 }; float trig_red_slowpath_f (float a, int *quadrant) { unsigned long long int p; unsigned int ia, hi, mid, lo, tmp, i; int e, q; float r; ia = (unsigned int)(fabsf (frexpf (a, &e)) * 0x1.0p32f); /* extract 96 relevant bits of 2/pi based on magnitude of argument */ i = (unsigned int)e >> 5; e = (unsigned int)e & 31; hi = two_over_pi_f [i+0]; mid = two_over_pi_f [i+1]; lo = two_over_pi_f [i+2]; tmp = two_over_pi_f [i+3]; if (e) { hi = (hi << e) | (mid >> (32 - e)); mid = (mid << e) | (lo >> (32 - e)); lo = (lo << e) | (tmp >> (32 - e)); } /* compute product x * 2/pi in 2.62 fixed-point format */ p = (unsigned long long int)ia * lo; p = (unsigned long long int)ia * mid + (p >> 32); p = ((unsigned long long int)(ia * hi) << 32) + p; /* round quotient to nearest */ q = (int)(p >> 62); // integral portion = quadrant<1:0> p = p & 0x3fffffffffffffffULL; // fraction if (p & 0x2000000000000000ULL) { // fraction >= 0.5 p = p - 0x4000000000000000ULL; // fraction - 1.0 q = q + 1; } /* compute remainder of x / (pi/2) */ double d; d = (double)(long long int)p; d = d * 0x1.921fb54442d18p-62; // 1.5707963267948966 * 0x1.0p-62 r = (float)d; if (a < 0.0f) { r = -r; q = -q; } *quadrant = q; return r; } /* Like rintf(), but -0.0f -> +0.0f, and |a| must be <= 0x1.0p+22 */ float quick_and_dirty_rintf (float a) { float cvt_magic = 0x1.800000p+23f; return (a + cvt_magic) - cvt_magic; } /* Argument reduction for trigonometric functions that reduces the argument to the interval [-PI/4, +PI/4] and also returns the quadrant. It returns -0.0f for an input of -0.0f */ float trig_red_f (float a, float switch_over, int *q) { float j, r; if (fabsf (a) > switch_over) { /* Payne-Hanek style reduction. M. Payne and R. Hanek, Radian reduction for trigonometric functions. SIGNUM Newsletter, 18:19-24, 1983 */ r = trig_red_slowpath_f (a, q); } else { /* FMA-enhanced Cody-Waite style reduction. WJ Cody and W. Waite, "Software Manual for the Elementary Functions", Prentice-Hall 1980 */ j = 0x1.45f306p-1f * a; // 2/pi j = quick_and_dirty_rintf (j); r = fmaf (j, -0x1.921fb0p+00f, a); // pio2_high r = fmaf (j, -0x1.5110b4p-22f, r); // pio2_mid r = fmaf (j, -0x1.846988p-48f, r); // pio2_low *q = (int)j; } return r; } /* Approximate sine on [-PI/4,+PI/4]. Maximum ulp error = 0.64721 Returns -0.0f for an argument of -0.0f Polynomial approximation based on unpublished work by T. Myklebust */ float sinf_poly (float a, float s) { float r; r = 0x1.7d3bbcp-19f; r = fmaf (r, s, -0x1.a06bbap-13f); r = fmaf (r, s, 0x1.11119ap-07f); r = fmaf (r, s, -0x1.555556p-03f); r = r * s + 0.0f; // ensure -0 is passed trough r = fmaf (r, a, a); return r; } /* Approximate cosine on [-PI/4,+PI/4]. Maximum ulp error = 0.87531 */ float cosf_poly (float s) { float r; r = 0x1.98e616p-16f; r = fmaf (r, s, -0x1.6c06dcp-10f); r = fmaf (r, s, 0x1.55553cp-05f); r = fmaf (r, s, -0x1.000000p-01f); r = fmaf (r, s, 0x1.000000p+00f); return r; } /* Map sine or cosine value based on quadrant */ float sinf_cosf_core (float a, int i) { float r, s; s = a * a; r = (i & 1) ? cosf_poly (s) : sinf_poly (a, s); if (i & 2) { r = 0.0f - r; // don't change "sign" of NaNs } return r; } /* maximum ulp error = 1.49241 */ float my_sinf (float a) { float r; int i; a = a * 0.0f + a; // inf -> NaN r = trig_red_f (a, 117435.992f, &i); r = sinf_cosf_core (r, i); return r; } /* maximum ulp error = 1.49510 */ float my_cosf (float a) { float r; int i; a = a * 0.0f + a; // inf -> NaN r = trig_red_f (a, 71476.0625f, &i); r = sinf_cosf_core (r, i + 1); return r; } 

    作为曾经试图实现它的人,我感受到了你的痛苦(没有双关语意)。

    在尝试之前,请确保您很好地理解浮点运算何时是精确的,并且知道双重算术如何工作。

    除此之外,我最好的建议是看看其他聪明人做了什么:

    需要了解更多c/c++开发分享Payne Hanek算法在C中的实现,也可以关注C/ C++技术分享栏目—计算机技术网(www.ctvol.com)

      以上就是c/c++开发分享Payne Hanek算法在C中的实现相关内容,想了解更多C/C++开发(异常处理)及C/C++游戏开发关注(计算机技术网)。

      本文来自网络收集,不代表计算机技术网立场,如涉及侵权请点击右边联系管理员删除。

      如若转载,请注明出处:https://www.ctvol.com/c-cdevelopment/978859.html

      (0)
      上一篇 2021年12月12日
      下一篇 2021年12月12日

      精彩推荐