[mathlib.git] / etc / realgen.c

/*
 * The purpose of this utility is to generate valid C arrays of exact
 * (x, ..., func(x, ...)) pairs for consumption by the rest of the test cases.
 *
 * The code is a bit ugly. Sorry, tried to do it as less ugly as possible.
 */
#define _XOPEN_SOURCE 600

#include <assert.h>
#include <ctype.h>
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include <gmp.h>
#include <mpfr.h>

#include "gen.h"
#include "subr_random.h"

/* Function prototypes */
#define DECL_GEN(type)							\
static void gen_##type(const char *fname, size_t total,			\
			type lower, type upper, type delta)
typedef long double long_double;
DECL_GEN(double);
DECL_GEN(long_double);
static void usage(const char *progname);
static char *strtoupper(const char *str);
static long double mystrtold(const char *str, const char *what);
static void ifdef_guard_open(const char *fname);
static void ifdef_guard_close(const char *fname);
static void struct_decl_open(const char *fname, const char *type);
static void struct_decl_close(void);

/*
 * Bite the bullet and don't protect macro arguments with surrounding
 * parentheses nor use intermediate variables. Trade safety, for pleasure.
 */
#define	COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, rndmode)		\
do {									\
	if (f->f_narg == 1)						\
		f->f_mpfr(mp_exact, mp_x, rndmode);			\
	else								\
		f->f_mpfr(mp_exact, mp_x, mp_y, rndmode);		\
} while(0)

/*
 * `d' represents the symmetric interval around zero that we exclude.
 * We do so because floating-point numbers tend to become more dense
 * around zero and less so towards the infinities. And we don't want
 * to be massively spammed with tiny values (I've seen it happening).
 *
 * A little bit suboptimal for the f_narg = 2 case, as we may be
 * discarding -say- valid x values, because y was out of bounds.
 */
#define GET_RANDOM_VAL(type, f, x, y, d)				\
do {									\
	if (f->f_narg == 1) {						\
		do {							\
			x = random_##type(FP_NORMAL);			\
		} while (x < lower || x > upper || fabs(x) < d ||	\
			 !f->f_u.fp1(x));				\
	}								\
									\
	if (f->f_narg == 2) {						\
		do {							\
			x = random_##type(FP_NORMAL);			\
			y = random_##type(FP_NORMAL);			\
		} while (x < lower || x > upper ||			\
			 y < lower || y > upper ||			\
			 fabsl(x) < d || fabsl(y) < d ||		\
			 !f->f_u.fp2(x, y));				\
	}								\
} while(0)

/*
 * Skip the inexact exception
 *
 * In case of a non-zero real rounded result, the error on the result is less
 * or equal to 1/2 ulp (unit in the last place) of that result in the rounding
 * to nearest mode, and less than 1 ulp of that result in the directed rounding
 * modes. We use much more bits than double / long double anyway, so a 0.5 ulp
 * or even 1 ulp doesn't make any difference.
 */
#define MPFR_ERROR(x)							\
	(mpfr_underflow_p() || mpfr_overflow_p() || mpfr_nanflag_p() 	\
	    || mpfr_erangeflag_p())

/* Used for debugging purposes only */
#define MPFR_PRINT_FLAGS(x)						\
do {									\
	if (mpfr_underflow_p())						\
		fprintf(stderr, "Underflow\n");				\
	if (mpfr_overflow_p())						\
		fprintf(stderr, "Overflow\n");				\
	if (mpfr_nanflag_p())						\
		fprintf(stderr, "NaN\n");				\
	if (mpfr_inexflag_p())						\
		fprintf(stderr, "Inexact\n");				\
	if (mpfr_erangeflag_p())					\
		fprintf(stderr, "Erange\n");				\
	fprintf(stderr, "\n");						\
} while(0)

#define MPFR_PRINT_FPCLASS(x)						\
do {									\
	if (mpfr_nan_p(x))						\
		fprintf(stderr, "Nan\n");				\
	if (mpfr_inf_p(x))						\
		fprintf(stderr, "Inf\n");				\
	if (mpfr_number_p(x))						\
		fprintf(stderr, "Number (neither NaN nor Inf)\n");	\
	if (mpfr_zero_p(x))						\
		fprintf(stderr, "Zero\n");				\
} while(0)

