关于算法：arctan如何实现？ |
您所在的位置：网站首页 › arctan怎么计算角度 › 关于算法：arctan如何实现？ |
|
该库的许多实现都深入到针对所有弧函数的FPATAN指导。 FPATAN如何实施？ 假设我们有1位符号，M位尾数和N位指数，那么获取该数字反正切的算法是什么？ 因为FPU会这样做，所以应该有这样的算法。 在x86处理器中FPATAN指令的实现通常是专有的。为了计算反正切函数或其他(反)三角函数，常用算法遵循三步过程： 用于将整个输入域映射到狭窄区间的参数减少 窄区间(一次近似区间)的核近似计算 基于参数约简的中间结果的扩展以产生最终结果参数减少通常基于众所周知的三角身份，可以在各种标准参考(例如MathWorld(http://mathworld.wolfram.com/InverseTangent.html))中进行查找。对于arctan的计算，常用的身份是 arctan(-x)= -arctan(x) arctan(1 / x)= 0.5 * pi-arctan(x)[x> 0] arctan(x)= arctan(c)+ arctan((x-c)/(1 + x * c))请注意，最后一个标识有助于构造值arctan(i / 2n)，i = 1 ... 2n的表，该表允许使用任意窄的主近似间隔，但要以增加表存储为代价。这是时空之间经典的编程折衷。 核心区间上的近似值通常是足够次数的极小多项式近似值。由于浮点除法的高成本，有理逼近通常在现代硬件上没有竞争力，并且由于两个多项式的计算加上除法带来的误差，有理逼近也遭受附加的数值误差。 通常使用Remez算法(http://en.wikipedia.org/wiki/Remez_algorithm)计算用于最小极大多项式逼近的系数。诸如Maple和Mathematica之类的工具具有内置工具来计算此类近似值。通过确保所有系数都是可精确表示的机器编号，可以提高多项式逼近的精度。我知道的唯一具有内置功能的工具是Sollya(http://sollya.gforge.inria.fr/)，它提供了fpminimax()函数。 多项式的评估通常利用有效且准确的霍纳方案(http://en.wikipedia.org/wiki/Horner%27s_method)或混合使用Estrin方案(http://en.wikipedia.org/wiki/ Estrin％27s_scheme)和霍纳氏症。 Estrin的方案允许人们很好地利用超标量处理器提供的指令级并行性，对整体指令数的影响很小，而对准确性的影响则常常(但不总是)良性的影响。 由于减少了舍入步骤的数量并提供了一些防止减法抵消的保护措施，因此使用FMA(融合乘加)可提高任一评估方案的准确性和性能。在许多处理器上都可以找到FMA，包括GPU和最近的x86 CPU。在标准C和标准C ++中，FMA操作作为fma()标准库函数公开，但是需要在不提供硬件支持的平台上进行仿真，这使其在这些平台上运行缓慢。 从编程的角度来看，当从文本表示形式转换为机器表示形式时，要避免将逼近和参数减少所需要的浮点常量转换时，希望避免转换错误的风险。 ASCII到浮点转换例程因包含棘手的错误而臭名昭著(例如http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/)。标准C(不是我所知的C ++，只能以专有扩展名提供)提供的一种机制是将浮点常量指定为直接表示基础位模式的十六进制文字，从而有效地避免了复杂的转换。 下面是计算双精度arctan()的C代码，该代码演示了上述许多设计原理和技术。这种快速构建的代码缺乏其他答案中指出的实现方式的复杂性，但应提供少于2 ul的错误结果，这在各种情况下都足够。我使用Remez算法的简单实现创建了一个自定义的minimax逼近，该算法对所有中间步骤都使用1024位浮点算法。我希望使用Sollya或类似工具可以产生数值上更好的近似值。 1234567891011121314151617181920212223242526272829303132333435double my_atan (double x) { double a, z, p, r, s, q, o; /* argument reduction: arctan (-x) = -arctan(x); arctan (1/x) = 1/2 * pi - arctan (x), when x > 0 */ z = fabs (x); a = (z > 1.0) ? 1.0 / z : z; /* evaluate minimax polynomial approximation */ s = a * a; // a**2 q = s * s; // a**4 o = q * q; // a**8 /* use Estrin's scheme for low-order terms */ p = fma (fma (fma (-0x1.53e1d2a25ff34p-16, s, 0x1.d3b63dbb65af4p-13), q, fma (-0x1.312788dde0801p-10, s, 0x1.f9690c82492dbp-9)), o, fma (fma (-0x1.2cf5aabc7cef3p-7, s, 0x1.162b0b2a3bfcep-6), q, fma (-0x1.a7256feb6fc5cp-6, s, 0x1.171560ce4a483p-5))); /* use Horner's scheme for high-order terms */ p = fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (p, s, -0x1.4f44d841450e1p-5), s, 0x1.7ee3d3f36bb94p-5), s, -0x1.ad32ae04a9fd1p-5), s, 0x1.e17813d66954fp-5), s, -0x1.11089ca9a5bcdp-4), s, 0x1.3b12b2db51738p-4), s, -0x1.745d022f8dc5cp-4), s, 0x1.c71c709dfe927p-4), s, -0x1.2492491fa1744p-3), s, 0x1.99999999840d2p-3), s, -0x1.555555555544cp-2) * s, a, a); /* back substitution based on argument reduction */ r = (z > 1.0) ? (0x1.921fb54442d18p+0 - p) : p; return copysign (r, x); }好。 相关讨论 出于好奇，是否存在使用弧度进行三角计算的情况比使用整数细分所能达到的精度高得多的情况？当然，使用以度，象限或整个圆为单位的角度，模量降低将更容易且更精确。我知道为什么弧度在微积分中很有用，但是无法精确表示整个圆的角度单位数似乎很棘手。 某些平台提供了sinpi()和cospi()函数，它们接受pi倍数的参数，这使参数简化变得容易。否则，很难精确地减少sin，cos，tan的参数，并且无论使用弧度还是度数，本质上都需要进行多精度中间计算。规范参考是：Mary H. Payne和Robert N. Hanek，《三角函数的弧度归约》，ACM SIGNUM通讯，第1卷。 18号1983年1月1日，第19-24页 有关度数参数约简的伴随论文是：Mary H. Payne和Robert N. Hanek，三角函数的度数约简，ACM SIGNUM通讯，第1卷。 18.不。 1983年4月2日，第18-19页 为什么在度数情况下需要多精度缩减？可以肯定的是，在乘数倍的情况下它更容易，但是fpmod(x，360.0)被指定为对于x的所有值都是绝对精确的，不是吗？顺便说一句，我不确定使用弧度时超精确参数减少有多大用处。如果尝试使用Math.Sin(x*2.0*Math.Pi)计算sin(2πx)，则如果对参数2.0*Math.Pi进行模数归约比对参数2x进行模数归一化，则结果将更为准确。 我可能对降级感到困惑(今天有点着急)。我从来不需要实施它，因此也没有考虑它。您所说的似乎是有道理的：如果可以使用本机精度的IEEE余数运算，那么这应该是精确缩减所需的全部。恕我直言，计算诸如sin(2πx)之类的术语的最佳解决方案是提供sinpi()函数，以便程序员可以编写sinpi(2*x)并获得与数学行为尽可能一致的结果。使用机器PI会引入??相位误差。 我不熟悉的用于计算触发函数的技术需要首先将角度转换为圆的2的幂次方。您知道不知道的任何技术吗？如果不是，您是否知道为什么sinpi等功能不应该普遍使用？让程序员将值按比例放大2pi，以便处理器可以按pi按比例缩小值，这似乎很疯狂。 我们正在远离主题(Stackoverflow并非为讨论而设计)。标准通常将现有用法编纂成文。 IEEE-754提到sinpi等作为推荐功能，一些C / C ++工具链将其作为扩展提供，而CUDA和OpenCL等GPU编程环境都包含了该功能。因此，如果程序员继续使用和要求它，那么我希望它成为标准的库函数，而且要持续几十年。 不赞成"对正弦，余弦，正切，棕褐色进行准确的参数减少很难...不管使用...或度数。reduced_degrees = fmod(raw_degrees, 360.0)是直接的范围减小。直接参考stackoverflow.com/questions/20928253/ @chux我同意按阶数简化trig函数参数很容易。不幸的是，当一个错误的发言时，除了在宽限期内，没有办法纠正评论。我建议使用remquo (angle,90.0)而不是fmod()。 同意remquo()更好-尽管我认为这是C99附加组件，但我在这里成功使用remquo()进行了改进，改进了sind()摘要：很难。另外，有时会在SO周围闲逛的Eric Postpischil和Stephen Canon也非常擅长。 许多特殊功能的常用方法如下： 处理NaN，无穷大和带符号的零作为特殊情况。 如果数字太大，结果将四舍五入到M_PI，则返回M_PI。将此阈值称为M。 如果存在某种形式的参数约简身份，请使用它来将参数带入更好的范围。 (这可能很棘手：对于sin和cos，这意味着您选择了2pi的精确值的倍数，以便落在正确的范围内。) 将[0,M)分成有限多个间隔。在每个间隔上使用Chebyshev近似法对较高阶的反正切。 (这是脱机完成的，通常是您在这些实现中看到的所有幻数的来源。此外，可以使用Remez的交换算法稍微加紧Chebyshev近似值，但是我不知道在任何情况下这样做有很大帮助) 找出参数所处的间隔(使用if和填充或仅使用表索引的技巧)，然后在该间隔上评估Chebyshev系列。这里特别需要一些属性： arctan实现应为单调；也就是说，如果x ，则arctan(x) 。 arctan实现应始终在正确答案的1 ulp之内返回答案。请注意，这是一个相对误差范围。评估切比雪夫级数以使这两个属性成立并不完全简单。这里使用两个double表示单个值的不同部分的技巧很常见。然后可能有一些案例工作表明该实现是单调的。同样，接近零时，对arctan的泰勒近似而不是Chebyshev近似-您处于相对误差范围之后，应该使用Horner规则评估级数。 如果您正在寻找要读取的atan实现，则fdlibm似乎不如glibc中的实现那么讨厌。参数减少似乎是基于触发标识tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a) tan(b))的，对tan(a)适当使用0.5，1或1.5。 相关讨论 既然我们在讨论这个主题，也许我应该在另一个问题中问这个问题，所以使用Pad近似而不是多项式的一个很好的理由是当近似函数(例如反正切)趋于+/-的有限极限时信息显然，度数大于1的多项式逼近永远不会有什么好处。现在我的问题是，由于我们无论如何都在进行参数约简，并且仅在[0…0.5]上使用过近似值，因此上述原因(我所听过的唯一原因)应该没有太大关系，应该是？ @PascalCuoq：我希望k的Chebyshev近似和总度(分子度+分母度)k的Pade-Chebyshev近似在紧凑区间上近似近似行为良好的函数。在没有这种减少参数的方案的情况下，我想您需要获得正确的度数差。 (我只需要编写特殊功能的低质量实现，因此在某些情况下使用合理逼近而不是多项式逼近可能会有更好的理由-我不知道。) 有理近似值很少有竞争力。浮点除法比FADD，FMUL或FMA要贵得多。同样，您必须处理两个多项式的误差加上除法的误差。在大多数情况下，您将需要直接多项式或表加多项式。就多项式而言，您需要针对目标精度优化系数，例如Sollyas fpminimax()函数提供的近似值。如果有FMA，它将有助于使评估误差较小。 Estrins方案可以帮助提高超标量体系结构的性能。三角函数确实有非常丑陋的实现，这些实现很hacky，并且有些麻烦。我认为在这里很难找到能够解释实际使用算法的人。 这是一个atan2实现：https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/e_atan2.c;h=a287ca6656b210c77367eec3c46d72f18476d61d;hb=HEAD 编辑：实际上，我找到了一个：http://www.netlib.org/fdlibm/e_atan2.c，它更容易理解，但可能会因为(？)而变慢。 FPU在某些电路中完成所有这些操作，因此CPU不必完成所有这些工作。 相关讨论 非常感谢。在第一个链接上，它还包含mpatan.h和mpatan.c，其中实现了atan-正是我想要的。 并非所有FPU都在硬件中执行此操作。可能有些架构没有三角指令。 SSE也不支持三角函数，因此MSVC 2013在矢量化代码时必须实现一种软件 x86 CPU中的FPATAN指令通常是通过微码实现的，即存储在处理器内部ROM中的一个小程序。尽管此类程序可能会使用可见ISA中不可用的专门操作，但通常不涉及任何特殊电路。 atan2的第二个实现要短得多，因为它使用了atan。 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |