为了以合理的精度简单有效地实现快速数学函数,多项式极小极大近似通常是选择的方法.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 = …我正在尝试将范围缩减作为实现正弦函数的第一步.
当使用0到20000的输入范围时,我得到的误差大到0.002339146.我的错误显然不应该那么大,我不知道如何减少它.我注意到误差幅度与输入θ幅度与余弦/正弦相关.
我能够获得本文提到的nearpi.c代码,但我不确定如何将代码用于单精度浮点.如果有人有兴趣,可以在以下链接找到nearpi.c文件:nearpi.c
这是我的MATLAB代码:
x = 0:0.1:20000;
% Perform range reduction
% Store constant 2/pi
twooverpi = single(2/pi);
% Compute y
y = (x.*twooverpi);
% Compute k (round to nearest integer
k = round(y);
% Solve for f
f = single(y-k);
% Solve for r
r = single(f*single(pi/2));
% Find last two bits of k
n = bitand(fi(k,1,32,0),fi(3,1,32,0));
n = single(n);
% Preallocate for speed
z(length(x)) = 0;
for i = 1:length(x)
switch(n(i)) …我已经在许多正弦/余弦的实现中看到了所谓的扩展模块化精度算法.但是它的用途是什么？例如,在cephes实现中,在缩减到[0,pi/4]范围后,他们正在进行这种模块化精度算法以提高精度.
下面的代码:
z = ((x - y * DP1) - y * DP2) - y * DP3;
其中DP1,DP2和DP3是一些硬编码系数.如何在数学上找到这些系数？我已经理解了大数字"模块化扩展算术"的目的,但这里的确切目的是什么？
我很难理解如何实现Payne和Hanek发布的范围缩减算法(三角函数的范围缩减)
我见过这个图书馆:http: //www.netlib.org/fdlibm/
但它看起来如此扭曲,我所创立的所有理论解释都太简单了,无法提供实现.
有一些好的......好的...很好的解释吗？