最佳的机器优化多项式极小极大近似于上的反正切？
为了以合理的精度简单有效地实现快速数学函数，多项式极小极大近似通常是选择的方法。 Minimax近似通常使用Remez算法的变体生成。 各种广泛使用的工具，如Maple和Mathematica，都具有内置function。 通常使用高精度算术来计算生成的系数。 众所周知，简单地将这些系数舍入到机器精度会导致所得实现中的次优精度。
相反，人们搜索密切相关的系数集合，这些系数可以精确地表示为机器编号，以生成机器优化的近似值。 两篇相关论文是：
Nicolas Brisebarre，Jean-Michel Muller和Arnaud Tisserand，“计算机高效多项式近似”，ACM数学软件交易，Vol。 32，第2期，2006年6月，第236-256页。
Nicolas Brisebarre和Sylvain Chevillard，“高效多项式L∞近似”，第18届IEEE计算机算术研讨会（ARITH-18），蒙彼利埃（法国），2007年6月，第169-176页。
来自后一篇论文的LLL算法的实现可以作为Sollya工具的fpminimax()命令获得。 我的理解是，所有用于生成机器优化近似的算法都是基于启发式算法，因此通常不知道通过最佳近似可以实现什么样的精度。 我不清楚用于评估近似的FMA（融合乘法 – 加法）的可用性是否对该问题的答案有影响。 在我看来它应该是天真的。
我目前正在研究[-1,1]上的反正切的简单多项式近似，其在IEEE-754单精度算法中使用Horner方案和FMA进行评估。 请参阅下面的C99代码中的函数atan_poly() 。 由于目前无法访问Linux机器，我没有使用Sollya来生成这些系数，而是使用了我自己的启发式算法，可以将其简单地描述为最陡的体面和模拟退火的混合（以避免陷入局部最小值） 。 机器优化多项式的最大误差非常接近1 ulp，但理想情况下我希望最大ulp误差低于1 ulp。
我知道我可以改变我的计算以提高精度，例如使用表示超过单精度精度的前导系数，但我想保持代码完全一样（即尽可能简单）仅调整系数以提供最准确的结果。
“经过validation的”最佳系数集将是理想的，欢迎指向相关文献的指针。 我进行了一次文献检索，但找不到任何有助于超越Sollya的fpminimax() ，并且没有任何文章能够检查FMA（如果有的话）在这个问题上的作用。
// max ulp err = 1.03143 float atan_poly (float a) { float r, s; s = a * a; r = 0x1.7ed1ccp-9f; r = fmaf (r, s, -0x1.0c2c08p-6f); r = fmaf (r, s, 0x1.61fdd0p-5f); r = fmaf (r, s, -0x1.3556b2p-4f); r = fmaf (r, s, 0x1.b4e128p-4f); r = fmaf (r, s, -0x1.230ad2p-3f); r = fmaf (r, s, 0x1.9978ecp-3f); r = fmaf (r, s, -0x1.5554dcp-2f); r = r * s; r = fmaf (r, a, a); return r; } // max ulp err = 1.52637 float my_atanf (float a) { float r, t; t = fabsf (a); r = t; if (t > 1.0f) { r = 1.0f / r; } r = atan_poly (r); if (t > 1.0f) { r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r } r = copysignf (r, a); return r; }
以下函数是[0, 1]上arctan的忠实舍入实现：
float atan_poly (float a) { float s = a * a, u = fmaf(a, -a, 0x1.fde90cp-1f); float r1 = 0x1.74dfb6p-9f; float r2 = fmaf (r1, u, 0x1.3a1c7cp-8f); float r3 = fmaf (r2, s, -0x1.7f24b6p-7f); float r4 = fmaf (r3, u, -0x1.eb3900p-7f); float r5 = fmaf (r4, s, 0x1.1ab95ap-5f); float r6 = fmaf (r5, u, 0x1.80e87cp-5f); float r7 = fmaf (r6, s, -0x1.e71aa4p-4f); float r8 = fmaf (r7, u, -0x1.b81b44p-3f); float r9 = r8 * s; float r10 = fmaf (r9, a, a); return r10; }
如果函数atan_poly未能在[1e-16, 1]上忠实地舍入并且打印“成功”，则以下测试工具将中止：
int checkit(float f) { double d = atan(f); float d1 = d, d2 = d; if (d1 < d) d2 = nextafterf(d1, 1.0/0.0); else d1 = nextafterf(d1, -1.0/0.0); float p = atan_poly(f); if (p != d1 && p != d2) return 0; return 1; } int main() { for (float f = 1; f > 1e-16; f = nextafterf(f, -1.0/0.0)) { if (!checkit(f)) abort(); } printf("success\n"); exit(0); }
在每次乘法中使用s的问题是多项式的系数不会快速衰减。 接近1的输入导致批次和大量取消几乎相等的数字，这意味着您试图找到一组系数，以便计算结束时的累积舍入非常接近arctan的残差。
常量0x1.fde90cp-1f是一个接近1的数字，其中(arctan(sqrt(x)) - x) / x^3非常接近最近的float。 也就是说，它是u的计算中的常数，因此几何完全确定了三次系数。 （对于此程序，立方系数必须为-0x1.b81b44p-3f或-0x1.b81b42p-3f 。）
通过s和u交替乘法具有将ri的舍入误差对r{i+2}的影响降低至多1/4的效果，因为s*u < 1/4无论是什么。 这为选择五阶及更高阶系数提供了相当大的余地。
我借助两个程序找到了系数：
- 一个程序插入一堆测试点，写下线性不等式系统，并计算不等式系统的系数界限。 请注意，给定
a，可以计算导致忠实舍入结果的r8范围。 为了得到线性不等式，我假装r8将被计算为floats中s多项式和实数运算中的u; 线性不等式将该实数r8约束在某个区间内。 我使用Parma Polyhedra库来处理这些约束系统。 - 另一个程序在一定范围内随机测试系数集，首先插入一组测试点，然后按降序插入
1到1e-8所有float，并检查atan_poly产生了atan((double)x)的忠实舍入。 如果某些x失败，它会打印出x以及失败的原因。
为了获得系数，我攻击了第一个程序来修复c3 ，在每个测试点的r7上计算出界限，然后在高阶系数上得到界限。 然后我将其修改为修复c3和c5并获得更高阶系数的界限。 我这样做直到除了三个最高阶系数c13 ， c15和c17我都有。
我在第二个程序中增加了一组测试点，直到它停止打印任何东西或打印出“成功”。 我几乎所有错误的多项式都需要很少的测试点 - 我在程序中计算了85个测试点。
在这里，我展示了一些选择系数的工作。 为了得到一个忠实的圆形arctan ，我的初始测试点集假设r1到r8是用实数算术计算的（并以某种方式不愉快地舍入，但是以一种我不记得的方式），但是r9和r10是用float算法计算的，我需要：
-0x1.b81b456625f15p-3 <= c3 <= -0x1.b81b416e22329p-3 -0x1.e71d48d9c2ca4p-4 <= c5 <= -0x1.e71783472f5d1p-4 0x1.80e063cb210f9p-5 <= c7 <= 0x1.80ed6efa0a369p-5 0x1.1a3925ea0c5a9p-5 <= c9 <= 0x1.1b3783f148ed8p-5 -0x1.ec6032f293143p-7 <= c11 <= -0x1.e928025d508p-7 -0x1.8c06e851e2255p-7 <= c13 <= -0x1.732b2d4677028p-7 0x1.2aff33d629371p-8 <= c15 <= 0x1.41e9bc01ae472p-8 0x1.1e22f3192fd1dp-9 <= c17 <= 0x1.d851520a087c2p-9
取c3 = -0x1.b81b44p-3，假设r8也在float运算中计算：
-0x1.e71df05b5ad56p-4 <= c5 <= -0x1.e7175823ce2a4p-4 0x1.80df529dd8b18p-5 <= c7 <= 0x1.80f00e8da7f58p-5 0x1.1a283503e1a97p-5 <= c9 <= 0x1.1b5ca5beeeefep-5 -0x1.ed2c7cd87f889p-7 <= c11 <= -0x1.e8c17789776cdp-7 -0x1.90759e6defc62p-7 <= c13 <= -0x1.7045e66924732p-7 0x1.27eb51edf324p-8 <= c15 <= 0x1.47cda0bb1f365p-8 0x1.f6c6b51c50b54p-10 <= c17 <= 0x1.003a00ace9a79p-8
取c5 = -0x1.e71aa4p-4，假设r7在float运算中完成：
0x1.80e3dcc972cb3p-5 <= c7 <= 0x1.80ed1cf56977fp-5 0x1.1aa005ff6a6f4p-5 <= c9 <= 0x1.1afce9904742p-5 -0x1.ec7cf2464a893p-7 <= c11 <= -0x1.e9d6f7039db61p-7 -0x1.8a2304daefa26p-7 <= c13 <= -0x1.7a2456ddec8b2p-7 0x1.2e7b48f595544p-8 <= c15 <= 0x1.44437896b7049p-8 0x1.396f76c06de2ep-9 <= c17 <= 0x1.e3bedf4ed606dp-9
取c7 = 0x1.80e87cp-5，假设r6在float运算中完成：
0x1.1aa86d25bb64fp-5 <= c9 <= 0x1.1aca48cd5caabp-5 -0x1.eb6311f6c29dcp-7 <= c11 <= -0x1.eaedb032dfc0cp-7 -0x1.81438f115cbbp-7 <= c13 <= -0x1.7c9a106629f06p-7 0x1.36d433f81a012p-8 <= c15 <= 0x1.3babb57bb55bap-8 0x1.5cb14e1d4247dp-9 <= c17 <= 0x1.84f1151303aedp-9
取c9 = 0x1.1ab95ap-5，假设r5在float运算中完成：
-0x1.eb51a3b03781dp-7 <= c11 <= -0x1.eb21431536e0dp-7 -0x1.7fcd84700f7cfp-7 <= c13 <= -0x1.7ee38ee4beb65p-7 0x1.390fa00abaaabp-8 <= c15 <= 0x1.3b100a7f5d3cep-8 0x1.6ff147e1fdeb4p-9 <= c17 <= 0x1.7ebfed3ab5f9bp-9
我为c11选择了接近中间范围的点，并随机选择了c13 ， c15和c17 。
编辑：我现在已经自动化了这个程序。 以下函数也是[0, 1]上arctan的忠实舍入实现：
float c5 = 0x1.997a72p-3; float c7 = -0x1.23176cp-3; float c9 = 0x1.b523c8p-4; float c11 = -0x1.358ff8p-4; float c13 = 0x1.61c5c2p-5; float c15 = -0x1.0b16e2p-6; float c17 = 0x1.7b422p-9; float juffa_poly (float a) { float s = a * a; float r1 = c17; float r2 = fmaf (r1, s, c15); float r3 = fmaf (r2, s, c13); float r4 = fmaf (r3, s, c11); float r5 = fmaf (r4, s, c9); float r6 = fmaf (r5, s, c7); float r7 = fmaf (r6, s, c5); float r8 = fmaf (r7, s, -0x1.5554dap-2f); float r9 = r8 * s; float r10 = fmaf (r9, a, a); return r10; }
令我惊讶的是，这个代码甚至存在。 对于这些附近的系数，由于当s接近1时该多项式的缓慢收敛，您可以得到r10与实数算术中评估的多项式的值之间的距离的界限，其大小为几个ulps。 我原本以为通过调整系数来预期舍入错误的行为方式基本上是“不可篡改的”。
这不是问题的答案，但是太长而不适合评论：
你的问题是关于系数C3，C5，…，C17在反正切的多项式近似中的最佳选择，其中你将C1固定为1，C2，C4，…，C16固定为0。
你的问题的标题说你正在寻找[-1,1]的近似值，并且将偶数系数固定为0的一个很好的理由是，近似精确到奇数函数就足够了。 通过仅在[0,1]上应用多项式逼近，您的问题中的代码“与”标题“矛盾”。
如果你使用Remez算法来寻找系数C2，C3，…，C8代替[0,1]上的反正切的多项式近似，你可能会得到类似下面的值：
#include #include float atan_poly (float a) { float r, s; s = a; // s = a * a; r = -3.3507930064626076153585890630056286726807491543578e-2; r = fmaf (r, s, 1.3859776280052980081098065189344699108643282883702e-1); r = fmaf (r, s, -1.8186361916440430105127602496688553126414578766147e-1); r = fmaf (r, s, -1.4583047494913656326643327729704639191810926020847e-2); r = fmaf (r, s, 2.1335202878219865228365738728594741358740655881373e-1); r = fmaf (r, s, -3.6801711826027841250774413728610805847547996647342e-3); r = fmaf (r, s, -3.3289852243978319173749528028057608377028846413080e-1); r = fmaf (r, s, -1.8631479933914856903459844359251948006605218562283e-5); r = fmaf (r, s, 1.2917291732886065585264586294461539492689296797761e-7); r = fmaf (r, a, a); return r; } int main() { for (float x = 0.0f; x < 1.0f; x+=0.1f) printf("x: %f\n%a\n%a\n\n", x, atan_poly(x), atan(x)); }
这与您问题中的代码大致相同 - 乘法的数量相似。 看看这个多项式，没有理由想要将任何系数固定为0.如果我们想要在[-1,1]上近似奇函数而不固定偶数系数，它们会自动出现非常小的主题吸收，然后我们想要将它们固定为0，但对于[0,1]的近似值，它们不会，所以我们不必固定它们。
它可能比你问题中的奇数多项更好或更差。 事实certificate它更糟糕（见下文）。 然而， LolRemez 0.2 （包含在问题底部的代码）的这种快速和肮脏的应用似乎足以提出系数选择的问题。 我特别好奇，如果你将这个答案中的系数置于同样的“最陡的体面和模拟退火的混合”优化步骤中，你应用于获得问题中的系数会发生什么。
所以，总结一下这个贴出来的答案，你确定你正在寻找最佳系数C3，C5，......，C17吗？ 在我看来，你正在寻找产生忠实近似于反正切的单精度浮点运算的最佳序列，并且这种近似不必是17阶多项式的Hornerforms。
x：0.000000 0x0p + 0 0x0p + 0 x：0.100000 0x1.983e2cp -4- 0x1.983e28938f9ecp -4- x：0.200000 0x1.94442p-3 0x1.94441ff1e8882p-3 x：0.300000 0x1.2a73a6p -2- 0x1.2a73a71dcec16p -2- x：0.400000 0x1.85a37ap -2- 0x1.85a3770ebe7aep -2- x：0.500000 0x1.dac67p -2- 0x1.dac670561bb5p -2- x：0.600000 0x1.14b1dcp-1 0x1.14b1ddf627649p-1 x：0.700000 0x1.38b116p-1 0x1.38b113eaa384ep-1 x：0.800000 0x1.5977a8p-1 0x1.5977a686e0ffbp-1 x：0.900000 0x1.773388p-1 0x1.77338c44f8faep-1
这是我链接到LolRemez 0.2的代码，以优化[0,1]上反正切的9次多项式逼近的相对精度：
#include "lol/math/real.h" #include "lol/math/remez.h" using lol::real; using lol::RemezSolver; real f(real const &y) { return (atan(y) - y) / y; } real g(real const &y) { return re (atan(y) / y); } int main(int argc, char **argv) { RemezSolver<8, real> solver; solver.Run("1e-1000", 1.0, f, g, 50); return 0; }
我思考了我在评论中收到的各种想法，并根据反馈进行了一些实验。 最后，我认为精确的启发式搜索是最好的前进方式。 我现在已经设法将atanf_poly()的最大错误减少到1.01036 ulps，只有三个参数超过了我声明的1 ulp错误界限的目标：
ulp = -1.00829 @ |a| = 9.80738342e-001 0x1.f62356p-1 (3f7b11ab) ulp = -1.01036 @ |a| = 9.87551928e-001 0x1.f9a068p-1 (3f7cd034) ulp = 1.00050 @ |a| = 9.99375939e-001 0x1.ffae34p-1 (3f7fd71a)
基于产生改进的近似的方式，不能保证这是最佳近似; 这里没有科学突破。 由于当前解决方案的ulp误差尚未完全平衡，并且由于继续搜索继续提供更好的近似值（尽管以指数方式增加时间间隔），我的猜测是可以实现1 ulp误差界限，但同时我们似乎已经非常接近最佳的机器优化近似。
新近似的更好质量是精确搜索过程的结果。 我观察到多项式中所有最大的ulp误差都接近于1，比如[0.75,1.0]是保守的。 这允许快速扫描有趣系数集，其最大误差小于某个界限，比如1.08 ulps。 然后，我可以详细测试并详尽地测试在该点锚定的启发式选择的超锥体内的所有系数集。 第二步搜索最小ulp错误作为主要目标，并将正确舍入结果的最大百分比作为次要目标。 通过在我的CPU的所有四个核心上使用这个两步过程，我能够显着加快搜索过程：到目前为止，我已经能够检查大约2 21个系数集。
基于所有“近似”解的每个系数的范围，我现在估计该近似问题的总有用搜索空间是> = 2 24个系数集，而不是之前抛出的2 20 I的更乐观数。 对于那些非常耐心或拥有大量计算能力的人来说，这似乎是一个可行的问题。
我的更新代码如下：
// max ulp err = 1.01036 float atanf_poly (float a) { float r, s; s = a * a; r = 0x1.7ed22cp-9f; r = fmaf (r, s, -0x1.0c2c2ep-6f); r = fmaf (r, s, 0x1.61fdf6p-5f); r = fmaf (r, s, -0x1.3556b4p-4f); r = fmaf (r, s, 0x1.b4e12ep-4f); r = fmaf (r, s, -0x1.230ae0p-3f); r = fmaf (r, s, 0x1.9978eep-3f); r = fmaf (r, s, -0x1.5554dap-2f); r = r * s; r = fmaf (r, a, a); return r; } // max ulp err = 1.51871 float my_atanf (float a) { float r, t; t = fabsf (a); r = t; if (t > 1.0f) { r = 1.0f / r; } r = atanf_poly (r); if (t > 1.0f) { r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r } r = copysignf (r, a); return r; }
更新 （在两年半后重新审视该问题后）
以T. Myklebust的草案出版物为起点，我发现[-1,1]上的反正切近似具有最小误差，最大误差为0.94528 ulp。
/* Based on: Tor Myklebust, "Computing accurate Horner form approximations to special functions in finite precision arithmetic", arXiv:1508.03211, August 2015. maximum ulp err = 0.94528 */ float atanf_poly (float a) { float r, s; s = a * a; r = 0x1.6d2086p-9f; // 2.78569828e-3 r = fmaf (r, s, -0x1.03f2ecp-6f); // -1.58660226e-2 r = fmaf (r, s, 0x1.5beebap-5f); // 4.24722321e-2 r = fmaf (r, s, -0x1.33194ep-4f); // -7.49753043e-2 r = fmaf (r, s, 0x1.b403a8p-4f); // 1.06448799e-1 r = fmaf (r, s, -0x1.22f5c2p-3f); // -1.42070308e-1 r = fmaf (r, s, 0x1.997748p-3f); // 1.99934542e-1 r = fmaf (r, s, -0x1.5554d8p-2f); // -3.33331466e-1 r = r * s; r = fmaf (r, a, a); return r; }
这不是一个答案，而是一个扩展的评论。
最近的Intel CPU和一些未来的AMD CPU都有AVX2。 在Linux中，在/proc/cpuinfo查找avx2标志以查看您的CPU是否支持这些。
AVX2是一个扩展，允许我们使用256位向量构建和计算 – 例如，八个单精度数，或四个双精度数 – 而不仅仅是标量。 它包括FMA3支持，意味着对这些向量进行融合乘法加法。 简单地说，AVX2允许我们并行评估八个多项，几乎与我们使用标量运算评估单个多项时相同。
函数error8()分析一组系数，使用预定义的x值，与atan(x)预先计算值进行比较，并返回ULP中的错误（分别低于和高于所需结果），以及结果数量完全匹配所需的浮点值。 这些不是简单地测试一组系数是否优于当前最知名的系数所需要的，而是允许不同的策略来测试哪些系数。 （基本上，ULP中的最大误差形成一个表面，我们试图找到该表面上的最低点;知道每个点的表面“高度”使我们能够对有关方向进行有根据的猜测 – – 如何更改系数。）
使用了四个预先计算的表： known_x用于参数， known_f用于正确舍入的单精度结果， known_a用于双精度“准确”值（我只是希望库atan()足够精确 – – 但是不应该在没有检查的情况下依赖它！），并且known_m将双精度差异缩放到ULP。 给定参数的期望范围， precalculate()函数将使用库atan()函数预先计算这些范围。 （它还依赖于IEEE-754浮点格式，浮点数和整数字节顺序相同，但在运行此代码的CPU上也是如此。）
注意， known_x ， known_f和known_a数组可以存储在二进制文件中; known_m内容通常来自known_a 。 使用库atan()而不validation它不是一个好主意 – 但因为我的匹配njuffa的结果，我没有费心去寻找更好的参考atan() 。
为简单起见，这里是示例程序forms的代码：
#define _POSIX_C_SOURCE 200809L #include #include #include #include #include #include /** poly8() - Compute eight polynomials in parallel. * @x - the arguments * @c - the coefficients. * * The first coefficients are for degree 17, the second * for degree 15, and so on, down to degree 3. * * The compiler should vectorize the expression using vfmaddXXXps * given an AVX2-capable CPU; for example, Intel Haswell, * Broadwell, Haswell E, Broadwell E, Skylake, or Cannonlake; * or AMD Excavator CPUs. Tested on Intel Core i5-4200U. * * Using GCC-4.8.2 and * gcc -O2 -march=core-avx2 -mtune=generic * this code produces assembly (AT&T syntax) * vmulps %ymm0, %ymm0, %ymm2 * vmovaps (%rdi), %ymm1 * vmovaps %ymm0, %ymm3 * vfmadd213ps 32(%rdi), %ymm2, %ymm1 * vfmadd213ps 64(%rdi), %ymm2, %ymm1 * vfmadd213ps 96(%rdi), %ymm2, %ymm1 * vfmadd213ps 128(%rdi), %ymm2, %ymm1 * vfmadd213ps 160(%rdi), %ymm2, %ymm1 * vfmadd213ps 192(%rdi), %ymm2, %ymm1 * vfmadd213ps 224(%rdi), %ymm2, %ymm1 * vmulps %ymm2, %ymm1, %ymm0 * vfmadd132ps %ymm3, %ymm3, %ymm0 * ret * if you omit the 'static inline'. */ static inline __v8sf poly8(const __v8sf x, const __v8sf *const c) { const __v8sf xx = x * x; return (((((((c[0]*xx + c[1])*xx + c[2])*xx + c[3])*xx + c[4])*xx + c[5])*xx + c[6])*xx + c[7])*xx*x + x; } /** error8() - Calculate maximum error in ULPs * @x - the arguments * @co - { C17, C15, C13, C11, C9, C7, C5, C3 } * @f - the correctly rounded results in single precision * @a - the expected results in double precision * @m - 16777216.0 raised to the same power of two as @a normalized * @n - number of vectors to test * @max_under - pointer to store the maximum underflow (negative, in ULPs) to * @max_over - pointer to store the maximum overflow (positive, in ULPs) to * Returns the number of correctly rounded float results, 0..8*n. */ size_t error8(const __v8sf *const x, const float *const co, const __v8sf *const f, const __v4df *const a, const __v4df *const m, const size_t n, float *const max_under, float *const max_over) { const __v8sf c[8] = { { co[0], co[0], co[0], co[0], co[0], co[0], co[0], co[0] }, { co[1], co[1], co[1], co[1], co[1], co[1], co[1], co[1] }, { co[2], co[2], co[2], co[2], co[2], co[2], co[2], co[2] }, { co[3], co[3], co[3], co[3], co[3], co[3], co[3], co[3] }, { co[4], co[4], co[4], co[4], co[4], co[4], co[4], co[4] }, { co[5], co[5], co[5], co[5], co[5], co[5], co[5], co[5] }, { co[6], co[6], co[6], co[6], co[6], co[6], co[6], co[6] }, { co[7], co[7], co[7], co[7], co[7], co[7], co[7], co[7] } }; __v4df min = { 0.0, 0.0, 0.0, 0.0 }; __v4df max = { 0.0, 0.0, 0.0, 0.0 }; __v8si eqs = { 0, 0, 0, 0, 0, 0, 0, 0 }; size_t i; for (i = 0; i < n; i++) { const __v8sf v = poly8(x[i], c); const __v4df d0 = { v[0], v[1], v[2], v[3] }; const __v4df d1 = { v[4], v[5], v[6], v[7] }; const __v4df err0 = (d0 - a[2*i+0]) * m[2*i+0]; const __v4df err1 = (d1 - a[2*i+1]) * m[2*i+1]; eqs -= (__v8si)_mm256_cmp_ps(v, f[i], _CMP_EQ_OQ); min = _mm256_min_pd(min, err0); max = _mm256_max_pd(max, err1); min = _mm256_min_pd(min, err1); max = _mm256_max_pd(max, err0); } if (max_under) { if (min[0] > min[1]) min[0] = min[1]; if (min[0] > min[2]) min[0] = min[2]; if (min[0] > min[3]) min[0] = min[3]; *max_under = min[0]; } if (max_over) { if (max[0] < max[1]) max[0] = max[1]; if (max[0] < max[2]) max[0] = max[2]; if (max[0] < max[3]) max[0] = max[3]; *max_over = max[0]; } return (size_t)((unsigned int)eqs[0]) + (size_t)((unsigned int)eqs[1]) + (size_t)((unsigned int)eqs[2]) + (size_t)((unsigned int)eqs[3]) + (size_t)((unsigned int)eqs[4]) + (size_t)((unsigned int)eqs[5]) + (size_t)((unsigned int)eqs[6]) + (size_t)((unsigned int)eqs[7]); } /** precalculate() - Allocate and precalculate tables for error8(). * @x0 - First argument to precalculate * @x1 - Last argument to precalculate * @xptr - Pointer to a __v8sf pointer for the arguments * @fptr - Pointer to a __v8sf pointer for the correctly rounded results * @aptr - Pointer to a __v4df pointer for the comparison results * @mptr - Pointer to a __v4df pointer for the difference multipliers * Returns the vector count if successful, * 0 with errno set otherwise. */ size_t precalculate(const float x0, const float x1, __v8sf **const xptr, __v8sf **const fptr, __v4df **const aptr, __v4df **const mptr) { const size_t align = 64; unsigned int i0, i1; size_t n, i, sbytes, dbytes; __v8sf *x = NULL; __v8sf *f = NULL; __v4df *a = NULL; __v4df *m = NULL; if (!xptr || !fptr || !aptr || !mptr) { errno = EINVAL; return (size_t)0; } memcpy(&i0, &x0, sizeof i0); memcpy(&i1, &x1, sizeof i1); i0 ^= (i0 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U; i1 ^= (i1 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U; if (i1 > i0) n = (((size_t)i1 - (size_t)i0) | (size_t)7) + (size_t)1; else if (i0 > i1) n = (((size_t)i0 - (size_t)i1) | (size_t)7) + (size_t)1; else { errno = EINVAL; return (size_t)0; } sbytes = n * sizeof (float); if (sbytes % align) sbytes += align - (sbytes % align); dbytes = n * sizeof (double); if (dbytes % align) dbytes += align - (dbytes % align); if (posix_memalign((void **)&x, align, sbytes)) { errno = ENOMEM; return (size_t)0; } if (posix_memalign((void **)&f, align, sbytes)) { free(x); errno = ENOMEM; return (size_t)0; } if (posix_memalign((void **)&a, align, dbytes)) { free(f); free(x); errno = ENOMEM; return (size_t)0; } if (posix_memalign((void **)&m, align, dbytes)) { free(a); free(f); free(x); errno = ENOMEM; return (size_t)0; } if (x1 > x0) { float *const xp = (float *)x; float curr = x0; for (i = 0; i < n; i++) { xp[i] = curr; curr = nextafterf(curr, HUGE_VALF); } i = n; while (i-->0 && xp[i] > x1) xp[i] = x1; } else { float *const xp = (float *)x; float curr = x0; for (i = 0; i < n; i++) { xp[i] = curr; curr = nextafterf(curr, -HUGE_VALF); } i = n; while (i-->0 && xp[i] < x1) xp[i] = x1; } { const float *const xp = (const float *)x; float *const fp = (float *)f; double *const ap = (double *)a; double *const mp = (double *)m; for (i = 0; i < n; i++) { const float curr = xp[i]; int temp; fp[i] = atanf(curr); ap[i] = atan((double)curr); (void)frexp(ap[i], &temp); mp[i] = ldexp(16777216.0, temp); } } *xptr = x; *fptr = f; *aptr = a; *mptr = m; errno = 0; return n/8; } static int parse_range(const char *const str, float *const range) { float fmin, fmax; char dummy; if (sscanf(str, " %f %f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %f:%f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %f,%f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %f/%f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff %ff %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff:%ff %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff,%ff %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff/%ff %c", &fmin, &fmax, &dummy) == 2) { if (range) { range[0] = fmin; range[1] = fmax; } return 0; } if (sscanf(str, " %f %c", &fmin, &dummy) == 1 || sscanf(str, " %ff %c", &fmin, &dummy) == 1) { if (range) { range[0] = fmin; range[1] = fmin; } return 0; } return errno = ENOENT; } static int fix_range(float *const f) { if (f && f[0] > f[1]) { const float tmp = f[0]; f[0] = f[1]; f[1] = tmp; } return f && isfinite(f[0]) && isfinite(f[1]) && (f[1] >= f[0]); } static const char *f2s(char *const buffer, const size_t size, const float value, const char *const invalid) { char format[32]; float parsed; int decimals, length; for (decimals = 0; decimals <= 16; decimals++) { length = snprintf(format, sizeof format, "%%.%df", decimals); if (length < 1 || length >= (int)sizeof format) break; length = snprintf(buffer, size, format, value); if (length < 1 || length >= (int)size) break; if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value) return buffer; decimals++; } for (decimals = 0; decimals <= 16; decimals++) { length = snprintf(format, sizeof format, "%%.%dg", decimals); if (length < 1 || length >= (int)sizeof format) break; length = snprintf(buffer, size, format, value); if (length < 1 || length >= (int)size) break; if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value) return buffer; decimals++; } length = snprintf(buffer, size, "%a", value); if (length < 1 || length >= (int)size) return invalid; if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value) return buffer; return invalid; } int main(int argc, char *argv[]) { float xrange[2] = { 0.75f, 1.00f }; float c17range[2], c15range[2], c13range[2], c11range[2]; float c9range[2], c7range[2], c5range[2], c3range[2]; float c[8]; __v8sf *known_x; __v8sf *known_f; __v4df *known_a; __v4df *known_m; size_t known_n; if (argc != 10 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) { fprintf(stderr, "\n"); fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]); fprintf(stderr, " %s C17 C15 C13 C11 C9 C7 C5 C3 x\n", argv[0]); fprintf(stderr, "\n"); fprintf(stderr, "Each of the coefficients can be a constant or a range,\n"); fprintf(stderr, "for example 0.25 or 0.75:1. x must be a non-empty range.\n"); fprintf(stderr, "\n"); return EXIT_FAILURE; } if (parse_range(argv[1], c17range) || !fix_range(c17range)) { fprintf(stderr, "%s: Invalid C17 range or constant.\n", argv[1]); return EXIT_FAILURE; } if (parse_range(argv[2], c15range) || !fix_range(c15range)) { fprintf(stderr, "%s: Invalid C15 range or constant.\n", argv[2]); return EXIT_FAILURE; } if (parse_range(argv[3], c13range) || !fix_range(c13range)) { fprintf(stderr, "%s: Invalid C13 range or constant.\n", argv[3]); return EXIT_FAILURE; } if (parse_range(argv[4], c11range) || !fix_range(c11range)) { fprintf(stderr, "%s: Invalid C11 range or constant.\n", argv[4]); return EXIT_FAILURE; } if (parse_range(argv[5], c9range) || !fix_range(c9range)) { fprintf(stderr, "%s: Invalid C9 range or constant.\n", argv[5]); return EXIT_FAILURE; } if (parse_range(argv[6], c7range) || !fix_range(c7range)) { fprintf(stderr, "%s: Invalid C7 range or constant.\n", argv[6]); return EXIT_FAILURE; } if (parse_range(argv[7], c5range) || !fix_range(c5range)) { fprintf(stderr, "%s: Invalid C5 range or constant.\n", argv[7]); return EXIT_FAILURE; } if (parse_range(argv[8], c3range) || !fix_range(c3range)) { fprintf(stderr, "%s: Invalid C3 range or constant.\n", argv[8]); return EXIT_FAILURE; } if (parse_range(argv[9], xrange) || xrange[0] == xrange[1] || !isfinite(xrange[0]) || !isfinite(xrange[1])) { fprintf(stderr, "%s: Invalid x range.\n", argv[9]); return EXIT_FAILURE; } known_n = precalculate(xrange[0], xrange[1], &known_x, &known_f, &known_a, &known_m); if (!known_n) { if (errno == ENOMEM) fprintf(stderr, "Not enough memory for precalculated tables.\n"); else fprintf(stderr, "Invalid (empty) x range.\n"); return EXIT_FAILURE; } fprintf(stderr, "Precalculated %lu arctangents to compare to.\n", 8UL * (unsigned long)known_n); fprintf(stderr, "\nC17 C15 C13 C11 C9 C7 C5 C3 max-ulps-under max-ulps-above correctly-rounded percentage cycles\n"); fflush(stderr); { const double percent = 12.5 / (double)known_n; size_t rounded; char c17buffer[64], c15buffer[64], c13buffer[64], c11buffer[64]; char c9buffer[64], c7buffer[64], c5buffer[64], c3buffer[64]; char minbuffer[64], maxbuffer[64]; float minulps, maxulps; unsigned long tsc_start, tsc_stop; for (c[0] = c17range[0]; c[0] <= c17range[1]; c[0] = nextafterf(c[0], HUGE_VALF)) for (c[1] = c15range[0]; c[1] <= c15range[1]; c[1] = nextafterf(c[1], HUGE_VALF)) for (c[2] = c13range[0]; c[2] <= c13range[1]; c[2] = nextafterf(c[2], HUGE_VALF)) for (c[3] = c11range[0]; c[3] <= c11range[1]; c[3] = nextafterf(c[3], HUGE_VALF)) for (c[4] = c9range[0]; c[4] <= c9range[1]; c[4] = nextafterf(c[4], HUGE_VALF)) for (c[5] = c7range[0]; c[5] <= c7range[1]; c[5] = nextafterf(c[5], HUGE_VALF)) for (c[6] = c5range[0]; c[6] <= c5range[1]; c[6] = nextafterf(c[6], HUGE_VALF)) for (c[7] = c3range[0]; c[7] <= c3range[1]; c[7] = nextafterf(c[7], HUGE_VALF)) { tsc_start = __builtin_ia32_rdtsc(); rounded = error8(known_x, c, known_f, known_a, known_m, known_n, &minulps, &maxulps); tsc_stop = __builtin_ia32_rdtsc(); printf("%-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %lu %.3f %lu\n", f2s(c17buffer, sizeof c17buffer, c[0], "?"), f2s(c15buffer, sizeof c15buffer, c[1], "?"), f2s(c13buffer, sizeof c13buffer, c[2], "?"), f2s(c11buffer, sizeof c11buffer, c[3], "?"), f2s(c9buffer, sizeof c9buffer, c[4], "?"), f2s(c7buffer, sizeof c7buffer, c[5], "?"), f2s(c5buffer, sizeof c5buffer, c[6], "?"), f2s(c3buffer, sizeof c3buffer, c[7], "?"), f2s(minbuffer, sizeof minbuffer, minulps, "?"), f2s(maxbuffer, sizeof maxbuffer, maxulps, "?"), rounded, (double)rounded * percent, (unsigned long)(tsc_stop - tsc_start)); fflush(stdout); } } return EXIT_SUCCESS; }
The code does compile using GCC-4.8.2 on Linux, but might have to be modified for other compilers and/or OSes. (I'd be happy to include/accept edits fixing those, though. I just don't have Windows or ICC myself so I could check.)
To compile this, I recommend
gcc -Wall -O3 -fomit-frame-pointer -march=native -mtune=native example.c -lm -o example
Run without arguments to see usage; 要么
./example 0x1.7ed24ap-9f -0x1.0c2c12p-6f 0x1.61fdd2p-5f -0x1.3556b0p-4f 0x1.b4e138p-4f -0x1.230ae2p-3f 0x1.9978eep-3f -0x1.5554dap-2f 0.75:1
to check what it reports for njuffa's coefficient set, compared against standard C library atan() function, with all possible x in [0.75, 1] considered.
Instead of a fixed coefficient, you can also use min:max to define a range to scan (scanning all unique single-precision floating-point values). Each possible combination of the coefficients is tested.
Because I prefer decimal notation, but need to keep the values exact, I use the f2s() function to display the floating-point values. It is a simple brute-force helper function, that uses the shortest formatting that yields the same value when parsed back to float.
例如，
./example 0x1.7ed248p-9f:0x1.7ed24cp-9f -0x1.0c2c10p-6f:-0x1.0c2c14p-6f 0x1.61fdd0p-5f:0x1.61fdd4p-5f -0x1.3556aep-4f:-0x1.3556b2p-4f 0x1.b4e136p-4f:0x1.b4e13ap-4f -0x1.230ae0p-3f:-0x1.230ae4p-3f 0x1.9978ecp-3f:0x1.9978f0p-3f -0x1.5554d8p-2f:-0x1.5554dcp-2f 0.75:1
computes all the 6561 (3 8 ) coefficient combinations ±1 ULP around njuffa's set for x in [0.75, 1]. (Indeed, it shows that decreasing C 17 by 1 ULP to 0x1.7ed248p-9f yields the exact same results.)
(That run took 90 seconds on Core i5-4200U at 2.6 GHz -- pretty much in line in my estimate of 30 coefficient sets per second per GHz per core. While this code is not threaded, the key functions are thread-safe, so threading should not be too difficult. This Core i5-4200U is a laptop, and gets pretty hot even when stressing just one core, so I didn't bother.)
(I consider the above code to be in public domain, or CC0-licensed where public domain dedication is not possible. In fact, I'm not sure if it is creative enough to be copyrightable at all. Anyway, feel free to use it anywhere in any way you wish, as long as you don't blame me if it breaks.)
有问题吗？ Enhancements? Edits to fix Linux/GCC'isms are welcome!