Свойства 80-битных вычислений с расширенной точностью, начиная с аргументов двойной точности

Question

Вот две реализации интерполяционных функций. аргументация u1 всегда между 0. а также 1.,

#include <stdio.h>

double interpol_64(double u1, double u2, double u3)
{ 
  return u2 * (1.0 - u1) + u1 * u3;  
}

double interpol_80(double u1, double u2, double u3)
{ 
  return u2 * (1.0 - (long double)u1) + u1 * (long double)u3;  
}

int main()
{
  double y64,y80,u1,u2,u3;
  u1 = 0.025;
  u2 = 0.195;
  u3 = 0.195;
  y64 = interpol_64(u1, u2, u3);
  y80 = interpol_80(u1, u2, u3);
  printf("u2: %a\ny64:%a\ny80:%a\n", u2, y64, y80);
}

На строгой платформе IEEE 754 с 80-битным long doubles, все вычисления в interpol_64() сделаны в соответствии с двойной точностью IEEE 754, и в interpol_80() в расширенной точности 80 бит. Программа печатает:

u2: 0x1.8f5c28f5c28f6p-3
y64:0x1.8f5c28f5c28f5p-3
y80:0x1.8f5c28f5c28f6p-3

Меня интересует свойство "результат, возвращаемый функцией, всегда находится между u2 а также u3". Это свойство ложно interpol_64(), как показано значениями в main() выше.

Есть ли у собственности шанс быть правдой interpol_80()? Если это не так, что является контрпримером? Поможет ли нам это знать? u2 != u3 или что между ними есть минимальное расстояние? Существует ли метод определения значения и ширины для промежуточных вычислений, при которых это свойство гарантированно будет истинным?

РЕДАКТИРОВАТЬ: для всех случайных значений, которые я пробовал, свойство сохранялось, когда промежуточные вычисления выполнялись внутри с расширенной точностью. Если interpol_80() взял long double аргументы, было бы относительно легко построить контрпример, но вопрос здесь конкретно о функции, которая принимает double аргументы. Это усложняет создание контрпримеров, если таковые имеются.

---

Примечание: компилятор, генерирующий инструкции x87, может генерировать тот же код для interpol_64() а также interpol_80(), но это касается моего вопроса.

Answer 1

Да, interpol_80() безопасен, давайте продемонстрируем это.

Проблема утверждает, что входные данные являются 64-битными плавающими

rnd64(ui) = ui

Результат точно (при условии * и + математические операции)

r = u2*(1-u1)+(u1*u3)

Оптимальное возвращаемое значение с округлением до 64 битов с плавающей запятой

r64 = rnd64(r)

Как у нас есть эти свойства

u2 <= r <= u3

Гарантируется, что

rnd64(u2) <= rnd64(r) <= rnd64(u3)
u2 <= r64 <= u3

Преобразование в 80 битов u1,u2,u3 также является точным.

rnd80(ui)=ui

Теперь давайте предположим, 0 <= u2 <= u3затем выполнение с неточными операциями с плавающей запятой приводит к не более 4 ошибкам округления:

rf = rnd(rnd(u2*rnd(1-u1)) + rnd(u1*u3))

При условии округления до ближайшего четного значения это будет не более 2 ULP от точного значения. Если округление выполняется с плавающей запятой 64 бит или плавающей запятой 80 бит:

r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r)
r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r)

rf64 может быть отключено 2 ulp, так что interpol-64() небезопасен, но как насчет rnd64( rf80 )?
Мы можем сказать, что:

rnd64(r - 2 ulp80(r)) <= rnd64(rf80) <= rnd64(r + 2 ulp80(r))

поскольку 0 <= u2 <= u3, затем

ulp80(u2) <= ulp80(r) <= ulp80(r3)
rnd64(u2 - 2 ulp80(u2)) <= rnd64(r - 2 ulp80(r)) <= rnd64(rf80)
rnd64(u3 + 2 ulp80(u3)) >= rnd64(r + 2 ulp80(r)) >= rnd64(rf80)

