ulps/ulp_real.c

   1 #define _XOPEN_SOURCE 600
   2
   3 #include <assert.h>
   4 #include <float.h>
   5 #include <math.h>
   6 #include <stdio.h>
   7 #include <string.h>
   8
   9 #include <gmp.h>
  10 #include <mpfr.h>
  11
  12 #include "gen.h"
  13 #include "subr_random.h"
  14 #include "ulp.h"
  15
  16 double
  17 calculp_double(double computed, double exact)
  18 {
  19         double xbefore, xafter;
  20         double ulps;
  21
  22         xbefore = nextafter(exact, -INFINITY);
  23         xafter = nextafter(exact, +INFINITY);
  24         assert(xafter != xbefore);
  25
  26         ulps = fabs( (exact - computed) / (xafter - xbefore) );
  27
  28         return (ulps);
  29 }
  30
  31 long double
  32 calculp_long_double(long double computedl, long double exactl)
  33 {
  34         long double xbeforel, xafterl;
  35         long double ulpsl;
  36
  37         xbeforel = nextafterl(exactl, -INFINITY);
  38         xafterl = nextafterl(exactl, +INFINITY);
  39         assert(xafterl != xbeforel);
  40
  41         ulpsl = fabsl( (exactl - computedl) / (xafterl - xbeforel) );
  42
  43         return (ulpsl);
  44 }
  45
  46
  47
  48 static void
  49 populate_vars(const struct fentry *f,
  50     double *x, double *y, long double *xl, long double *yl,
  51     mpfr_ptr mp_x, mpfr_ptr mp_y, mpfr_ptr mp_xl, mpfr_ptr mp_yl)
  52 {
  53         const mpfr_rnd_t tonearest = GMP_RNDN;
  54         long double txl, tyl;
  55         double tx, ty;
  56
  57         assert(f && x && y && xl && yl);
  58         txl = tx = 0.0;
  59         tyl = ty = 0.0;
  60
  61         /* Generate random arguments */
  62         if (f->f_narg == 1) {
  63                 do {
  64                         tx  = random_double(FP_NORMAL);
  65                         txl = random_long_double(FP_NORMAL);
  66                 } while (!f->f_u.fp1(tx) || !f->f_u.fp1(txl));
  67         }
  68         if (f->f_narg == 2) {
  69                 do {
  70                         tx  = random_double(FP_NORMAL);
  71                         ty  = random_double(FP_NORMAL);
  72                         txl = random_long_double(FP_NORMAL);
  73                         tyl = random_long_double(FP_NORMAL);
  74                 } while (!f->f_u.fp2(tx, ty) || !f->f_u.fp2(txl, tyl));
  75         }
  76
  77         /* Hack, yikes */
  78         if (!strcmp(f->f_name, "yn")) {
  79                 txl = tx = 1.0;
  80         }
  81
  82         *x  = tx;
  83         *y  = ty;
  84         *xl = txl;
  85         *yl = tyl;
  86
  87         /* Copy arguments to mpfr variables */
  88         mpfr_set_d (mp_x,  tx,  tonearest);
  89         mpfr_set_d (mp_y,  ty,  tonearest);
  90         mpfr_set_ld(mp_xl, txl, tonearest);
  91         mpfr_set_ld(mp_yl, tyl, tonearest);
  92 }
  93
  94 void
  95 getfunctionulp_real(const struct fentry *f, struct ulp *u)
  96 {
  97         const mpfr_rnd_t tonearest = GMP_RNDN;
  98         mpfr_t mp_exactl, mp_xl, mp_yl;
  99         mpfr_t mp_exact, mp_x, mp_y;
 100         long double xl, yl, computedl, exactl, myulpl;
 101         double x, y, computed, exact, myulp;
 102         size_t i;
 103
 104         /* Initialize high precision variables */
 105         mpfr_inits2(200, mp_exact,  mp_x,  mp_y,  (mpfr_ptr)NULL);
 106         mpfr_inits2(200, mp_exactl, mp_xl, mp_yl, (mpfr_ptr)NULL);
 107
 108         ULP_INIT(u);
 109
 110         for (i = 0; i < NITERATIONS; i++) {
 111                 if (i % 100 == 0) {
 112                         printf("Percentage complete: %2.2f\r",
 113                             (100.0 * i)/NITERATIONS);
 114                         fflush(stdout);
 115                 }
 116
 117                 /* Generate random arguments */
 118                 populate_vars(f, &x, &y, &xl, &yl, mp_x, mp_y, mp_xl, mp_yl);
 119
 120                 /*
 121                  * Ready to call the mpfr*()
 122                  * The double version ALWAYS exist!
 123                  */
 124                 if(f->f_narg == 1) {
 125                         f->f_mpfr(mp_exact, mp_x, tonearest);
 126                 } else {
 127                         /* Hack, yikes */
 128                         if (!strcmp(f->f_name, "yn"))
 129                                 f->f_mpfr(mp_exact, (long)x, mp_y, tonearest);
 130                         else
 131                                 f->f_mpfr(mp_exact, mp_x, mp_y, tonearest);
 132                 }
 133                 exact = mpfr_get_d(mp_exact,  tonearest);
 134
 135                 /* We can't tell the same for long double functions though */
 136                 if (f->f_libml_real) {
 137                         if(f->f_narg == 1) {
 138                                 f->f_mpfr(mp_exactl, mp_xl, tonearest);
 139                         } else {
 140                                 f->f_mpfr(mp_exactl, mp_xl, mp_yl, tonearest);
 141                         }
 142                         exactl = mpfr_get_ld(mp_exactl, tonearest);
 143                 }
 144
 145                 /* Ready to call the libm*() */
 146                 if (f->f_narg == 1) {
 147                         computed = f->f_libm_real(x);
 148                 } else {
 149                         computed = f->f_libm_real(x, y);
 150                 }
 151
 152                 if (f->f_libml_real) {
 153                         if (f->f_narg == 1) {
 154                                 computedl = f->f_libml_real(xl);
 155                         } else {
 156                                 computedl = f->f_libml_real(xl, yl);
 157                         }
 158                 }
 159
 160                 /* Skip bogus results */
 161                 if (fpclassify(computed) == FP_NORMAL &&
 162                     fpclassify(exact) == FP_NORMAL) {
 163                         myulp = calculp_double(computed, exact);
 164                         ULP_UPDATE(u, myulp);
 165                 } else {
 166                         u->ulp_skipped++;
 167                 }
 168
 169                 if (f->f_libml_real) {
 170                         if (fpclassify(computedl) == FP_NORMAL &&
 171                             fpclassify(exactl) == FP_NORMAL) {
 172                                 myulpl = calculp_long_double(computedl, exactl);
 173                                 ULP_UPDATEL(u, myulpl);
 174                         }
 175                 } else {
 176                         u->ulp_skippedl++;
 177                 }
 178         }
 179
 180         u->ulp_avg  /= (i - u->ulp_skipped);
 181         u->ulp_avgl /= (i - u->ulp_skippedl);
 182
 183         /* Free resources */
 184         mpfr_clears(mp_exact,  mp_x,  mp_y,  (mpfr_ptr)NULL);
 185         mpfr_clears(mp_exactl, mp_xl, mp_yl, (mpfr_ptr)NULL);
 186 }
 187
 188 void
 189 printulps_double(struct ulp u)
 190 {
 191         if (u.ulp_max > 9.9 || u.ulp_min > 9.9) {
 192                 printf("%-10.4e %-10.e %-10.4e   ",
 193                     u.ulp_max, u.ulp_min, u.ulp_avg);
 194         } else {
 195                 printf("%-10.4f %-10.4f %-10.4f   ",
 196                     u.ulp_max, u.ulp_min, u.ulp_avg);
 197         }
 198         printf("%5u\n", u.ulp_skipped);
 199 }
 200
 201 void
 202 printulps_long_double(struct ulp u)
 203 {
 204         if (u.ulp_maxl > 9.9 || u.ulp_minl > 9.9) {
 205                 printf("%-10.4e %-10.4e %-10.4e   ",
 206                     (double)u.ulp_maxl,
 207                     (double)u.ulp_minl,
 208                     (double)u.ulp_avgl);
 209         } else {
 210                 printf("%-10.4f %-10.4f %-10.4f   ",
 211                     (double)u.ulp_maxl,
 212                     (double)u.ulp_minl,
 213                     (double)u.ulp_avgl);
 214         }
 215         printf("%5u\n", u.ulp_skippedl);
 216 }