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

Вот две реализации функций интерполяции. Аргумент 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 double s все вычисления в 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(), но это касается моего вопроса.

Ответы

Ответ 1

Да, interol_80() безопасен, пусть демонстрирует его.

Проблема заключается в том, что входы имеют 64 байта float

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, тогда выполнение с неточными операциями float приведет к не более 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, поэтому интерполь-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 (big-large)/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

КЭД

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

ИЗМЕНИТЬ

Вот продолжение, поскольку следующее утверждение было немного приближенным и вызвало некоторые комментарии, когда 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

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

Ответ 2

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

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

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

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