Вот две реализации функций интерполяции. Аргумент 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 все вычисления в interpol_64() выполняются в соответствии с двойной точностью IEEE 754, а в interpol_80() - с 80-битной расширенной точностью. Программа печатает:
u2: 0x1.8f5c28f5c28f6p-3
y64:0x1.8f5c28f5c28f5p-3
y80:0x1.8f5c28f5c28f6p-3
Меня интересует свойство «результат, возвращаемый функцией, всегда находится между u2 и u3». Это свойство имеет значение false для interpol_64(), как показано значениями в main() выше.
Есть ли у собственности шанс соответствовать interpol_80()? Если нет, то какой контрпример? Помогает ли нам знать, что u2 != u3 или что между ними существует минимальное расстояние? Есть ли метод определения ширины значимой для промежуточных вычислений, при которой свойство гарантированно будет истинным?
РЕДАКТИРОВАТЬ: для всех случайных значений, которые я пробовал, свойство сохранялось, когда промежуточные вычисления выполнялись с повышенной точностью внутренне. Если бы interpol_80() принимало long double аргументов, было бы относительно легко построить контрпример, но здесь вопрос конкретно касается функции, которая принимает double аргументов. Это значительно усложняет построение контрпримера, если он есть.
Примечание: компилятор, генерирующий инструкции x87, может генерировать один и тот же код для interpol_64() и interpol_80(), но это не относится к моему вопросу.
