Точность c_k = a + ( N + k) * b

Question

a, b - 32-разрядные значения с плавающей запятой, N - 32-разрядное целое число, а k может принимать значения 0, 1, 2,... M. Необходимо вычислить c_k = a + (N + k) * b; Операции должны быть 32-битными (не с двойной точностью). Проблема заключается в точности - что из следующего является более точным?:

I) c_k = a + (N + k) * b

II) сначала рассчитать: c_0 = a + N * b
Затем рассчитайте c_1, c_2 и т. Д. Итеративно путем сложения:
c_1 = c_0 + b;
c_2 = c_1 + b;

Answer 1

Цепное сложение является одной из худших операций, которые вы можете сделать, потому что ошибка округления в последнем результате будет чистой суммой ошибки округления одной операции для каждого сложения в цепочке. Было бы точнее либо использовать первый способ, либо использовать c_i = c_0 + b*i,

Answer 2

Поскольку вам кажется, что вас не волнует количество операций, при условии модели IEEE 754 вы можете выполнить ее точно с 32-битными операциями.
См. Арифметические и быстрые робастные геометрические предикаты с плавающей точкой Шевчука - http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf или http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

Вы определяете две точные операции (см. Статью)

(product,residue) = twoproduct(a,b)
(sum,residue) = twosum(a,b)

Затем вы должны разложить N+k на два 24-битных значения, например

NkH = (N+k) / 256;
NkL = (N+K) % 256;

Тогда у вас есть два потенциально неточных умножения

( HH , HL ) = twoproduct( NkH , b)
( LH , LL ) = twoproduct( NkL , b)

Затем вы можете суммировать эти ( HH, HL) + ( LH, LL) + a

Это может быть выполнено точно с помощью суммы быстрого расширения (см. Статью снова)

(c1,c2,c3,c4,c5) = sort_increasing_magnitude(HH,HL,LH,LL,a)
(s2,s1) = twosum( c2,c1 )
(s3,s2) = twosum( c3,s2 )
(s4,s3) = twosum( c4,s3 )
(s5,s4) = twosum( c5,s4 )

Затем вы получите точно округленный результат в s5, как если бы операции выполнялись с арифметикой с бесконечной точностью.