稳健快速浮点平方和

我想创造一个强大和快速 . 一个简单的解决方案是(在C++中):

double sumOfSquares(const double *v, size_t n) {
    double r = 0;
    for (size_t i=0; i<n; i++) {
        r += v[i]*v[i];
    }
    return r;
}

这个程序不健壮,因为当 v[i] 小的或大的(在/溢出下计算)。而且它也不是最精确的程序,因为我们可以用更好的方式对数字求和(比如在求和之前排序)。

DLASSQ 为了获得灵感(我在问题的最后加入了源代码)。

德拉斯克 维护一个scale变量(这个变量是动态计算的,它是已经看到的最大的元素),数组元素由这个变量进行除法缩放。这样例程就避免了不足/溢出。

• 有 n
• 这些部门 可以

我在考虑一个可能的改进:规模永远是2的幂(已经看到的最大元素四舍五入为2的幂)。这种方式:

• 因为2次方的倒数是2次方,所以我可以完美地存储刻度的倒数
• 所以按比例除法可以用乘法来代替(所以程序会快得多)
• 这个乘法不会失去任何精度,因为它只是一个指数调整

我的问题是:这是个好主意吗?有没有我没看到的陷阱?

这是 (fortran代码):

      SUBROUTINE DLASSQ( N, X, INCX, SCALE, SUMSQ )
*
*  -- LAPACK auxiliary routine (version 3.7.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     December 2016
*
*     .. Scalar Arguments ..
      INTEGER            INCX, N
      DOUBLE PRECISION   SCALE, SUMSQ
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   X( * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            IX
      DOUBLE PRECISION   ABSXI
*     ..
*     .. External Functions ..
      LOGICAL            DISNAN
      EXTERNAL           DISNAN
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS
*     ..
*     .. Executable Statements ..
*
      IF( N.GT.0 ) THEN
         DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX
            ABSXI = ABS( X( IX ) )
            IF( ABSXI.GT.ZERO.OR.DISNAN( ABSXI ) ) THEN
               IF( SCALE.LT.ABSXI ) THEN
                  SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2
                  SCALE = ABSXI
               ELSE
                  SUMSQ = SUMSQ + ( ABSXI / SCALE )**2
               END IF
            END IF
   10    CONTINUE
      END IF
      RETURN
*
*     End of DLASSQ
*
      END