FreeBSD libm и точность функций

[netch80]

Программа:

#include <math.h>
#include <mpfr.h>
#include <stdio.h>

double fpu_sin(double x) {
  double r;
  r = sin(x);
  return r;
}

int main() 
{
  int i;
  int cd = 0;
  for (i = 0; i < 3000000; ++i) {
    double x = i * 1.0e-6;
    double r1 = fpu_sin(x);
    double r2;
    mpfr_t mx, mr2;
    mpfr_init2(mx, 53);
    mpfr_init2(mr2, 53);
    mpfr_set_d(mx, x, MPFR_RNDN);
    mpfr_sin(mr2, mx, MPFR_RNDN);
    r2 = mpfr_get_d(mr2, MPFR_RNDN);
    mpfr_clear(mr2);
    if (r1 != r2) {
      printf("Mismatched: x=%.20lg(%a)", x, x);
      printf(" r1=%.18lg(%a)", r1, r1);
      printf(" r2=%.18lg(%a)", r2, r2);
      printf("\n");
      if (++cd > 10) {
        break;
      }
    }
  }
  return cd ? 1 : 0;
}

На glibc 2.23 расхождений нет, программа выдаёт пустой выхлоп.

На FreeBSD (10.3, amd64) их много, вот первые несколько:

Mismatched: x=0.00260599999999999981(0x1.5592d98bf7f06p-9) r1=0.00260599705034083228(0x1.5592c0359972bp-9) r2=0.00260599705034083185(0x1.5592c0359972ap-9)
Mismatched: x=0.0098189999999999996(0x1.41bfbdf090f73p-7) r1=0.00981884222127722001(0x1.41be6b1ccbce4p-7) r2=0.00981884222127721827(0x1.41be6b1ccbce3p-7)
Mismatched: x=0.0187539999999999998(0x1.3343fa2ad3e92p-6) r1=0.0187529006832448421(0x1.333f5dc8f0659p-6) r2=0.0187529006832448386(0x1.333f5dc8f0658p-6)
Mismatched: x=0.0191819999999999977(0x1.3a4723aafff36p-6) r1=0.0191808236882919231(0x1.3a42349ce64acp-6) r2=0.0191808236882919196(0x1.3a42349ce64abp-6)

везде на 1 бит, но заметно.

Хм, я думал, что уж требования IEEE754 будут выполнять.
Или не верить MPFR? Это, кажется, уже перебор будет.

Comments

[netch80]

Есть отдельная numerics@. Но я его ещё не пересмотрел вглубь, там много умных слов ;)

Вообще же тема началась для меня, что в дискуссии напомнили давно забытое, что трансцедентные функции на FPU в принципе неточны (до 1.5ulp). Решил проверить на SSE, где они все софтовые, и ограничение на точность только в объёме используемого кода.