nju*_*ffa 22 c algorithm floating-point
为了以合理的精度简单有效地实现快速数学函数,多项式极小极大近似通常是选择的方法.Minimax近似通常使用Remez算法的变体生成.各种广泛使用的工具,如Maple和Mathematica,都具有内置功能.通常使用高精度算法计算生成的系数.众所周知,简单地将这些系数舍入到机器精度导致所得实施中的次优精度.
相反,人们搜索密切相关的系数集合,这些系数可以精确地表示为机器编号,以生成机器优化的近似值.两篇相关论文是:
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]上的反正切的简单多项式近似,使用Horner方案和FMA在IEEE-754单精度算法中进行评估.请参阅atan_poly()下面的C99代码中的功能.由于目前无法访问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;
}
下面的函数是一个忠实全面的落实arctan上[0, 1]:
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非常接近最近的浮点数.也就是说,它是计算中的常数,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在某个区间内.我使用Parma Polyhedra库来处理这些约束系统.floats 1到1e-8降序,并检查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.
编辑:我现在已经自动化了这个程序.下面的函数也是一个忠实全面的落实arctan上[0, 1]:
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;
}
我觉得这个代码甚至存在令人惊讶.对于这些附近的系数,r10由于该多项式在s接近时的缓慢收敛,您可以得到在实际算术中以几个ulps的数量级计算的多项式之间的距离和界限之间的界限1.我原本以为通过调整系数来预期舍入错误的行为方式基本上是"不可篡改的".
我思考了我在评论中收到的各种想法，并根据这些反馈进行了一些实验。最后，我决定精炼的启发式搜索是最好的前进方式。我现在已经设法将最大误差降低atanf_poly()到 1.01036 ulp，只有三个参数超过了我设定的 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 误差都接近统一，比如在 [0.75,1.0] 中是保守的。这允许对最大误差小于某个界限（例如 1.08 ulps）的有趣系数集进行快速扫描。然后，我可以详细和详尽地测试在该点锚定的启发式选择的超锥内的所有系数集。这第二步搜索最小 ulp 错误作为主要目标，并将正确舍入结果的最大百分比作为次要目标。通过在我的 CPU 的所有四个内核上使用这个两步过程，我能够显着加快搜索过程：到目前为止，我已经能够检查大约 2 21 个系数集。
基于所有“接近”解决方案中每个系数的范围，我现在估计这个近似问题的总有用搜索空间是 >= 2 24 个系数集，而不是我之前抛出的更乐观的 2 20 个。对于非常有耐心或拥有大量计算能力的人来说，这似乎是一个可行的问题。
我更新的代码如下：
// 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;
}