int
main(int argc, char *argv[])
{
	const char *progname, *fname, *type;
	long double min, max, delta;
	long double ldlower, ldupper;
	double dlower, dupper;
	int total, opt;

	/*
	 * Save program name, function name and floating-point type.
	 * We will be using them later.
	 */
	progname = argv[0];
	fname = argv[1];
	type = argv[2];
	optind += 2;

	/* Parse arguments */
	total = -1;
	while ((opt = getopt(argc, argv, "d:m:M:t:")) != -1) {
		switch (opt) {
		case 'd':
			delta = mystrtold(optarg, "delta");
			break;
		case 'm':
			min = mystrtold(optarg, "min");
			break;
		case 'M':
			max = mystrtold(optarg, "Max");
			break;
		case 't':
			total = atoi(optarg);
			break;
		default:
			usage(progname);
			/* NEVER REACHED */
		}
	}

	/* optind is the argv[] index of the first non-option element */
	if (optind < argc) {
		fprintf(stderr, "non-option argv[]-elements: ");
		while (optind < argc)
			fprintf(stderr, "%s ", argv[optind++]);
		fprintf(stderr, "\n");
		usage(progname);
		/* NEVER REACHED */
	}

	/* Validate input -- make sure all arguments were given */
	if (total < 0 || min > max || argc != 7) {
		usage(progname);
		/* NEVER REACHED */
	}

	if (strcmp(type, "float") && strcmp(type, "double")
	    && (strcmp(type, "ldouble") && strcmp(type, "all"))) {
		usage(progname);
		/* NEVER REACHED */
	}

	/* Initialize random number generator */
	init_randgen();

	/* Open #ifdef guard */
	ifdef_guard_open(fname);

	/* Ready to go */
	if (!strcmp(type, "double") || !strcmp(type, "all")) {
		dlower = min;
		dupper = max;
		struct_decl_open(fname, "double");
		gen_double(fname, total, dlower, dupper, delta);
		struct_decl_close();
	}

	if (!strcmp(type, "ldouble") || !strcmp(type, "all")) {
		ldlower = min;
		ldupper = max;
		struct_decl_open(fname, "long double");
		gen_long_double(fname, total, ldlower, ldupper, delta);
		struct_decl_close();
	}

	/* Close #ifdef guard */
	ifdef_guard_close(fname);

	return (EXIT_SUCCESS);
}

static void
gen_double(const char *fname, size_t total,
    double lower, double upper, double delta)
{
	const struct fentry *f;
	const mpfr_rnd_t tonearest = GMP_RNDN;
	mpfr_t mp_x, mp_y, mp_exact;
	double x, y, exact;
	size_t i;

	/* Look up the function */
	f = getfunctionbyname(fname);
	assert(f);

	/* Initialize high precision variables */
	mpfr_inits2(97, mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);

	for (i = 0; i < total; i++) {
		/* Generate a random input for function f */
		GET_RANDOM_VAL(double, f, x, y, delta);

		/* Set the mpfr variables */
		if (f->f_narg >= 1)
			mpfr_set_d(mp_x, x, tonearest);
		if (f->f_narg == 2)
			mpfr_set_d(mp_y, y, tonearest);
		mpfr_set_d(mp_exact, 0.0, tonearest);

		/* Compute exact value, mp_exact = f(mp_x, ...)  */
		mpfr_clear_flags();
		COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, tonearest);
		if (MPFR_ERROR() || !mpfr_number_p(mp_exact)) {
                        fprintf(stderr, "ERROR: exactf value overflowed\n");
			i--;
			continue;
		}

                /* mpfr_fprintf(stderr, "-> %.35RNg\n", mp_exact); */

		/* Extract exact value */
		exact = mpfr_get_d(mp_exact, tonearest);
		if (fpclassify(exact) != FP_NORMAL) {
			i--;
			continue;
		}

		if (f->f_narg == 1)
			printf("\t{ % .16e, % .16e }", x, exact);
		else
			printf("\t{ % .16e,\n\t  % .16e,\n\t  % .16e }", x, y, exact);

		if (i < total - 1)
			printf(",\n");
		else
			printf("\n");
	}

	/* Free resources */
	mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
}

