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

Question

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

I) c_k = а + (N + K) * б

II), первый подсчет: c_0 = а + Н * Ь
Затем вычислите c_1, c_2 и т. Д. Итеративно путем добавления:
c_1 = c_0 + b;
c_2 = c_1 + b;

Answer 1

Поскольку вы, как представляется, не заботитесь о количестве операций, предполагая модель IEEE 754, вы можете выполнить ее точно с помощью операций с 32 битами.
См Shewchuck Адаптивная Точность арифметики с плавающей точкой и Fast Прочные Геометрические Предикаты - 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 + к надвое 24 бит significands, например

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

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

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

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

Это может быть выполнено точно с быстрой расширительной суммой (см бумаги снова)

(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, как если бы операции выполнялись с арифметикой с бесконечной точностью.

Answer 2

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