ulps/ulp.h

   1 #ifndef __ULPH_H__
   2 #define __ULPH_H__
   3
   4 #include <assert.h>
   5
   6 struct ulp {
   7         double ulp_max;
   8         double ulp_min;
   9         double ulp_avg;
  10         size_t ulp_skipped;
  11
  12         long double ulp_maxl;
  13         long double ulp_minl;
  14         long double ulp_avgl;
  15         size_t ulp_skippedl;
  16 };
  17
  18 struct ulp_complex {
  19         struct ulp ulp_real;
  20         struct ulp ulp_imag;
  21         struct ulp ulp_total;
  22 };
  23
  24 /*******************************************************************************
  25  *                              Real arithmetic
  26  ******************************************************************************/
  27 #define ULP_INIT(u)                             \
  28         assert(u);                              \
  29         (u)->ulp_max      = -DBL_MAX;           \
  30         (u)->ulp_min      =  DBL_MAX;           \
  31         (u)->ulp_avg      =  0.0;               \
  32         (u)->ulp_skipped  =  0;                 \
  33         (u)->ulp_maxl     = -LDBL_MAX;          \
  34         (u)->ulp_minl     =  LDBL_MAX;          \
  35         (u)->ulp_avgl     =  0.0;               \
  36         (u)->ulp_skippedl =  0;
  37
  38 #define ULP_UPDATE(u, n)                        \
  39 assert(u);                                      \
  40 do {                                            \
  41         if ((n) > (u)->ulp_max)                 \
  42                   (u)->ulp_max = n;             \
  43         if ((n) < (u)->ulp_min)                 \
  44                   (u)->ulp_min = n;             \
  45         (u)->ulp_avg += n;                      \
  46 } while(0)
  47
  48 #define ULP_UPDATEL(u, nl)                      \
  49 assert(u);                                      \
  50 do {                                            \
  51         if ((nl) > (u)->ulp_maxl)               \
  52                    (u)->ulp_maxl = nl;          \
  53         if ((nl) < (u)->ulp_minl)               \
  54                    (u)->ulp_minl = nl;          \
  55         (u)->ulp_avgl += nl;                    \
  56 } while(0)
  57
  58
  59 /*******************************************************************************
  60  *                              Complex arithmetic
  61  ******************************************************************************/
  62 #define ULP_COMPLEX_INIT(uc)                                            \
  63         ULP_INIT(&uc->ulp_real);                                        \
  64         ULP_INIT(&uc->ulp_imag);                                        \
  65         ULP_INIT(&uc->ulp_total);
  66
  67 #define ULP_COMPLEX_UPDATE(uc, n)                                       \
  68 assert(uc);                                                             \
  69 ULP_UPDATE(&uc->ulp_real, creal(n));                                    \
  70 ULP_UPDATE(&uc->ulp_imag, cimag(n));                                    \
  71 ULP_UPDATE(&uc->ulp_total, fabs(creal(n)) + fabs(cimag(n)));
  72
  73 #define ULP_COMPLEX_UPDATEL(uc, n)                                      \
  74 assert(uc);                                                             \
  75 ULP_UPDATEL(&uc->ulp_real, creall(n));                                  \
  76 ULP_UPDATEL(&uc->ulp_imag, cimagl(n));                                  \
  77 ULP_UPDATEL(&uc->ulp_total, fabsl(creall(n)) + fabsl(cimagl(n)));
  78
  79 double calculp_double(double computed, double exact);
  80 long double calculp_long_double(long double computedl, long double exactl);
  81
  82 void getfunctionulp_real(const struct fentry *f, struct ulp *u);
  83 void getfunctionulp_complex(const struct fentry *f, struct ulp_complex *u);
  84
  85 void printulps_double(struct ulp u);
  86 void printulps_long_double(struct ulp u);
  87 void printulps_double_complex(struct ulp_complex u);
  88 void printulps_long_double_complex(struct ulp_complex u);
  89
  90 #define NITERATIONS (1 * 1000)
  91
  92 #endif  /* ! __ULP_H__ */