Oto dwie implementacje funkcji interpolacji. Argumentu1 jest zawsze pomiędzy0. i1..

#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);
}

Na ścisłej platformie IEEE 754 z 80-bitowąlong doubles, wszystkie obliczenia winterpol_64() są wykonywane zgodnie z podwójną precyzją IEEE 754 i ininterpol_80() w 80-bitowej precyzji. Program drukuje:

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

Interesuje mnie właściwość „wynik zwracany przez funkcję jest zawsze pomiędzyu2 iu3” Ta właściwość jest fałszywainterpol_64(), jak pokazują wartości wmain() powyżej.

Czy nieruchomość ma szansę być prawdąinterpol_80()? Jeśli tak nie jest, jaki jest kontrprzykład? Czy to pomaga, jeśli to wiemyu2 != u3 czy istnieje między nimi minimalna odległość? Czy istnieje metoda określania szerokości istotności dla obliczeń pośrednich, w których właściwość byłaby gwarantowana jako prawdziwa?

EDYTOWAĆ: na wszystkich losowych wartościach, które wypróbowałem, właściwość utrzymywana, gdy wewnętrzne obliczenia były wykonywane z rozszerzoną precyzją wewnętrznie. Jeśliinterpol_80() wzięłalong double argumenty, względnie łatwo byłoby zbudować kontrprzykład, ale pytanie dotyczy konkretnie funkcji, która to przyjmujedouble argumenty. To znacznie utrudnia zbudowanie kontrprzykładu, jeśli taki istnieje.

Uwaga: kompilator generujący instrukcje x87 może wygenerować ten sam kod dlainterpol_64() iinterpol_80(), ale to jest styczne do mojego pytania.