Надежная и быстрая сумма квадратов с плавающей точкой - PullRequest
Новая
       8

Надежная и быстрая сумма квадратов с плавающей точкой

0 голосов
/ 13 января 2019

Я хотел бы создать надежную и быструю сумму квадратов . Простое решение будет (на 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] является маленьким или большим (вычисление в / переполнение). И это не самая точная процедура, так как мы можем лучше суммировать числа (например, сортировать перед суммированием).

Я проверил LAPACK DLASSQ, чтобы получить вдохновение (я добавил исходный код в конце моего вопроса).

По сути, DLASSQ поддерживает переменную масштаба (эта переменная рассчитывается на лету, это самый большой из уже увиденных элементов), а элементы массива масштабируются по этой переменной с делением. Таким образом, подпрограмма позволяет избежать недополнения / переполнения.

Однако мне не нравится это решение, так как:

  • есть, по крайней мере, n делений, что означает, что процедура медленная
  • эти деления могут снизить точность

Я думаю о возможном улучшении: шкала всегда будет степенью 2 (самый большой из уже увиденных элементов, округленный до степени 2). Таким образом:

  • Так как обратная величина числа 2 равна степени 2, я могу прекрасно хранить обратную величину шкалы
  • Таким образом, деление по шкале может быть заменено умножением (так что процедура будет намного быстрее)
  • Это умножение не потеряет никакой точности, потому что это всего лишь поправка на показатель степени

Мой вопрос: это хорошая идея? Есть ли ловушка, которую я не вижу?

Вот исходный код DLASSQ (код Фортрана): 

      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
...