优化的多项式极大极小近似反正切[
为了简单高效地实现具有合理精度的快速数学函数，多项式极大极小近似常常是选择的方法。 Minimax近似值通常由Remez算法的变体生成。 各种广泛使用的工具，如Maple和Mathematica都有内置的功能。 所生成的系数通常使用高精度算术来计算。 众所周知，简单地将这些系数舍入到机器精度会导致最终实现中的次最佳精度。
相反，人们搜索紧密相关的系数集合，这些系数可以精确地表示为机器编号，以生成机器优化的近似值。 两篇相关论文是：
Nicolas Brisebarre，Jean-Michel Muller和Arnaud Tisserand，“Computing Machine-Efficient Polynomial Approximations”，ACM Transactions on Mathematical Software，Vol。 第32卷，第2期，2006年6月，第236-256页。
Nicolas Brisebarre和Sylvain Chevillard，“高效多项式L∞逼近”，第18届IEEE计算机算术研讨会（ARITH-18），蒙彼利埃（法国），2007年6月，第169-176页。
后一篇文章中的LLL算法的实现可用作Sollya工具的fpminimax()命令。 我的理解是，所提出的用于生成机器优化近似的所有算法都基于启发式算法，因此一般不知道通过最优近似可达到什么精度。 目前还不清楚用于评估近似值的FMA（融合乘法 - 加法）的可用性是否会对该问题的答案产生影响。 天真地认为它应该。
我目前正在使用Horner方案和FMA在IEEE-754单精度算法中评估[-1,1]上的反正切的简单多项式近似。 请参阅下面的C99代码中的函数atan_poly() 。 由于目前缺乏对Linux机器的访问，我没有使用Sollya来生成这些系数，而是使用了我自己的启发式，它可以被宽松地描述为最陡的体面和模拟退火的混合体（以避免陷入局部最小值） 。 我的机器优化多项式的最大误差非常接近1 ulp，但理想情况下我希望最大ulp误差低于1 ulp。
我知道我可以改变我的计算来提高精度，例如通过使用超过单精度精度的超前系数，但我想保持代码完全按照原样（即尽可能简单）只调整系数以提供最准确的结果。
一个“经过验证的”最佳系数组是理想的，对相关文献的指引是受欢迎的。 我做了一个文献检索，但找不到任何有助于超越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]上进行如实的舍入并打印出“success”，则以下测试工具将中止：
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("successn");
exit(0);
}
在每次乘法中使用s的问题是多项式的系数不会快速衰减。 接近1的输入会导致大量和大量相同数字的取消，这意味着您正在尝试查找一组系数，以便计算结束时的累计舍入与arctan的残差接近。
(arctan(sqrt(x)) - x) / x^3非常接近最接近的浮点数，常量0x1.fde90cp-1f是一个接近于1的数字。 也就是说，这是一个常数，进入u的计算，以便三次系数几乎完全确定。 （对于此程序，立方系数必须是-0x1.b81b44p-3f或-0x1.b81b42p-3f 。）
通过交替乘法s和u具有减小舍入误差的影响的效果ri在r{i+2}至多1/4的一个因素，因为s*u < 1/4任何a是。 这给选择五阶及以上的系数带来了相当大的余地。
我通过两个程序找到了系数：
a ，可以计算r8的范围，导致一个忠实的四舍五入结果。 为了得到线性不等式，我假设r8将被计算为float s和u中的一个多项式，用于实数运算; 线性不等式约束了这个实数r8位于某个区间内。 我使用帕尔马多面体库来处理这些约束系统。 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 s值。 我曾经预料，舍入误差的行为方式只是通过调整系数而从根本上“不可撼动”。
这不是对这个问题的回答，但是太长以至于不适合评论：
您的问题是关于系数C3，C5，...，C17在多项式近似反正切的情况下的最佳选择，其中C1固定为1，C2，C4，...，C16固定为0。
问题的标题说明您正在寻找[-1,1]上的近似值，将偶数系数固定为0的一个很好的理由是，逼近恰好是一个奇函数是充分必要的。 您的问题中的代码通过仅在[0，1]上应用多项式近似来“与”标题相矛盾。
如果使用Remez算法来查找系数C2，C3，...，C8而不是[0，1]上的反正切函数的多项式近似值，则可能会得到类似下面的值的结果：
#include <stdio.h>
#include <math.h>
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: %fn%an%ann", x, atan_poly(x), atan(x));
}
这与你问题中的代码具有大致相同的复杂度 - 乘法的次数是相似的。 看看这个多项式，没有任何理由特别想要将任何系数固定为0.如果我们想在[-1,1]上近似一个奇数函数而不固定偶数系数，它们会自动变得非常小，以吸收，然后我们想将它们固定为0，但是对于[0,1]上的近似值，它们不是，所以我们不必固定它们。
它可能比你的问题中的奇数多项式更好或更差。 事实证明，情况更糟（见下文）。 这个LolRemez 0.2的快速和肮脏的应用（代码包含在问题的底部）似乎足以提出选择系数的问题。 特别要好奇的是，如果您将这个答案中的系数置于相同的“最陡下降和模拟退火的混合”优化步骤中，您将应用这些步骤来获得问题中的系数，会发生什么情况。
因此，为了总结这个评论 - 张贴为答案，你确定你正在寻找最佳系数C3，C5，...，C17？ 在我看来，你正在寻找能够产生忠实的反正切逼近的单精度浮点运算的最佳序列，并且这种近似不一定是17次奇次多项式的霍尔形式。
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的所有四个内核中使用这两步处理，我能够显着加快搜索过程：到目前为止，我已经能够检查大约221个系数集。
基于所有“接近”解决方案中每个系数的范围，我现在估计这个近似问题的总有用搜索空间大于或等于224个系数集，而不是之前抛出的更乐观的220个数。 这似乎是一个可行的问题，以解决谁是非常耐心的人或有大量的计算马力在他们的处置。
我的更新代码如下：
// 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;
}
上一篇: optimized polynomial minimax approximation to arctangent on [