ulps: Fix long lines
[mathlib.git] / etc / realgen.c
CommitLineData
bc9df8f8 1/*
4bc214f1 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.
f8501da8
SK
4 *
5 * The code is a bit ugly. Sorry, tried to do it as less ugly as possible.
bc9df8f8 6 */
3f789a7a 7#define _XOPEN_SOURCE 600
8
6d871857 9#include <assert.h>
7be71e3c 10#include <ctype.h>
fdbf2829 11#include <errno.h>
3f789a7a 12#include <math.h>
13#include <stdio.h>
6d871857 14#include <stdlib.h>
15#include <string.h>
16#include <unistd.h>
3f789a7a 17
18#include <gmp.h>
19#include <mpfr.h>
20
e4b6fe8f
SK
21#include "gen.h"
22#include "subr_random.h"
3f789a7a 23
6d871857 24/* Function prototypes */
820e96f5 25#define DECL_GEN(type) \
26static void gen_##type(const char *fname, size_t total, \
d7f40c9d 27 type lower, type upper, type delta)
b090cc4b 28typedef long double long_double;
993de363 29DECL_GEN(double);
b090cc4b 30DECL_GEN(long_double);
6d871857 31static void usage(const char *progname);
7be71e3c 32static char *strtoupper(const char *str);
fdbf2829 33static long double mystrtold(const char *str, const char *what);
040bdf30
SK
34static void ifdef_guard_open(const char *fname);
35static void ifdef_guard_close(const char *fname);
36static void struct_decl_open(const char *fname, const char *type);
37static void struct_decl_close(void);
7be71e3c 38
0f546176 39/*
4bc214f1 40 * Bite the bullet and don't protect macro arguments with surrounding
41 * parentheses nor use intermediate variables. Trade safety, for pleasure.
0f546176 42 */
4bc214f1 43#define COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, rndmode) \
44do { \
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); \
0f546176 49} while(0)
b090cc4b 50
51/*
d7f40c9d 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 *
b090cc4b 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 */
d7f40c9d 60#define GET_RANDOM_VAL(type, f, x, y, d) \
b090cc4b 61do { \
62 if (f->f_narg == 1) { \
63 do { \
64 x = random_##type(FP_NORMAL); \
d7f40c9d 65 } while (x < lower || x > upper || fabs(x) < d || \
66 !f->f_u.fp1(x)); \
b090cc4b 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 || \
d7f40c9d 74 y < lower || y > upper || \
75 fabsl(x) < d || fabsl(y) < d || \
76 !f->f_u.fp2(x, y)); \
b090cc4b 77 } \
78} while(0)
79
d97dce9c 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 */
820e96f5 94#define MPFR_PRINT_FLAGS(x) \
95do { \
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"); \
d97dce9c 107} while(0)
108
109#define MPFR_PRINT_FPCLASS(x) \
110do { \
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
3f789a7a 121int
6d871857 122main(int argc, char *argv[])
3f789a7a 123{
4bc214f1 124 const char *progname, *fname, *type;
d7f40c9d 125 long double min, max, delta;
126 long double ldlower, ldupper;
f40ce656 127 double dlower, dupper;
d7f40c9d 128 int total, opt;
6d871857 129
130 /*
4bc214f1 131 * Save program name, function name and floating-point type.
132 * We will be using them later.
6d871857 133 */
134 progname = argv[0];
135 fname = argv[1];
136 type = argv[2];
6d871857 137 optind += 2;
138
4bc214f1 139 /* Parse arguments */
fdbf2829 140 total = -1;
d7f40c9d 141 while ((opt = getopt(argc, argv, "d:m:M:t:")) != -1) {
6d871857 142 switch (opt) {
d7f40c9d 143 case 'd':
144 delta = mystrtold(optarg, "delta");
145 break;
6d871857 146 case 'm':
fdbf2829 147 min = mystrtold(optarg, "min");
6d871857 148 break;
149 case 'M':
fdbf2829 150 max = mystrtold(optarg, "Max");
6d871857 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");
4bc214f1 167 usage(progname);
168 /* NEVER REACHED */
6d871857 169 }
170
fdbf2829 171 /* Validate input -- make sure all arguments were given */
d7f40c9d 172 if (total < 0 || min > max || argc != 7) {
f8501da8
SK
173 usage(progname);
174 /* NEVER REACHED */
175 }
6d871857 176
f8501da8 177 if (strcmp(type, "float") && strcmp(type, "double")
7be71e3c 178 && (strcmp(type, "ldouble") && strcmp(type, "all"))) {
f8501da8
SK
179 usage(progname);
180 /* NEVER REACHED */
181 }
3f789a7a 182
183 /* Initialize random number generator */
184 init_randgen();
185
3ab38300 186 /* Open #ifdef guard */
040bdf30 187 ifdef_guard_open(fname);
7be71e3c 188
6d871857 189 /* Ready to go */
7be71e3c 190 if (!strcmp(type, "double") || !strcmp(type, "all")) {
d7f40c9d 191 dlower = min;
192 dupper = max;
77aba2f5 193 struct_decl_open(fname, "double");
d7f40c9d 194 gen_double(fname, total, dlower, dupper, delta);
040bdf30 195 struct_decl_close();
f8501da8 196 }
6d871857 197
7be71e3c 198 if (!strcmp(type, "ldouble") || !strcmp(type, "all")) {
d7f40c9d 199 ldlower = min;
200 ldupper = max;
77aba2f5 201 struct_decl_open(fname, "long double");
d7f40c9d 202 gen_long_double(fname, total, ldlower, ldupper, delta);
040bdf30 203 struct_decl_close();
f40ce656 204 }
205
3ab38300 206 /* Close #ifdef guard */
040bdf30 207 ifdef_guard_close(fname);
7be71e3c 208
6d871857 209 return (EXIT_SUCCESS);
210}
3f789a7a 211
6d871857 212static void
d7f40c9d 213gen_double(const char *fname, size_t total,
214 double lower, double upper, double delta)
6d871857 215{
216 const struct fentry *f;
f8501da8 217 const mpfr_rnd_t tonearest = GMP_RNDN;
b090cc4b 218 mpfr_t mp_x, mp_y, mp_exact;
219 double x, y, exact;
220 size_t i;
6d871857 221
993de363 222 /* Look up the function */
f8501da8 223 f = getfunctionbyname(fname);
6d871857 224 assert(f);
6d871857 225
f8501da8 226 /* Initialize high precision variables */
b0341cf8 227 mpfr_inits2(97, mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
6d871857 228
229 for (i = 0; i < total; i++) {
4bc214f1 230 /* Generate a random input for function f */
d7f40c9d 231 GET_RANDOM_VAL(double, f, x, y, delta);
3f789a7a 232
6d871857 233 /* Set the mpfr variables */
f8501da8
SK
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);
3f789a7a 239
4bc214f1 240 /* Compute exact value, mp_exact = f(mp_x, ...) */
d97dce9c 241 mpfr_clear_flags();
993de363 242 COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, tonearest);
d97dce9c 243 if (MPFR_ERROR() || !mpfr_number_p(mp_exact)) {
abc118ca 244 fprintf(stderr, "ERROR: exactf value overflowed\n");
993de363 245 i--;
246 continue;
247 }
248
5e599d5a
SK
249 /* mpfr_fprintf(stderr, "-> %.35RNg\n", mp_exact); */
250
d97dce9c 251 /* Extract exact value */
252 exact = mpfr_get_d(mp_exact, tonearest);
5e599d5a 253 if (fpclassify(exact) != FP_NORMAL) {
abc118ca
SK
254 i--;
255 continue;
2a5d8617 256 }
d97dce9c 257
f8501da8 258 if (f->f_narg == 1)
7be71e3c 259 printf("\t{ % .16e, % .16e }", x, exact);
f8501da8
SK
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");
3f789a7a 267 }
268
f8501da8 269 /* Free resources */
b0341cf8 270 mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
6d871857 271}
272
f40ce656 273static void
b090cc4b 274gen_long_double(const char *fname, size_t total,
d7f40c9d 275 long double lower, long double upper, long double delta)
f40ce656 276{
277 const struct fentry *f;
c866c81e 278 const mpfr_rnd_t tonearest = GMP_RNDN;
b090cc4b 279 mpfr_t mp_x, mp_y, mp_exact;
280 long double x, y, exact;
281 size_t i;
f40ce656 282
993de363 283 /* Lookup the function */
f8501da8
SK
284 f = getfunctionbyname(fname);
285 assert(f);
f40ce656 286
f8501da8 287 /* Initialize high precision variables */
b0341cf8 288 mpfr_inits2(97, mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
f40ce656 289
f8501da8 290 for (i = 0; i < total; i++) {
4bc214f1 291 /* Generate a random input for function f */
d7f40c9d 292 GET_RANDOM_VAL(long_double, f, x, y, delta);
f40ce656 293
f8501da8
SK
294 /* Set the mpfr variables */
295 if (f->f_narg >= 1)
3e29be9d 296 mpfr_set_ld(mp_x, x, tonearest);
f8501da8 297 if (f->f_narg == 2)
3e29be9d
SK
298 mpfr_set_ld(mp_y, y, tonearest);
299 mpfr_set_ld(mp_exact, 0.0, tonearest);
f40ce656 300
4bc214f1 301 /* Compute exact value, mp_exact = f(mp_x, ...) */
d97dce9c 302 mpfr_clear_flags();
993de363 303 COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, tonearest);
d97dce9c 304 if (MPFR_ERROR() || !mpfr_number_p(mp_exact)) {
993de363 305 i--;
306 continue;
307 }
308
5e599d5a
SK
309 /* mpfr_fprintf(stderr, "-> %.35RNg\n", mp_exact); */
310
d97dce9c 311 /* Extract exact value */
312 exact = mpfr_get_ld(mp_exact, tonearest);
5e599d5a 313 if (fpclassify(exact) != FP_NORMAL) {
abc118ca
SK
314 fprintf(stderr, "ERROR: exactl value overflowed\n");
315 i--;
316 continue;
2a5d8617 317 }
d97dce9c 318
993de363 319 /*
320 * Don't forget the L suffix in long double floating-point
321 * constants, as by default they are treated as double.
322 */
f8501da8 323 if (f->f_narg == 1)
4eda6c9f 324 printf("\t{ % .35LeL, % .35LeL }", x, exact);
f8501da8 325 else
4eda6c9f 326 printf("\t{ % .35LeL, % .35LeL,\n\t% .35LeL }", x, y, exact);
f40ce656 327
f8501da8
SK
328 if (i < total - 1)
329 printf(",\n");
330 else
331 printf("\n");
332 }
f40ce656 333
f8501da8 334 /* Free resources */
b0341cf8 335 mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
f40ce656 336}
6d871857 337
338static void
339usage(const char *progname)
340{
341 fprintf(stderr,
d7f40c9d 342 "usage: %s fname fptype -t total -m min -M Max -d delta\n"
993de363 343 "\t`fname' is the function name, as declared in math.h\n"
344 "\t`fptype' is one of `double', `ldouble' or `all'\n"
d7f40c9d 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"
6fd49e2c
SK
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);
6d871857 352
353 exit(EXIT_FAILURE);
3f789a7a 354}
7be71e3c 355
356/*
357 * The caller must free up the memory
358 */
359static char *
360strtoupper(const char *str)
361{
f19a43ad 362 const char *oldp;
363 char *newp, *newstr;
7be71e3c 364
365 assert(str);
366
11876926 367 newstr = malloc(strlen(str)+1);
7be71e3c 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}
fdbf2829 379
380static long double
381mystrtold(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}
040bdf30
SK
405
406static void
407ifdef_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
419static void
420ifdef_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
430static void
431struct_decl_open(const char *fname, const char *type)
432{
433 const struct fentry *f;
fe5fe624 434 const char *str;
040bdf30
SK
435
436 f = getfunctionbyname(fname);
437 assert(f);
438
fe5fe624
SK
439 str = strcmp(type, "double") == 0 ? "" : "l";
440
040bdf30 441 if (f->f_narg == 1) {
fe5fe624 442 printf("const struct\nt1%sdentry {\n"
040bdf30
SK
443 "\t%s x; /* Input */\n"
444 "\t%s y; /* Output */\n"
fe5fe624
SK
445 "} t1%sdtable[] = {\n",
446 str, type, type, str);
040bdf30 447 } else {
fe5fe624 448 printf("const struct\nt1%sdentry {\n"
040bdf30
SK
449 "\t%s x; /* Input */\n"
450 "\t%s y; /* Input */\n"
451 "\t%s z; /* Output */\n"
fe5fe624
SK
452 "} t1%sdtable[] = {\n",
453 str, type, type, type, str);
040bdf30
SK
454 }
455}
456
457static void
458struct_decl_close(void)
459{
fe5fe624 460 printf("};\n");
040bdf30 461}