ulps: Fix long lines
[mathlib.git] / 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 }