static void
gen_long_double(const char *fname, size_t total,
    long double lower, long double upper, long double delta)
{
	const struct fentry *f;
	const mpfr_rnd_t tonearest = GMP_RNDN;
	mpfr_t mp_x, mp_y, mp_exact;
	long double x, y, exact;
	size_t i;

	/* Lookup the function */
	f = getfunctionbyname(fname);
	assert(f);

	/* Initialize high precision variables */
	mpfr_inits2(97, mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);

	for (i = 0; i < total; i++) {
		/* Generate a random input for function f */
		GET_RANDOM_VAL(long_double, f, x, y, delta);

		/* Set the mpfr variables */
		if (f->f_narg >= 1)
			mpfr_set_ld(mp_x, x, tonearest);
		if (f->f_narg == 2)
			mpfr_set_ld(mp_y, y, tonearest);
		mpfr_set_ld(mp_exact, 0.0, tonearest);

		/* Compute exact value, mp_exact = f(mp_x, ...) */
		mpfr_clear_flags();
		COMPUTE_EXACT_VAL(f, mp_exact, mp_x, mp_y, tonearest);
		if (MPFR_ERROR() || !mpfr_number_p(mp_exact)) {
			i--;
			continue;
		}

		/* mpfr_fprintf(stderr, "-> %.35RNg\n", mp_exact); */

		/* Extract exact value */
		exact = mpfr_get_ld(mp_exact, tonearest);
		if (fpclassify(exact) != FP_NORMAL) {
			fprintf(stderr, "ERROR: exactl value overflowed\n");
			i--;
			continue;
		}

		/*
		 * Don't forget the L suffix in long double floating-point
		 * constants, as by default they are treated as double.
		 */
		if (f->f_narg == 1)
			printf("\t{ % .35LeL, % .35LeL }", x, exact);
		else
			printf("\t{ % .35LeL, % .35LeL,\n\t% .35LeL }", x, y, exact);

		if (i < total - 1)
			printf(",\n");
		else
			printf("\n");
	}

	/* Free resources */
	mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
}

static void
usage(const char *progname)
{
	fprintf(stderr,
	    "usage: %s fname fptype -t total -m min -M Max -d delta\n"
	    "\t`fname' is the function name, as declared in math.h\n"
	    "\t`fptype' is one of `double', `ldouble' or `all'\n"
	    "\t`total' is the total number of (x, ..., f(x, ...)) pairs\n"
	    "\t`delta' is the interval around zero to exclude\n"
	    "\t`min' and `Max' are long double constants\n"
	    "\t(+-INF or INFINITY are acceptable input)\n"
	    "Example:\n"
	    "%s sin all -t10 -m-Inf -MInf -d1E-10\n",
	    progname, progname);

	exit(EXIT_FAILURE);
}

/*
 * The caller must free up the memory
 */
static char *
strtoupper(const char *str)
{
	const char *oldp;
	char *newp, *newstr;

	assert(str);

	newstr = malloc(strlen(str)+1);
	assert(newstr);

	oldp = str;
	newp = newstr;
	do {
		*newp = toupper(*oldp);
		newp++;
	} while(*oldp++);

	return (newstr);
}

static long double
mystrtold(const char *str, const char *what)
{
	long double rv;
	char *endp;

	errno = 0;

	rv = strtold(str, &endp);

	if (*endp)
		fprintf(stderr,
		    "WARNING: Conversion stopped at %s (was given = %s)\n"
		    "WARNING: Continuing with %s = %.35Le\n",
		    endp, str, what, rv);

	/* We only allow over/under-flows to happen (ERANGE) */
	if (errno) {
		perror("strtold");
		if (errno != ERANGE)
			exit(EXIT_FAILURE);
	}

	return (rv);
}

static void
ifdef_guard_open(const char *fname)
{
	char *FNAME;

	FNAME = strtoupper(fname);
        printf("#ifndef __T_%s_H__\n", FNAME);
	printf("#define __T_%s_H__\n", FNAME);
        printf("\n");

	free(FNAME);
}

static void
ifdef_guard_close(const char *fname)
{
	char *FNAME;

	FNAME = strtoupper(fname);
	printf("#endif  /* ! __T_%s_H__ */\n", FNAME);

	free(FNAME);
}

static void
struct_decl_open(const char *fname, const char *type)
{
	const struct fentry *f;
	const char *str;

	f = getfunctionbyname(fname);
	assert(f);

	str = strcmp(type, "double") == 0 ? "" : "l";

	if (f->f_narg == 1) {
		printf("const struct\nt1%sdentry {\n"
		    "\t%s x;    /* Input */\n"
		    "\t%s y;    /* Output */\n"
		    "} t1%sdtable[] = {\n",
		    str, type, type, str);
	} else {
		printf("const struct\nt1%sdentry {\n"
		    "\t%s x;    /* Input */\n"
                    "\t%s y;    /* Input */\n"
		    "\t%s z;    /* Output */\n"
                    "} t1%sdtable[] = {\n",
		    str, type, type, type, str);
	}
}

