ulps/ulp_complex.c

   1 #define _XOPEN_SOURCE 600
   2
   3 #include <assert.h>
   4 #include <complex.h>
   5 #include <float.h>
   6 #include <math.h>
   7 #include <stdio.h>
   8 #include <string.h>
   9
  10 #include <gmp.h>
  11 #include <mpc.h>
  12 #include <mpfr.h>
  13
  14 #include "gen.h"
  15 #include "mytypes.h"
  16 #include "subr_mpc.h"
  17 #include "subr_random.h"
  18 #include "ulp.h"
  19
  20 static double complex
  21 calculp_double_complex(double complex computed, double complex exact)
  22 {
  23         double_complex z;
  24         double ulp_real;
  25         double ulp_imag;
  26
  27         ulp_real = calculp_double(creal(computed), creal(exact));
  28         ulp_imag = calculp_double(cimag(computed), cimag(exact));
  29
  30         z.parts[0] = ulp_real;
  31         z.parts[1] = ulp_imag;
  32
  33         return (z.z);
  34 }
  35
  36 static long double complex
  37 calculp_long_double_complex(long double complex computedl,
  38                             long double complex exactl)
  39 {
  40         long_double_complex z;
  41         long double ulp_reall;
  42         long double ulp_imagl;
  43
  44         ulp_reall = calculp_long_double(creall(computedl), creall(exactl));
  45         ulp_imagl = calculp_long_double(cimagl(computedl), cimagl(exactl));
  46
  47         z.parts[0] = ulp_reall;
  48         z.parts[1] = ulp_imagl;
  49
  50         return (z.z);
  51 }
  52
  53 static void
  54 populate_complex_vars(const struct fentry *f,
  55     double complex *x, double complex *y,
  56     long double complex *xl, long double complex *yl,
  57     mpc_t mp_x, mpc_t mp_y, mpc_t mp_xl, mpc_t mp_yl)
  58 {
  59         const mpc_rnd_t tonearest = MPC_RNDNN;
  60         long double complex txl, tyl;
  61         double complex tx, ty;
  62
  63         assert(f && x && y && xl && yl);
  64         txl = tx = 0.0;
  65         tyl = ty = 0.0;
  66
  67         /* Generate random arguments */
  68         if (f->f_narg == 1) {
  69                 do {
  70                         tx  = random_double_complex(FP_NORMAL);
  71                         txl = random_long_double_complex(FP_NORMAL);
  72                 } while (!f->f_uc.fp1(tx) || !f->f_uc.fp1(txl));
  73         }
  74         if (f->f_narg == 2) {
  75                 do {
  76                         tx  = random_double_complex(FP_NORMAL);
  77                         ty  = random_double_complex(FP_NORMAL);
  78                         txl = random_long_double_complex(FP_NORMAL);
  79                         tyl = random_long_double_complex(FP_NORMAL);
  80                 } while (!f->f_uc.fp2(tx, ty) || !f->f_uc.fp2(txl, tyl));
  81         }
  82
  83         *x  = tx;
  84         *y  = ty;
  85         *xl = txl;
  86         *yl = tyl;
  87
  88         DPRINTF(( "x  = % .16e + I *% .16e\n",  creal (*x),  cimag (*x)));
  89         DPRINTF(( "y  = % .16e + I *% .16e\n",  creal (*y),  cimag (*y)));
  90         DPRINTF(("xl = % .16Le + I *% .16Le\n", creall(*xl), cimagl(*xl)));
  91         DPRINTF(("yl = % .16Le + I *% .16Le\n", creall(*yl), cimagl(*yl)));
  92
  93         /* Copy arguments to mpc variables */
  94         mpc_set_d_d  (mp_x,  creal (tx),  cimag (tx),  tonearest);
  95         mpc_set_d_d  (mp_y,  creal (ty),  cimag (ty),  tonearest);
  96         mpc_set_ld_ld(mp_xl, creall(txl), cimagl(txl), tonearest);
  97         mpc_set_ld_ld(mp_yl, creall(tyl), cimagl(tyl), tonearest);
  98 }
  99
 100 static int
 101 isbogus(long double complex z)
 102 {
 103         long double real = creall(z);
 104         long double imag = cimagl(z);
 105         int bogus = 0;
 106
 107         /* Real part */
 108         if ((fpclassify(real) != FP_NORMAL) &&
 109             (fpclassify(real) != FP_SUBNORMAL))
 110                 bogus = 1;
 111
 112         /* Imaginary part */
 113         if ((fpclassify(imag) != FP_NORMAL) &&
 114             (fpclassify(imag) != FP_SUBNORMAL))
 115                 bogus = 1;
 116
 117         return (bogus);
 118 }
 119
 120 void
 121 getfunctionulp_complex(const struct fentry *f, struct ulp_complex *uc)
 122 {
 123         const mpc_rnd_t tonearest = MPC_RNDNN;
 124         const mpfr_rnd_t mpfr_tonearest = MPFR_RNDN;
 125         mpc_t mp_exactl, mp_xl, mp_yl;
 126         mpc_t mp_exact, mp_x, mp_y;
 127         long double complex xl, yl, computedl, exactl, myulpl;
 128         double complex x, y, computed, exact, myulp;
 129         size_t i;
 130
 131         /* Initialize high precision variables */
 132         mpc_init2(mp_exact, 200);
 133         mpc_init2(mp_x,     200);
 134         mpc_init2(mp_y,     200);
 135         mpc_init2(mp_exactl,200);
 136         mpc_init2(mp_xl,    200);
 137         mpc_init2(mp_yl,    200);
 138
 139         ULP_COMPLEX_INIT(uc);
 140
 141         for (i = 0; i < NITERATIONS; i++) {
 142                 if (i % 100 == 0) {
 143                         printf("Percentage complete: %2.2f\r",
 144                             (100.0 * i)/NITERATIONS);
 145                         fflush(stdout);
 146                 }
 147
 148                 /* Generate random arguments */
 149                 DPRINTF(("populate_complex_vars: entering\n"));
 150                 populate_complex_vars(f, &x, &y, &xl, &yl, mp_x, mp_y, mp_xl, mp_yl);
 151                 DPRINTF(("populate_complex_vars: left\n"));
 152
 153                 /*
 154                  * Ready to call the mpc*()
 155                  * The double version ALWAYS exist!
 156                  */
 157                 char *str_x = mpc_get_str(10, 16, mp_x, tonearest);
 158                 char *str_y = mpc_get_str(10, 16, mp_y, tonearest);
 159                 DPRINTF(("f->f_mpc: enter x=%s, y=%s\n", str_x, str_y));
 160                 if(f->f_narg == 1) {
 161                         f->f_mpc(mp_exact, mp_x, tonearest);
 162                 } else {
 163                         f->f_mpc(mp_exact, mp_x, mp_y, tonearest);
 164                 }
 165                 exact = mpc_get_dc(mp_exact, mpfr_tonearest);
 166                 DPRINTF(("f->f_mpc: left\n"));
 167                 mpc_free_str(str_x);
 168                 mpc_free_str(str_y);
 169
 170                 /* We can't tell the same for long double functions though */
 171                 char *str_xl = mpc_get_str(10, 35, mp_xl, tonearest);
 172                 char *str_yl = mpc_get_str(10, 35, mp_yl, tonearest);
 173                 DPRINTF(("f->f_mpcl: enter x=%s, y=%s\n", str_xl, str_yl));
 174                 if (f->f_libml_complex) {
 175                         if(f->f_narg == 1) {
 176                                 f->f_mpc(mp_exactl, mp_xl, tonearest);
 177                         } else {
 178                                 f->f_mpc(mp_exactl, mp_xl, mp_yl, tonearest);
 179                         }
 180                         exactl = mpc_get_ldc(mp_exactl, mpfr_tonearest);
 181                 }
 182                 DPRINTF(("f->f_mpcl: left\n"));
 183                 mpc_free_str(str_xl);
 184                 mpc_free_str(str_yl);
 185
 186                 /* Ready to call the libm*() */
 187                 if (f->f_narg == 1) {
 188                         computed = f->f_libm_complex(x);
 189                 } else {
 190                         computed = f->f_libm_complex(x, y);
 191                 }
 192
 193                 if (f->f_libml_complex) {
 194                         if (f->f_narg == 1) {
 195                                 computedl = f->f_libml_complex(xl);
 196                         } else {
 197                                 computedl = f->f_libml_complex(xl, yl);
 198                         }
 199                 }
 200
 201                 /* Skip bogus results */
 202                 if (!isbogus(computed) && !isbogus(exact)) {
 203                         myulp = calculp_double_complex(computed, exact);
 204                         ULP_COMPLEX_UPDATE(uc, myulp);
 205                 } else {
 206                         uc->ulp_total.ulp_skipped++;
 207                 }
 208
 209                 if (f->f_libml_complex) {
 210                         if (!isbogus(computed) && !isbogus(exact)) {
 211                                 myulpl = calculp_long_double_complex(computedl, exactl);
 212                                 ULP_COMPLEX_UPDATEL(uc, myulpl);
 213                         }
 214                 } else {
 215                         uc->ulp_total.ulp_skippedl++;
 216                 }
 217         }
 218
 219         uc->ulp_total.ulp_avg  /= (i - uc->ulp_total.ulp_skipped);
 220
 221         /* Free resources */
 222         mpc_clear(mp_exact);
 223         mpc_clear(mp_x);
 224         mpc_clear(mp_y);
 225         mpc_clear(mp_exactl);
 226         mpc_clear(mp_xl);
 227         mpc_clear(mp_yl);
 228 }
 229
 230 void
 231 printulps_double_complex(struct ulp_complex uc)
 232 {
 233         if (uc.ulp_total.ulp_max > 9.9 || uc.ulp_total.ulp_min > 9.9) {
 234                 printf("%-10.4e %-10.4e %-10.4e   ",
 235                     uc.ulp_total.ulp_max,
 236                     uc.ulp_total.ulp_min,
 237                     uc.ulp_total.ulp_avg);
 238         } else {
 239                 printf("%-10.4f %-10.4f %-10.4f   ",
 240                     uc.ulp_total.ulp_max,
 241                     uc.ulp_total.ulp_min,
 242                     uc.ulp_total.ulp_avg);
 243         }
 244         printf("%5u\n", uc.ulp_total.ulp_skipped);
 245 }
 246
 247 void
 248 printulps_long_double_complex(struct ulp_complex uc)
 249 {
 250         if (uc.ulp_total.ulp_maxl > 9.9 || uc.ulp_total.ulp_minl > 9.9) {
 251                 printf("%-10.4e %-10.4e %-10.4e   ",
 252                     (double)uc.ulp_real.ulp_maxl,
 253                     (double)uc.ulp_real.ulp_minl,
 254                     (double)uc.ulp_real.ulp_avgl);
 255         } else {
 256         }
 257         printf("%5u\n", uc.ulp_total.ulp_skippedl);
 258 }