对[-1,1]上反正切的最佳机器优化多项式极小极大逼近？

Question

多项式极小极大逼近是快速数学函数的一种简单有效的实现方法，具有合理的精度。Minimax近似通常是用Remez算法的一个变体生成的。各种广泛可用的工具，如Maple和Mathematica，都有这方面的内置功能。生成的系数通常使用高精度算法计算。众所周知，简单地将这些系数舍入到机器精度会导致最终实现的精度不佳。

取而代之的是，人们搜索密切相关的系数集，这些系数集可以精确地表示为机器数，以生成机器优化的近似值。两篇相关论文是：

Nicolas Brisebarre,Jean-Michel Muller,Arnaud Tisserand,“计算机机器效率多项式逼近”，ACM数学软件交易，卷。32,第2号，2006年6月，第236-256页。

我目前正在研究一个简单的多项式逼近[-1,1]上的反正切，是在IEEE-754单精度算法中，使用Horner格式和FMA计算的。请参见下面C99代码中的函数atan_poly()。由于目前无法访问Linux机器，我没有使用Sollya来生成这些系数，而是使用了我自己的启发式，可以粗略地描述为最陡体面和模拟退火的混合（以避免陷入局部极小值）。我的机器优化多项式的最大误差非常接近1ulp，但理想情况下我希望最大ulp误差低于1ulp。

我知道我可以改变我的计算来提高精度，例如通过使用一个表示超过单精度精度的前导系数，但我希望保持代码完全原样（即尽可能简单），只调整系数以提供尽可能精确的结果。

一个“被证明的”最优系数集将是理想的，相关文献的指针是受欢迎的。我做了一个文献检索，但找不到任何一篇论文在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;
}

Answer 1

以下函数是[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("success\n");
  exit(0);
}

在每次乘法中使用s的问题是多项式的系数不会迅速衰减。接近1的输入会导致大量几乎相等的数字的抵消，这意味着您试图找到一组系数，以便计算结束时的累计舍入接近arctan的残差。

我通过两个程序找到了系数：

• 一个程序插入一堆测试点，写下一个线性不等式系统，并从该不等式系统中计算系数的界。注意，给定，可以计算R8的范围，从而得到忠实的舍入结果。为了得到线性不等式，我假设R8将在实数算术中以floatSS和U中的多项式计算；线性不等式约束这个实数r8位于某个区间。我使用Parma Polyhedra库来处理这些约束系统。
• 另一个程序随机测试某些范围内的系数集，首先插入一组测试点，然后按降序插入从1到1e-8的所有float，并检查atan_poly是否产生了atan(double)x)的忠实舍入。如果某个X失败，它会打印出X以及失败的原因。

为了获得系数，我黑了第一个程序来修复C3，为每个测试点计算R7的边界，然后获得高阶系数的边界。然后我对它进行了黑客攻击，以修复C3和C5并获得高阶系数的边界。我一直这样做，直到获得了除三个最高阶系数C13、C15和C17之外的所有系数。

我在第二个程序中增加了测试点集，直到它停止打印任何内容或打印出“成功”。我需要很少的测试点来拒绝几乎所有错误的多项式--我在程序中计算了85个测试点。

在这里，我展示了一些我的工作，选择系数。假设R1到R8是用实算术求取的，但R9和R10是用float算术求取的，为了得到对初始测试点集的忠实舍入的arctan，我需要：

-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.b81b44p3，假设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;
}