Here are two implementations of interpolation functions. Argument u1
is always between 0.
and 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);
}
On a strict IEEE 754 platform with 80-bit long double
s, all computations in interpol_64()
are done according to IEEE 754 double precision, and in interpol_80()
in 80-bit extended precision.
The program prints:
u2: 0x1.8f5c28f5c28f6p-3
y64:0x1.8f5c28f5c28f5p-3
y80:0x1.8f5c28f5c28f6p-3
I am interested in the property “the result returned by the function is always in-between u2
and u3
”. This property is false of interpol_64()
, as shown by the values in the main()
above.
Does the property have a chance to be true of interpol_80()
? If it isn't, what is a counter-example? Does it help if we know that u2 != u3
or that there is a minimum distance between them? Is there a method to determine a significand width for intermediate computations at which the property would be guaranteed to be true?
EDIT: on all the random values I tried, the property held when intermediate computations were done in extended precision internally. If interpol_80()
took long double
arguments, it would be relatively easy to build a counter-example, but the question here is specifically about a function that takes double
arguments. This makes it much harder to build a counter-example, if there is one.
Note: a compiler generating x87 instructions may generate the same code for interpol_64()
and interpol_80()
, but this is tangential to my question.