К счастью, как и все числа в диапазоне (u2-ulp64(u2)/2 , u2+ulp64(u2)/2) мы получаем

rnd64(u2 - 2 ulp80(u2)) = u2
rnd64(u3 + 2 ulp80(u3)) = u3

поскольку ulp80(x)=ulp62(x)/2^(64-53)

Таким образом, мы получаем доказательство

u2 <= rnd64(rf80) <= u3

Для u2 <= u3 <= 0 мы можем легко применить это же доказательство.

Последний случай, который нужно изучить, это u2 <= 0 <= u3. Если мы вычтем 2 больших значения, то результат может быть вплоть до выключения ulp(большой)/2, а не ulp(большой-большой)/2...
Таким образом, это утверждение, которое мы сделали, больше не имеет места:

r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r)

К счастью, u2 <= u2*(1-u1) <= 0 <= u1*u3 <= u3 и это сохраняется после округления

u2 <= rnd(u2*rnd(1-u1)) <= 0 <= rnd(u1*u3) <= u3

Таким образом, поскольку добавленные количества имеют противоположный знак:

u2 <= rnd(u2*rnd(1-u1)) + rnd(u1*u3) <= u3

то же самое происходит после округления, так что мы можем еще раз гарантировать

u2 <= rnd64( rf80 ) <= u3

QED

Чтобы быть полными, мы должны позаботиться о ненормальных затратах (постепенный недостаток), но я надеюсь, что вы не будете так уж злы со стресс-тестами. Я не буду демонстрировать, что происходит с этими...

РЕДАКТИРОВАТЬ:

Вот продолжение, так как следующее утверждение было немного приблизительным и генерировало некоторые комментарии, когда 0 <= u2 <= u3

r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r)

Мы можем написать следующие неравенства:

rnd(1-u1) <= 1
rnd(1-u1) <= 1-u1+ulp(1)/4
u2*rnd(1-u1) <= u2 <= r
u2*rnd(1-u1) <= u2*(1-u1)+u2*ulp(1)/4
u2*ulp(1) < 2*ulp(u2) <= 2*ulp(r)
u2*rnd(1-u1) < u2*(1-u1)+ulp(r)/2

Для следующей операции округления мы используем

ulp(u2*rnd(1-u1)) <= ulp(r)
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(u2*rnd(1-u1))/2
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(r)/2
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)

Для второй части суммы имеем:

u1*u3 <= r
rnd(u1*u3) <= u1*u3 + ulp(u1*u3)/2
rnd(u1*u3) <= u1*u3 + ulp(r)/2

rnd(u2*rnd(1-u1))+rnd(u1*u3) < u2*(1-u1)+u1*u3 + 3*ulp(r)/2
rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 3*ulp(r)/2 + ulp(r+3*ulp(r)/2)/2
ulp(r+3*ulp(r)/2) <= 2*ulp(r)
rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 5*ulp(r)/2

Я не доказал первоначальное утверждение, но не так далеко...

Answer 2

Основным источником потери точности в interpol_64 это умножения. Умножение двух 53-битных мантисс позволяет получить 105- или 106-битную (в зависимости от того, несет ли старший бит) мантиссу. Это слишком большое значение, чтобы соответствовать 80-разрядному значению с расширенной точностью, поэтому в целом у вас также будет потеря точности в 80-разрядной версии. Определить, когда именно это происходит, очень сложно; самое простое сказать, что это происходит, когда накапливаются ошибки округления. Обратите внимание, что при добавлении двух терминов также есть небольшой шаг округления.

Большинство людей, вероятно, просто решат эту проблему с помощью такой функции:

double interpol_64(double u1, double u2, double u3)
{ 
  return u2 + u1 * (u3 - u2);
}

Но, похоже, вы ищете понимание вопросов округления, а не лучшей реализации.