src/math/util.c

   1 #include <stdio.h>
   2 #include <stdint.h>
   3 #include "util.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 && signbit(got) == signbit(want))
  47                 return dwant;
  48         if (isinf(got)) {
  49                 got = copysignf(0x1p127, got);
  50                 want *= 0.5;
  51         }
  52         return scalbn(got - want, -eulpf(want)) + dwant;
  53 }
  54
  55 float ulperr(double got, double want, float dwant)
  56 {
  57         if (isnan(got) && isnan(want))
  58                 return 0;
  59         if (got == want && signbit(got) == signbit(want))
  60                 return dwant;
  61         if (isinf(got)) {
  62                 got = copysign(0x1p1023, got);
  63                 want *= 0.5;
  64         }
  65         return scalbn(got - want, -eulp(want)) + dwant;
  66 }
  67
  68 float ulperrl(long double got, long double want, float dwant)
  69 {
  70 #if LDBL_MANT_DIG == 53
  71         return ulperr(got, want, dwant);
  72 #elif LDBL_MANT_DIG == 64
  73         if (isnan(got) && isnan(want))
  74                 return 0;
  75         if (got == want && signbit(got) == signbit(want))
  76                 return dwant;
  77         if (isinf(got)) {
  78                 got = copysignl(0x1p16383L, got);
  79                 want *= 0.5;
  80         }
  81         return scalbn(got - want, -eulpl(want)) + dwant;
  82 #else
  83         // TODO
  84         return inf;
  85 #endif
  86 }
  87
  88 #define length(a) (sizeof(a)/sizeof*(a))
  89 #define flag(x) {x, #x}
  90 static struct {
  91         int flag;
  92         char *s;
  93 } eflags[] = {
  94         flag(INEXACT),
  95         flag(INVALID),
  96         flag(DIVBYZERO),
  97         flag(UNDERFLOW),
  98         flag(OVERFLOW)
  99 };
 100
 101 char *estr(int f)
 102 {
 103         static char buf[256];
 104         char *p = buf;
 105         int i, all = 0;
 106
 107         for (i = 0; i < length(eflags); i++)
 108                 if (f & eflags[i].flag) {
 109                         p += sprintf(p, "%s%s", all ? "|" : "", eflags[i].s);
 110                         all |= eflags[i].flag;
 111                 }
 112         if (all != f) {
 113                 p += sprintf(p, "%s%d", all ? "|" : "", f & ~all);
 114                 all = f;
 115         }
 116         p += sprintf(p, "%s", all ? "" : "0");
 117         return buf;
 118 }
 119
 120 char *rstr(int r)
 121 {
 122         switch (r) {
 123         case RN: return "RN";
 124         case RZ: return "RZ";
 125         case RU: return "RU";
 126         case RD: return "RD";
 127         }
 128         return "R?";
 129 }
 130
 131 void setupfenv(int r)
 132 {
 133         fesetround(r);
 134         feclearexcept(FE_ALL_EXCEPT);
 135 }
 136
 137 int getexcept(void)
 138 {
 139         return fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW);
 140 }
 141