etc/realgen.c

   1 /*
   2  * The purpose of this utility is to generate valid C arrays of exact
   3  * (x, ..., func(x, ...)) pairs for consumption by the rest of the test cases.
   4  *
   5  * The code is a bit ugly. Sorry, tried to do it as less ugly as possible.
   6  */
   7 #define _XOPEN_SOURCE 600
   8
   9 #include <assert.h>
  10 #include <ctype.h>
  11 #include <errno.h>
  12 #include <math.h>
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include <unistd.h>
  17
  18 #include <gmp.h>
  19 #include <mpfr.h>
  20
  21 #include "gen.h"
  22 #include "subr_random.h"
  23
  24 /* Function prototypes */
  25 #define DECL_GEN(type)                                                  \
  26 static void gen_##type(const char *fname, size_t total,                 \
  27                         type lower, type upper, type delta)
  28 typedef long double long_double;
  29 DECL_GEN(double);
  30 DECL_GEN(long_double);
  31 static void usage(const char *progname);
  32 static char *strtoupper(const char *str);
  33 static long double mystrtold(const char *str, const char *what);
  34 static void ifdef_guard_open(const char *fname);
  35 static void ifdef_guard_close(const char *fname);
  36 static void struct_decl_open(const char *fname, const char *type);
  37 static void struct_decl_close(void);
  38
  39 /*
  40  * Bite the bullet and don't protect macro arguments with surrounding
  41  * parentheses nor use intermediate variables. Trade safety, for pleasure.
  42  */
  43 #define COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, rndmode)             \
  44 do {                                                                    \
  45         if (f->f_narg == 1)                                             \
  46                 f->f_mpfr(mp_exact, mp_x, rndmode);                     \
  47         else                                                            \
  48                 f->f_mpfr(mp_exact, mp_x, mp_y, rndmode);               \
  49 } while(0)
  50
  51 /*
  52  * `d' represents the symmetric interval around zero that we exclude.
  53  * We do so because floating-point numbers tend to become more dense
  54  * around zero and less so towards the infinities. And we don't want
  55  * to be massively spammed with tiny values (I've seen it happening).
  56  *
  57  * A little bit suboptimal for the f_narg = 2 case, as we may be
  58  * discarding -say- valid x values, because y was out of bounds.
  59  */
  60 #define GET_RANDOM_VAL(type, f, x, y, d)                                \
  61 do {                                                                    \
  62         if (f->f_narg == 1) {                                           \
  63                 do {                                                    \
  64                         x = random_##type(FP_NORMAL);                   \
  65                 } while (x < lower || x > upper || fabs(x) < d ||       \
  66                          !f->f_u.fp1(x));                               \
  67         }                                                               \
  68                                                                         \
  69         if (f->f_narg == 2) {                                           \
  70                 do {                                                    \
  71                         x = random_##type(FP_NORMAL);                   \
  72                         y = random_##type(FP_NORMAL);                   \
  73                 } while (x < lower || x > upper ||                      \
  74                          y < lower || y > upper ||                      \
  75                          fabsl(x) < d || fabsl(y) < d ||                \
  76                          !f->f_u.fp2(x, y));                            \
  77         }                                                               \
  78 } while(0)
  79
  80 /*
  81  * Skip the inexact exception
  82  *
  83  * In case of a non-zero real rounded result, the error on the result is less
  84  * or equal to 1/2 ulp (unit in the last place) of that result in the rounding
  85  * to nearest mode, and less than 1 ulp of that result in the directed rounding
  86  * modes. We use much more bits than double / long double anyway, so a 0.5 ulp
  87  * or even 1 ulp doesn't make any difference.
  88  */
  89 #define MPFR_ERROR(x)                                                   \
  90         (mpfr_underflow_p() || mpfr_overflow_p() || mpfr_nanflag_p()    \
  91             || mpfr_erangeflag_p())
  92
  93 /* Used for debugging purposes only */
  94 #define MPFR_PRINT_FLAGS(x)                                             \
  95 do {                                                                    \
  96         if (mpfr_underflow_p())                                         \
  97                 fprintf(stderr, "Underflow\n");                         \
  98         if (mpfr_overflow_p())                                          \
  99                 fprintf(stderr, "Overflow\n");                          \
 100         if (mpfr_nanflag_p())                                           \
 101                 fprintf(stderr, "NaN\n");                               \
 102         if (mpfr_inexflag_p())                                          \
 103                 fprintf(stderr, "Inexact\n");                           \
 104         if (mpfr_erangeflag_p())                                        \
 105                 fprintf(stderr, "Erange\n");                            \
 106         fprintf(stderr, "\n");                                          \
 107 } while(0)
 108
 109 #define MPFR_PRINT_FPCLASS(x)                                           \
 110 do {                                                                    \
 111         if (mpfr_nan_p(x))                                              \
 112                 fprintf(stderr, "Nan\n");                               \
 113         if (mpfr_inf_p(x))                                              \
 114                 fprintf(stderr, "Inf\n");                               \
 115         if (mpfr_number_p(x))                                           \
 116                 fprintf(stderr, "Number (neither NaN nor Inf)\n");      \
 117         if (mpfr_zero_p(x))                                             \
 118                 fprintf(stderr, "Zero\n");                              \
 119 } while(0)
 120
 121 int
 122 main(int argc, char *argv[])
 123 {
 124         const char *progname, *fname, *type;
 125         long double min, max, delta;
 126         long double ldlower, ldupper;
 127         double dlower, dupper;
 128         int total, opt;
 129
 130         /*
 131          * Save program name, function name and floating-point type.
 132          * We will be using them later.
 133          */
 134         progname = argv[0];
 135         fname = argv[1];
 136         type = argv[2];
 137         optind += 2;
 138
 139         /* Parse arguments */
 140         total = -1;
 141         while ((opt = getopt(argc, argv, "d:m:M:t:")) != -1) {
 142                 switch (opt) {
 143                 case 'd':
 144                         delta = mystrtold(optarg, "delta");
 145                         break;
 146                 case 'm':
 147                         min = mystrtold(optarg, "min");
 148                         break;
 149                 case 'M':
 150                         max = mystrtold(optarg, "Max");
 151                         break;
 152                 case 't':
 153                         total = atoi(optarg);
 154                         break;
 155                 default:
 156                         usage(progname);
 157                         /* NEVER REACHED */
 158                 }
 159         }
 160
 161         /* optind is the argv[] index of the first non-option element */
 162         if (optind < argc) {
 163                 fprintf(stderr, "non-option argv[]-elements: ");
 164                 while (optind < argc)
 165                         fprintf(stderr, "%s ", argv[optind++]);
 166                 fprintf(stderr, "\n");
 167                 usage(progname);
 168                 /* NEVER REACHED */
 169         }
 170
 171         /* Validate input -- make sure all arguments were given */
 172         if (total < 0 || min > max || argc != 7) {
 173                 usage(progname);
 174                 /* NEVER REACHED */
 175         }
 176
 177         if (strcmp(type, "float") && strcmp(type, "double")
 178             && (strcmp(type, "ldouble") && strcmp(type, "all"))) {
 179                 usage(progname);
 180                 /* NEVER REACHED */
 181         }
 182
 183         /* Initialize random number generator */
 184         init_randgen();
 185
 186         /* Open #ifdef guard */
 187         ifdef_guard_open(fname);
 188
 189         /* Ready to go */
 190         if (!strcmp(type, "double") || !strcmp(type, "all")) {
 191                 dlower = min;
 192                 dupper = max;
 193                 struct_decl_open(fname, "double");
 194                 gen_double(fname, total, dlower, dupper, delta);
 195                 struct_decl_close();
 196         }
 197
 198         if (!strcmp(type, "ldouble") || !strcmp(type, "all")) {
 199                 ldlower = min;
 200                 ldupper = max;
 201                 struct_decl_open(fname, "long double");
 202                 gen_long_double(fname, total, ldlower, ldupper, delta);
 203                 struct_decl_close();
 204         }
 205
 206         /* Close #ifdef guard */
 207         ifdef_guard_close(fname);
 208
 209         return (EXIT_SUCCESS);
 210 }
 211
 212 static void
 213 gen_double(const char *fname, size_t total,
 214     double lower, double upper, double delta)
 215 {
 216         const struct fentry *f;
 217         const mpfr_rnd_t tonearest = GMP_RNDN;
 218         mpfr_t mp_x, mp_y, mp_exact;
 219         double x, y, exact;
 220         size_t i;
 221
 222         /* Look up the function */
 223         f = getfunctionbyname(fname);
 224         assert(f);
 225
 226         /* Initialize high precision variables */
 227         mpfr_inits2(97, mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
 228
 229         for (i = 0; i < total; i++) {
 230                 /* Generate a random input for function f */
 231                 GET_RANDOM_VAL(double, f, x, y, delta);
 232
 233                 /* Set the mpfr variables */
 234                 if (f->f_narg >= 1)
 235                         mpfr_set_d(mp_x, x, tonearest);
 236                 if (f->f_narg == 2)
 237                         mpfr_set_d(mp_y, y, tonearest);
 238                 mpfr_set_d(mp_exact, 0.0, tonearest);
 239
 240                 /* Compute exact value, mp_exact = f(mp_x, ...)  */
 241                 mpfr_clear_flags();
 242                 COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, tonearest);
 243                 if (MPFR_ERROR() || !mpfr_number_p(mp_exact)) {
 244                         fprintf(stderr, "ERROR: exactf value overflowed\n");
 245                         i--;
 246                         continue;
 247                 }
 248
 249                 /* mpfr_fprintf(stderr, "-> %.35RNg\n", mp_exact); */
 250
 251                 /* Extract exact value */
 252                 exact = mpfr_get_d(mp_exact, tonearest);
 253                 if (fpclassify(exact) != FP_NORMAL) {
 254                         i--;
 255                         continue;
 256                 }
 257
 258                 if (f->f_narg == 1)
 259                         printf("\t{ % .16e, % .16e }", x, exact);
 260                 else
 261                         printf("\t{ % .16e,\n\t  % .16e,\n\t  % .16e }", x, y, exact);
 262
 263                 if (i < total - 1)
 264                         printf(",\n");
 265                 else
 266                         printf("\n");
 267         }
 268
 269         /* Free resources */
 270         mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
 271 }
 272
 273 static void
 274 gen_long_double(const char *fname, size_t total,
 275     long double lower, long double upper, long double delta)
 276 {
 277         const struct fentry *f;
 278         const mpfr_rnd_t tonearest = GMP_RNDN;
 279         mpfr_t mp_x, mp_y, mp_exact;
 280         long double x, y, exact;
 281         size_t i;
 282
 283         /* Lookup the function */
 284         f = getfunctionbyname(fname);
 285         assert(f);
 286
 287         /* Initialize high precision variables */
 288         mpfr_inits2(97, mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
 289
 290         for (i = 0; i < total; i++) {
 291                 /* Generate a random input for function f */
 292                 GET_RANDOM_VAL(long_double, f, x, y, delta);
 293
 294                 /* Set the mpfr variables */
 295                 if (f->f_narg >= 1)
 296                         mpfr_set_ld(mp_x, x, tonearest);
 297                 if (f->f_narg == 2)
 298                         mpfr_set_ld(mp_y, y, tonearest);
 299                 mpfr_set_ld(mp_exact, 0.0, tonearest);
 300
 301                 /* Compute exact value, mp_exact = f(mp_x, ...) */
 302                 mpfr_clear_flags();
 303                 COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, tonearest);
 304                 if (MPFR_ERROR() || !mpfr_number_p(mp_exact)) {
 305                         i--;
 306                         continue;
 307                 }
 308
 309                 /* mpfr_fprintf(stderr, "-> %.35RNg\n", mp_exact); */
 310
 311                 /* Extract exact value */
 312                 exact = mpfr_get_ld(mp_exact, tonearest);
 313                 if (fpclassify(exact) != FP_NORMAL) {
 314                         fprintf(stderr, "ERROR: exactl value overflowed\n");
 315                         i--;
 316                         continue;
 317                 }
 318
 319                 /*
 320                  * Don't forget the L suffix in long double floating-point
 321                  * constants, as by default they are treated as double.
 322                  */
 323                 if (f->f_narg == 1)
 324                         printf("\t{ % .35LeL, % .35LeL }", x, exact);
 325                 else
 326                         printf("\t{ % .35LeL, % .35LeL,\n\t% .35LeL }", x, y, exact);
 327
 328                 if (i < total - 1)
 329                         printf(",\n");
 330                 else
 331                         printf("\n");
 332         }
 333
 334         /* Free resources */
 335         mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
 336 }
 337
 338 static void
 339 usage(const char *progname)
 340 {
 341         fprintf(stderr,
 342             "usage: %s fname fptype -t total -m min -M Max -d delta\n"
 343             "\t`fname' is the function name, as declared in math.h\n"
 344             "\t`fptype' is one of `double', `ldouble' or `all'\n"
 345             "\t`total' is the total number of (x, ..., f(x, ...)) pairs\n"
 346             "\t`delta' is the interval around zero to exclude\n"
 347             "\t`min' and `Max' are long double constants\n"
 348             "\t(+-INF or INFINITY are acceptable input)\n"
 349             "Example:\n"
 350             "%s sin all -t10 -m-Inf -MInf -d1E-10\n",
 351             progname, progname);
 352
 353         exit(EXIT_FAILURE);
 354 }
 355
 356 /*
 357  * The caller must free up the memory
 358  */
 359 static char *
 360 strtoupper(const char *str)
 361 {
 362         const char *oldp;
 363         char *newp, *newstr;
 364
 365         assert(str);
 366
 367         newstr = malloc(strlen(str)+1);
 368         assert(newstr);
 369
 370         oldp = str;
 371         newp = newstr;
 372         do {
 373                 *newp = toupper(*oldp);
 374                 newp++;
 375         } while(*oldp++);
 376
 377         return (newstr);
 378 }
 379
 380 static long double
 381 mystrtold(const char *str, const char *what)
 382 {
 383         long double rv;
 384         char *endp;
 385
 386         errno = 0;
 387
 388         rv = strtold(str, &endp);
 389
 390         if (*endp)
 391                 fprintf(stderr,
 392                     "WARNING: Conversion stopped at %s (was given = %s)\n"
 393                     "WARNING: Continuing with %s = %.35Le\n",
 394                     endp, str, what, rv);
 395
 396         /* We only allow over/under-flows to happen (ERANGE) */
 397         if (errno) {
 398                 perror("strtold");
 399                 if (errno != ERANGE)
 400                         exit(EXIT_FAILURE);
 401         }
 402
 403         return (rv);
 404 }
 405
 406 static void
 407 ifdef_guard_open(const char *fname)
 408 {
 409         char *FNAME;
 410
 411         FNAME = strtoupper(fname);
 412         printf("#ifndef __T_%s_H__\n", FNAME);
 413         printf("#define __T_%s_H__\n", FNAME);
 414         printf("\n");
 415
 416         free(FNAME);
 417 }
 418
 419 static void
 420 ifdef_guard_close(const char *fname)
 421 {
 422         char *FNAME;
 423
 424         FNAME = strtoupper(fname);
 425         printf("#endif  /* ! __T_%s_H__ */\n", FNAME);
 426
 427         free(FNAME);
 428 }
 429
 430 static void
 431 struct_decl_open(const char *fname, const char *type)
 432 {
 433         const struct fentry *f;
 434         const char *str;
 435
 436         f = getfunctionbyname(fname);
 437         assert(f);
 438
 439         str = strcmp(type, "double") == 0 ? "" : "l";
 440
 441         if (f->f_narg == 1) {
 442                 printf("const struct\nt1%sdentry {\n"
 443                     "\t%s x;    /* Input */\n"
 444                     "\t%s y;    /* Output */\n"
 445                     "} t1%sdtable[] = {\n",
 446                     str, type, type, str);
 447         } else {
 448                 printf("const struct\nt1%sdentry {\n"
 449                     "\t%s x;    /* Input */\n"
 450                     "\t%s y;    /* Input */\n"
 451                     "\t%s z;    /* Output */\n"
 452                     "} t1%sdtable[] = {\n",
 453                     str, type, type, type, str);
 454         }
 455 }
 456
 457 static void
 458 struct_decl_close(void)
 459 {
 460         printf("};\n");
 461 }