static void
struct_decl_close(void)
{
	printf("};\n");
}
Commit	Line	Data
bc9df8f8	1	/*
4bc214f1	2	* The purpose of this utility is to generate valid C arrays of exact
4bc214f1	3	* (x, ..., func(x, ...)) pairs for consumption by the rest of the test cases.
f8501da8 SK	4	*
f8501da8 SK	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
3f789a7a	8
6d871857	9	#include <assert.h>
7be71e3c	10	#include <ctype.h>
fdbf2829	11	#include <errno.h>
3f789a7a	12	#include <math.h>
3f789a7a	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"
e4b6fe8f SK	22	#include "subr_random.h"
3f789a7a	23
6d871857	24	/* Function prototypes */
820e96f5	25	#define DECL_GEN(type) \
820e96f5	26	static void gen_##type(const char *fname, size_t total, \
d7f40c9d	27	type lower, type upper, type delta)
b090cc4b	28	typedef long double long_double;
993de363	29	DECL_GEN(double);
b090cc4b	30	DECL_GEN(long_double);
6d871857	31	static void usage(const char *progname);
7be71e3c	32	static char strtoupper(const char str);
fdbf2829	33	static long double mystrtold(const char str, const char what);
040bdf30 SK	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);
7be71e3c	38
0f546176	39	/*
4bc214f1	40	* Bite the bullet and don't protect macro arguments with surrounding
4bc214f1	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) \
	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); \
0f546176	49	} while(0)
b090cc4b	50
b090cc4b	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	61	do { \
	62	if (f->f_narg == 1) { \
	63	do { \
	64	x = random_##type(FP_NORMAL); \
d7f40c9d	65	} while (x < lower \|\| x > upper \|\| fabs(x) < d \|\| \
d7f40c9d	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) \
	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"); \
d97dce9c	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
3f789a7a	121	int
6d871857	122	main(int argc, char *argv[])
3f789a7a	123	{
4bc214f1	124	const char progname, fname, *type;
d7f40c9d	125	long double min, max, delta;
d7f40c9d	126	long double ldlower, ldupper;
f40ce656	127	double dlower, dupper;
d7f40c9d	128	int total, opt;
6d871857	129
6d871857	130	/*
4bc214f1	131	* Save program name, function name and floating-point type.
4bc214f1	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;
6d871857	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;
6d871857	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);
4bc214f1	168	/* NEVER REACHED */
6d871857	169	}
6d871857	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;
d7f40c9d	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;
d7f40c9d	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	}
f40ce656	205
3ab38300	206	/* Close #ifdef guard */
040bdf30	207	ifdef_guard_close(fname);
7be71e3c	208
6d871857	209	return (EXIT_SUCCESS);
6d871857	210	}
3f789a7a	211
6d871857	212	static void
d7f40c9d	213	gen_double(const char *fname, size_t total,
d7f40c9d	214	double lower, double upper, double delta)
6d871857	215	{
6d871857	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
6d871857	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); */
5e599d5a SK	250
d97dce9c	251	/* Extract exact value */
d97dce9c	252	exact = mpfr_get_d(mp_exact, tonearest);
5e599d5a	253	if (fpclassify(exact) != FP_NORMAL) {
abc118ca SK	254	i--;
abc118ca SK	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	}
3f789a7a	268
f8501da8	269	/* Free resources */
b0341cf8	270	mpfr_clears(mp_x, mp_y, mp_exact, (mpfr_ptr)NULL);
6d871857	271	}
6d871857	272
f40ce656	273	static void
b090cc4b	274	gen_long_double(const char *fname, size_t total,
d7f40c9d	275	long double lower, long double upper, long double delta)
f40ce656	276	{
f40ce656	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);
f8501da8 SK	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 */
f8501da8 SK	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);
3e29be9d SK	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); */
5e599d5a SK	310
d97dce9c	311	/* Extract exact value */
d97dce9c	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
	338	static void
	339	usage(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"
993de363	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
6d871857	353	exit(EXIT_FAILURE);
3f789a7a	354	}
7be71e3c	355
	356	/*
	357	* The caller must free up the memory
	358	*/
	359	static char *
	360	strtoupper(const char *str)
	361	{
f19a43ad	362	const char *oldp;
f19a43ad	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
	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	}
040bdf30 SK	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;
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";
fe5fe624 SK	440
040bdf30	441	if (f->f_narg == 1) {
fe5fe624	442	printf("const struct\nt1%sdentry {\n"
040bdf30 SK	443	"\t%s x; /* Input */\n"
040bdf30 SK	444	"\t%s y; /* Output */\n"
fe5fe624 SK	445	"} t1%sdtable[] = {\n",
fe5fe624 SK	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",
fe5fe624 SK	453	str, type, type, type, str);
040bdf30 SK	454	}
	455	}
	456
	457	static void
	458	struct_decl_close(void)
	459	{
fe5fe624	460	printf("};\n");
040bdf30	461	}