boost :: math :: erf的算法

Question

有没有关于boost的erf函数背后算法的详细信息？该模块的文档不是很精确.我发现的是几种方法混合在一起.对我来说,它看起来像Abramowitz和Stegun的变种.

• 哪种方法混合？
• 这些方法是如何混合的？
• erf函数的复杂性是多少(恒定时间)？

塞巴斯蒂安

Answer 1

Boost Math Toolkit的文档有很多参考文献,其中包括Abramowitz和Stegun.的叔-功能接口包含一个策略可用于控制数值精度(以及因此它的运行时间复杂度)的模板的参数.

#include <boost/math/special_functions/erf.hpp>
namespace boost{ namespace math{

template <class T>
calculated-result-type erf(T z);

template <class T, class Policy>
calculated-result-type erf(T z, const Policy&);

template <class T>
calculated-result-type erfc(T z);

template <class T, class Policy>
calculated-result-type erfc(T z, const Policy&);

}} // namespaces

更新:

在下面提供的erf-function参考的"Implementation"部分的逐字副本下面:

履行

这些函数的所有版本首先使用通常的反射公式使其参数为正:

erf(-z) = 1 - erf(z);

erfc(-z) = 2 - erfc(z);  // preferred when -z < -0.5

erfc(-z) = 1 + erf(z);   // preferred when -0.5 <= -z < 0

这些函数的通用版本是根据不完整的伽玛函数实现的.

当识别有效数(尾数)大小时(当前用于53,64和113位实数,加上通过提升加倍的单精度24位处理),则使用由JM设计的一系列有理近似.

对于z <= 0.5,则使用erf的有理逼近,基于观察到erf是奇函数,因此erf使用以下公式计算:

erf(z) = z * (C + R(z*z));

其中有理逼近R(z*z)针对绝对误差进行了优化:只要其绝对误差与常数C相比足够小,那么在计算R(z*z)期间产生的任何舍入误差都将有效从结果中消失.结果,该区域中erf和erfc的误差非常低:最后一位在极少数情况下是不正确的.

对于z> 0.5,我们观察到在一小段时间[a,b],然后:

erfc(z) * exp(z*z) * z ~ c

对于某些常数c.

因此,对于z> 0.5,我们使用以下公式计算erfc:

erfc(z) = exp(-z*z) * (C + R(z - B)) / z;

R(z-B)再次针对绝对误差进行优化,并且常数C是在该范围的端点处取得的erfc(z)*exp(z*z)*z的平均值.再一次,只要R(z - B)中的绝对误差与c相比较小,那么c + R(z - B)将被正确舍入,结果中的误差将仅取决于exp的精度功能.实际上,除了极少数情况外,错误仅限于结果的最后一位.选择常数B使得有理逼近范围的左手端为0.

对于范围[a,+∞]上的大z,上述近似值被修改为:

erfc(z) = exp(-z*z) * (C + R(1 / z)) / z;

理性近似以难以理解的细节来解释.如果您需要更多详细信息,可以随时查看源代码.