src/common/mtest.c

   1 #include <stdio.h>
   2 #include <stdint.h>
   3 #include "mtest.h"
   4
   5 int eulpf(float x)
   6 {
   7         union { float f; uint32_t i; } u = { x };
   8         int e = u.i>>23 & 0xff;
   9
  10         if (!e)
  11                 e++;
  12         return e - 0x7f - 23;
  13 }
  14
  15 int eulp(double x)
  16 {
  17         union { double f; uint64_t i; } u = { x };
  18         int e = u.i>>52 & 0x7ff;
  19
  20         if (!e)
  21                 e++;
  22         return e - 0x3ff - 52;
  23 }
  24
  25 int eulpl(long double x)
  26 {
  27 #if LDBL_MANT_DIG == 53
  28         return eulp(x);
  29 #elif LDBL_MANT_DIG == 64
  30         union { long double f; struct {uint64_t m; uint16_t e; uint16_t pad;} i; } u = { x };
  31         int e = u.i.e & 0x7fff;
  32
  33         if (!e)
  34                 e++;
  35         return e - 0x3fff - 63;
  36 #else
  37         // TODO
  38         return 0;
  39 #endif
  40 }
  41
  42 float ulperrf(float got, float want, float dwant)
  43 {
  44         if (isnan(got) && isnan(want))
  45                 return 0;
  46         if (got == want) {
  47                 if (signbit(got) == signbit(want))
  48                         return dwant;
  49                 return inf;
  50         }
  51         if (isinf(got)) {
  52                 got = copysignf(0x1p127, got);
  53                 want *= 0.5;
  54         }
  55         return scalbn(got - want, -eulpf(want)) + dwant;
  56 }
  57
  58 float ulperr(double got, double want, float dwant)
  59 {
  60         if (isnan(got) && isnan(want))
  61                 return 0;
  62         if (got == want) {
  63                 if (signbit(got) == signbit(want))
  64                         return dwant;
  65                 return inf; // treat 0 sign errors badly
  66         }
  67         if (isinf(got)) {
  68                 got = copysign(0x1p1023, got);
  69                 want *= 0.5;
  70         }
  71         return scalbn(got - want, -eulp(want)) + dwant;
  72 }
  73
  74 float ulperrl(long double got, long double want, float dwant)
  75 {
  76 #if LDBL_MANT_DIG == 53
  77         return ulperr(got, want, dwant);
  78 #elif LDBL_MANT_DIG == 64
  79         if (isnan(got) && isnan(want))
  80                 return 0;
  81         if (got == want) {
  82                 if (signbit(got) == signbit(want))
  83                         return dwant;
  84                 return inf;
  85         }
  86         if (isinf(got)) {
  87                 got = copysignl(0x1p16383L, got);
  88                 want *= 0.5;
  89         }
  90         return scalbnl(got - want, -eulpl(want)) + dwant;
  91 #else
  92         // TODO
  93         return inf;
  94 #endif
  95 }
  96
  97 #define length(a) (sizeof(a)/sizeof*(a))
  98 #define flag(x) {x, #x}
  99 static struct {
 100         int flag;
 101         char *s;
 102 } eflags[] = {
 103         flag(INEXACT),
 104         flag(INVALID),
 105         flag(DIVBYZERO),
 106         flag(UNDERFLOW),
 107         flag(OVERFLOW)
 108 };
 109
 110 char *estr(int f)
 111 {
 112         static char buf[256];
 113         char *p = buf;
 114         int i, all = 0;
 115
 116         for (i = 0; i < length(eflags); i++)
 117                 if (f & eflags[i].flag) {
 118                         p += sprintf(p, "%s%s", all ? "|" : "", eflags[i].s);
 119                         all |= eflags[i].flag;
 120                 }
 121         if (all != f) {
 122                 p += sprintf(p, "%s%d", all ? "|" : "", f & ~all);
 123                 all = f;
 124         }
 125         p += sprintf(p, "%s", all ? "" : "0");
 126         return buf;
 127 }
 128
 129 char *rstr(int r)
 130 {
 131         switch (r) {
 132         case RN: return "RN";
 133 #ifdef FE_TOWARDZERO
 134         case RZ: return "RZ";
 135 #endif
 136 #ifdef FE_UPWARD
 137         case RU: return "RU";
 138 #endif
 139 #ifdef FE_DOWNWARD
 140         case RD: return "RD";
 141 #endif
 142         }
 143         return "R?";
 144 }