etc/profmath.c

   1 #include <assert.h>
   2 #include <math.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6 #include <time.h>
   7 #include <sys/time.h>
   8
   9 #include "gen.h"
  10 #include "subr_random.h"
  11
  12 #define NSAMPLE (1*1000)
  13
  14 /*
  15  * They are esentially the same, but we still define them as if
  16  * they were different, for conceptual consistency and clarity.
  17  */
  18 #define MARK_START(tv1)                                         \
  19         assert(gettimeofday(tv1, NULL) == 0);
  20
  21 #define MARK_END(tv2)                                           \
  22         assert(gettimeofday(tv2, NULL) == 0);
  23
  24 #define USECS(tv1, tv2)                                         \
  25         (((tv2).tv_sec  - (tv1).tv_sec ) * 1000000 +            \
  26          ((tv2).tv_usec - (tv1).tv_usec))
  27
  28 /*
  29  * Make these globals, so that we allocate memory only once
  30  * during the program's lifetime.
  31  */
  32 double *dx = NULL;
  33 double *dy = NULL;
  34 double *dz = NULL;
  35 double complex *cdx = NULL;
  36 double complex *cdy = NULL;
  37 double complex *cdz = NULL;
  38 struct timeval *start = NULL;
  39 struct timeval *end   = NULL;
  40
  41 static void
  42 init_arrays(void)
  43 {
  44         dx    = malloc(NSAMPLE * sizeof(double));
  45         dy    = malloc(NSAMPLE * sizeof(double));
  46         dz    = malloc(NSAMPLE * sizeof(double));
  47
  48         cdx    = malloc(NSAMPLE * sizeof(double complex));
  49         cdy    = malloc(NSAMPLE * sizeof(double complex));
  50         cdz    = malloc(NSAMPLE * sizeof(double complex));
  51
  52         start = malloc(NSAMPLE * sizeof(struct timeval));
  53         end   = malloc(NSAMPLE * sizeof(struct timeval));
  54
  55         assert(dx && dy && dz && start && end);
  56 }
  57
  58 static void
  59 reset(void)
  60 {
  61         assert(dx && dy && dz && start && end);
  62
  63         memset(dx,    0, NSAMPLE * sizeof(double));
  64         memset(dy,    0, NSAMPLE * sizeof(double));
  65         memset(dz,    0, NSAMPLE * sizeof(double));
  66
  67         memset(cdx,   0, NSAMPLE * sizeof(double complex));
  68         memset(cdy,   0, NSAMPLE * sizeof(double complex));
  69         memset(cdz,   0, NSAMPLE * sizeof(double complex));
  70
  71         memset(start, 0, NSAMPLE * sizeof(struct timeval));
  72         memset(end,   0, NSAMPLE * sizeof(struct timeval));
  73 }
  74
  75 static void
  76 cleanup(void)
  77 {
  78         assert(dx && dy && dz && start && end);
  79
  80         free(dx);  free(dy);  free(dz);
  81         free(cdx); free(cdy); free(cdz);
  82         free(start);
  83         free(end);
  84
  85         dx = dy = dz = NULL;
  86         cdx = cdy = cdz = NULL;
  87         start = end = NULL;
  88 }
  89
  90 static int
  91 proffunc(const char *fname)
  92 {
  93         const struct fentry *f;
  94         FILE *fp;
  95         size_t i, j;
  96
  97         /* Lookup the function */
  98         f = getfunctionbyname(fname);
  99         if (f == NULL)
 100                 return (-1);
 101
 102         /* Zero out the arrays we will be writing at */
 103         reset();
 104
 105         /* Generate random input -- do it before hand */
 106         for (i = 0; i < NSAMPLE; i++) {
 107                 if (f->f_mpfr) {
 108                         if (f->f_narg == 1) {
 109                                 do {
 110                                         dx[i] = random_double(FP_NORMAL);
 111                                 } while (!f->f_u.fp1(dx[i]));
 112                         }
 113                         if (f->f_narg >= 2) {
 114                                 do {
 115                                         dx[i] = random_double(FP_NORMAL);
 116                                         dy[i] = random_double(FP_NORMAL);
 117                                 } while (!f->f_u.fp2(dx[i], dy[i]));
 118                         }
 119                 } else {
 120                         if (f->f_narg == 1) {
 121                                 do {
 122                                         cdx[i] = random_double_complex(FP_NORMAL);
 123                                 } while (!f->f_uc.fp1(cdx[i]));
 124                         }
 125                         if (f->f_narg >= 2) {
 126                                 do {
 127                                         cdx[i] = random_double_complex(FP_NORMAL);
 128                                         cdy[i] = random_double_complex(FP_NORMAL);
 129                                 } while (!f->f_uc.fp2(cdx[i], cdy[i]));
 130                         }
 131                 }
 132         }
 133
 134         /* Ready to go */
 135         for (i = 0; i < NSAMPLE; i++) {
 136                 MARK_START(&start[i]);
 137
 138                 for (j = 0; j < 1000; j++) {
 139                         if (f->f_mpfr) {
 140                                 if (f->f_narg == 1) {
 141                                         dz[i] = f->f_libm_real(dx[i]);
 142                                 } else {
 143                                         dz[i] = f->f_libm_real(dx[i], dy[i]);
 144                                 }
 145                         } else {
 146                                 if (f->f_narg == 1) {
 147                                         cdz[i] = f->f_libm_complex(cdx[i]);
 148                                 } else {
 149                                         cdz[i] = f->f_libm_complex(cdx[i], cdy[i]);
 150                                 }
 151                         }
 152                 }
 153
 154                 MARK_END(&end[i]);
 155         }
 156
 157         /* Calculate diffs and dump them to the file */
 158         char buf[256];
 159         snprintf(buf, sizeof(buf), "%s.csv", fname);
 160         fp = fopen(buf, "w");
 161         assert(fp);
 162
 163         if (f->f_mpfr)
 164                 fprintf(fp, "#%d\n", f->f_narg);
 165         else
 166                 fprintf(fp, "#%d\n", f->f_narg + 1);
 167
 168         for (i = 0; i < NSAMPLE; i++) {
 169                 if (f->f_mpfr) {
 170                         if (fpclassify(dz[i]) != FP_NORMAL)
 171                                 continue;
 172
 173                         if (f->f_narg == 1) {
 174                                 fprintf(fp, "% .16e   % ld\n",
 175                                     dx[i], USECS(start[i], end[i]));
 176                         } else {
 177                                 fprintf(fp, "%.16e   %.16e   %ld\n",
 178                                     dx[i], dy[i], USECS(start[i], end[i]));
 179                         }
 180                 } else {
 181                         if ((fpclassify(creal(cdz[i])) != FP_NORMAL) ||
 182                              fpclassify(cimag(cdz[i])) != FP_NORMAL)
 183                                 continue;
 184
 185                         if (f->f_narg == 1) {
 186                                 fprintf(fp, "% .16e % .16e   % ld\n",
 187                                     creal(cdx[i]), cimag(cdx[i]),
 188                                     USECS(start[i], end[i]));
 189                         } else {
 190                                 fprintf(fp, "%.16e   %.16e   %.16e   %.16e %ld\n",
 191                                     creal(cdx[i]), cimag(cdx[i]),
 192                                     creal(cdy[i]), cimag(cdy[i]),
 193                                     USECS(start[i], end[i]));
 194                         }
 195                 }
 196         }
 197
 198         /* Success */
 199         return (0);
 200 }
 201
 202 int
 203 main(int argc, char *argv[])
 204 {
 205         int i, rv, total, all;
 206         const char *target;
 207         const struct fentry *f;
 208
 209         /* Skip program name */
 210         argv++;
 211         argc--;
 212
 213         /*
 214          * If none argument or "all" is given, then we assume
 215          * that all functions should be probed.
 216          */
 217         all = 0;
 218         if (argc == 0 || (argc == 1 && !strcmp(argv[0], "all"))) {
 219                 argv++;
 220                 argc--;
 221
 222                 all = 1;
 223         }
 224
 225         /* Initialize random number generator */
 226         init_randgen();
 227
 228         /* Allocate memory */
 229         init_arrays();
 230
 231         total = all ? fsize : argc;
 232
 233         for (i = 0; i < total; i++) {
 234                 f = getfunctionbyidx(i);
 235                 target = all ? f->f_name : argv[i];
 236                 printf("[%2d/%2d] %-10s : ", i+1, total, target);
 237                 fflush(stdout);
 238                 rv = proffunc(target);
 239                 if (rv != -1) {
 240                         printf("OK\n");
 241                 } else {
 242                         fprintf(stderr, "function: %s not probed -- skipping\n",
 243                             argv[i]);
 244                         continue;
 245                 }
 246         }
 247
 248         /* Free up resources */
 249         cleanup();
 250
 251         return (EXIT_SUCCESS);
 252 }