嵌套产品的准确求和

Accurate summation of nested products

本文关键字：求和 嵌套      更新时间：2023-10-16

我想减少以下计算中的数值浮点误差。

我有一个以下形式的等式：

b_3+w_3*(b_2+w_2*(b_1+w_1*(b_0+w_0)))

其中变量 w 表示范围[0,1]中的某个浮点数，b 表示范围[1,~1000000]中的浮点常量。 b与下标单调增加（尽管这可能并不重要）。当然，这可以扩展到任意数量的术语：

b_4+w_4*(c_3+w_3*(b_2+w_2*(b_1+w_1*(b_0+w_0))))

这可以递归定义为： 

func(x,n):
   if(n==MAX)
      return x
   else
      return func(b[n]+x*w[n],n+1)
func(1,0)

如果我正在做一个在线求和，我可以使用Kahan求和算法（Kahan 1965），或者其他几种方法之一，ala Higham 1993或McNamee 2004，来限制我的错误大小。如果我做在线重复产品，我可以使用某种转换技术将问题简化为求和。

事实上，我不确定如何处理这个特定问题。有没有人有想法（以及随之而来的引用）？

谢谢！

海厄姆 1993."浮点求和的准确性"。SIAM 科学计算杂志。

卡汉 1965。"实用语：关于减少截断错误的进一步评论"。中美洲共同市场。"10.1145/363707.363723"。

麦克纳米 2004."准确求和方法的比较"。西格萨姆公牛。"10.1145/980175.980177"。

您的计算看起来类似于 Horner 方案，只是每个阶段都使用不同的权重 w[i]，而不是单个变量 x。

有补偿霍纳方案的算法，我认为你可以根据你的目的进行调整。例如，请参阅以下文章中的定理 3 和算法 2。

P. Langlois，如何使用补偿霍纳算法确保忠实的多项式计算。 第18届IEEE计算机算术研讨会，2007年6月25日至27日，ARITH '07，第141-149页，http://www.acsel-lab.com/arithmetic/papers/ARITH18/ARITH18_Langlois.pdf

如果在算法 2 中将 TwoProd （s[i+1]， x） 替换为 TwoProd （s[i+1]， w[i+1]），似乎你会得到想要的结果，但我还没有尝试过。

您定义func的方式，它的计算结果为以下表达式：

For MAX = n+1, func(1,0) ==
     n     n
   ---  -----
1 + >     | |  w[j]
   /---   | | 
    i=0  j=n-i

因此，我解决总和的方法是：

double s = 0.0;
double a = 1.0;
for (int i = 1; i <= MAX; ++i) {
    a *= w[MAX-i];
    s += a;
}
return 1.0 + s;

即使我们将func x输入值视为变量，它也只影响最终项。但由于它的范围，您应该小心计算它。

double s = 0.0;
double a = 1.0;
double ax = x;
for (int i = 1; i < MAX; ++i) {
    a *= w[MAX-i];
    ax *= w[MAX-i];
    s += a;
}
ax *= w[0];
s += ax;
return 1.0 + s;
相关文章：