1#!/usr/bin/nickle
2
3# Use nickle's extended precision floating point implementation
4# to generate some simple test vectors for long double math functions
5
6typedef struct  {
7	real(real a) f;
8	string name;
9} func_f_f_t;
10
11typedef struct  {
12	real(real a, real b) f;
13	string name;
14} func_f_ff_t;
15
16typedef struct  {
17	real(real a, int b) f;
18	string name;
19} func_f_fi_t;
20
21typedef struct  {
22	int(real a) f;
23	string name;
24} func_i_f_t;
25
26string[] limited_funcs = {
27	"ceill",
28	"copysignl",
29	"fabsl",
30	"floorl",
31	"fmaxl",
32	"fminl",
33	"frexpl",
34	"hypotl",
35	"ilogbl",
36	"ldexpl",
37	"logbl",
38	"llrintl",
39	"lrintl",
40	"lroundl",
41	"llroundl",
42	"nanl",
43	"nearbyintl",
44	"rintl",
45	"roundl",
46	"scalbnl",
47	"scalblnl",
48	"truncl",
49	"logbl",
50	"sqrtl",
51};
52
53bool
54is_full_func(string name)
55{
56	for (int i = 0; i < dim(limited_funcs); i++)
57		if (limited_funcs[i] == name)
58			return false;
59	return true;
60}
61
62exception infinity(real v);
63exception nan();
64
65int prec = 512;
66int out_prec = 192;
67
68string
69toupper(string s)
70{
71	string o = "";
72	for (int i = 0; i < String::length(s); i++) {
73		int c = s[i];
74		if ('a' <= c && c <= 'z')
75			c = c - 'a' + 'A';
76		o = o + String::new(c);
77	}
78	return o;
79}
80
81string
82make_prec(string name)
83{
84	string prec = toupper(name) + "_PREC";
85	printf("#ifndef %s\n", prec);
86	printf("#define %s DEFAULT_PREC\n", prec);
87	printf("#endif\n");
88	return prec;
89}
90
91void
92gen_real_f_f(func_f_f_t f)
93{
94	real x, y;
95	string vec = sprintf("%s_vec", f.name);
96	printf("\n");
97	if (is_full_func(f.name))
98		printf("#ifdef FULL_LONG_DOUBLE\n");
99	string prec_name = make_prec(f.name);
100	printf("long_double_test_f_f_t %s[] = {\n", vec);
101	for (x = -10; x <= 10; x += .1) {
102		try {
103			string sy;
104			try {
105				try {
106					y = imprecise(f.f(imprecise(x, prec)), out_prec);
107					sy = sprintf("%.-eL", y);
108				} catch divide_by_zero(real x, real y) {
109					if (x == 0)
110						raise invalid_argument(f.name, 0, x);
111					raise infinity(x);
112				}
113			} catch infinity(real v) {
114				sy = "(long double) INFINITY";
115				if (v < 0)
116					sy = "-" + sy;
117			} catch nan() {
118				sy = "(long double) NAN";
119			}
120			printf("    { .line = __LINE__, .x = %.-eL, .y = %s },\n", x, sy);
121		} catch invalid_argument(string s, int i, poly x) {
122		}
123	}
124	printf("};\n\n");
125	printf("int test_%s(void) {\n", f.name);
126	printf("    unsigned int i;\n");
127	printf("    int result = 0;\n");
128	printf("    for (i = 0; i < sizeof(%s)/sizeof(%s[0]); i++) {\n", vec, vec);
129	printf("        long double y = %s(%s[i].x);\n", f.name, vec);
130	printf("        result += check_long_double(\"%s\", %s[i].line, %s, %s[i].y, y);\n", f.name, vec, prec_name, vec);
131	printf("    }\n");
132	printf("    return result;\n");
133	printf("}\n");
134	if (is_full_func(f.name))
135		printf("#endif /* FULL_LONG_DOUBLE */\n");
136}
137
138real cbrt(real x) { return x**(1/3); }
139real exp10(real x) { return 10**x; }
140real exp2(real x) { return 2**x; }
141real expm1(real x) {
142	x = imprecise(x);
143	int bits = precision(x);
144	int obits = bits;
145
146	if (0 < x && x < 1)
147		obits -= exponent(x);
148
149	x = imprecise(x, obits);
150
151	return imprecise(exp(x) - 1, bits);
152}
153real lgamma(real x) {
154	if (x < 0 && x == floor(x))
155		raise infinity(1);
156	return log(gamma(x));
157}
158real log1p(real x) { return log(1+x); }
159real logb(real x) {
160	if (x == 0)
161		raise infinity(-1);
162	return exponent(imprecise(x)) - 1;
163}
164real pow10(real x) { return 10**x; }
165
166real round(x) { if (x < 0) return -round(-x); return floor(x+0.5); }
167real trunc(x) { if (x < 0) return -trunc(-x); return floor(x); }
168
169real acosh(x) {
170	if (x < 1)
171		raise nan();
172	return log(x + sqrt(x*x-1));
173}
174
175real asinh(x) {
176	if (x == 0) return 0;
177	real sign = 1;
178	if (x < 0) {
179		sign = -1;
180		x = -x;
181	}
182	return sign * log(x + sqrt(x*x+1));
183}
184
185real atanh(x) {
186	if (abs(x) > 1)
187		raise nan();
188	if (abs(x) == 1)
189		raise infinity(x);
190	return 0.5 * log((1 + x) / (1 - x));
191}
192
193real cosh(x) {
194	return (exp(x) + exp(-x)) / 2;
195}
196
197real sinh(x) {
198	return (exp(x) - exp(-x)) / 2;
199}
200
201real tanh(x) {
202	return sinh(x) / cosh(x);
203}
204
205real tgamma(real x) {
206	if (x == 0)
207		raise infinity(1);
208	if (x < 0 && x == floor(x))
209		raise nan();
210	return gamma(x);
211}
212
213real nearbyint(real x) {
214	real y;
215
216	if (x < 0)
217		y = ceil(x-0.5);
218	else
219		y = floor(x+0.5);
220	if (abs(x-y) == 0.5) {
221		if (y % 2 != 0) {
222			if (y > 0)
223				y--;
224			else
225				y++;
226		}
227	}
228	return y;
229}
230
231real _erf(real x, real off)
232{
233	x = imprecise(x);
234	int bits = precision(x);
235	int obits = bits + 512;
236	real factor = 2 / sqrt(pi_value(obits));
237
238	x = imprecise(x, obits);
239	off = imprecise(off, obits) / factor;
240	real val = x - off;
241
242	for (int n = 1; ; n++) {
243		int f = 2 * n + 1;
244		real a = ((-1)**n * x**f) / (n! * f);
245		val += a;
246		if (exponent(val) - exponent(a) > obits)
247			break;
248	}
249	return imprecise(val * factor, bits);
250}
251
252real erf(real x)
253{
254	return _erf(x, 0);
255}
256
257real erfc(real x)
258{
259	return -_erf(x, 1);
260}
261
262real jn(real x, int n)
263{
264	x = imprecise(x);
265	int bits = precision(x);
266	int obits = bits + 512;
267
268	x = imprecise(x, obits);
269	real val = imprecise(0, obits);
270
271	for (int m = 0; ; m++) {
272		real a = ((-1)**m / (m! * gamma(m + n + 1))) * (x/2)**(2 * m + n);
273		val += a;
274		if (exponent(val) - exponent(a) > obits)
275			break;
276	}
277	return imprecise(val, bits);
278}
279
280real scalbnl(real x, int exp)
281{
282	return x * (2 ** exp);
283}
284
285real ldexpl(real x, int exp)
286{
287	return x * (2 ** exp);
288}
289
290real rintl(real x) {
291	return nearbyint(x);
292}
293
294real j0(real x) = jn(x,0);
295real j1(real x) = jn(x,1);
296
297real default_prec = 1e-20;
298
299func_f_f_t[] funcs_f_f = {
300	{ .f = acosh, .name = "acoshl" },
301	{ .f = acos, .name = "acosl" },
302	{ .f = asinh, .name = "asinhl" },
303	{ .f = asin, .name = "asinl" },
304	{ .f = atanh, .name = "atanhl" },
305	{ .f = atan, .name = "atanl" },
306	{ .f = cbrt, .name = "cbrtl" },
307	{ .f = ceil, .name = "ceill" },
308	{ .f = cosh, .name = "coshl" },
309	{ .f = cos, .name = "cosl" },
310	{ .f = erfc, .name = "erfcl" },
311	{ .f = erf, .name = "erfl" },
312	{ .f = exp10, .name = "exp10l" },
313	{ .f = exp2, .name = "exp2l" },
314	{ .f = exp, .name = "expl" },
315	{ .f = expm1, .name = "expm1l" },
316	{ .f = floor, .name = "floorl" },
317#	{ .f = j0, .name = "j0l" },
318#	{ .f = j1, .name = "j1l" },
319#	{ .f = jn, .name = "jnl" },
320	{ .f = lgamma, .name = "lgammal" },
321	{ .f = log10, .name = "log10l" },
322	{ .f = log1p, .name = "log1pl" },
323	{ .f = log2, .name = "log2l" },
324	{ .f = logb, .name = "logbl" },
325	{ .f = log, .name = "logl" },
326	{ .f = nearbyint, .name = "nearbyintl" },
327#	{ .f = pow10, .name = "pow10l" },	/* an alias for exp10 */
328        { .f = rintl, .name = "rintl" },
329	{ .f = round, .name = "roundl" },
330	{ .f = sinh, .name = "sinhl" },
331	{ .f = sin, .name = "sinl" },
332	{ .f = sqrt, .name = "sqrtl" },
333	{ .f = tanh, .name = "tanhl" },
334	{ .f = tan, .name = "tanl" },
335	{ .f = tgamma, .name = "tgammal" },
336	{ .f = trunc, .name = "truncl" },
337#	{ .f = y0, .name = "y0l" },
338#	{ .f = y1, .name = "y1l" },
339#	{ .f = yn, .name = "ynl" },
340};
341
342void
343gen_real_f_ff(func_f_ff_t f)
344{
345	real x0, x1, y;
346	string vec = sprintf("%s_vec", f.name);
347
348	printf("\n");
349	if (is_full_func(f.name))
350		printf("#ifdef FULL_LONG_DOUBLE\n");
351	string prec_name = make_prec(f.name);
352	printf("long_double_test_f_ff_t %s[] = {\n", vec);
353	for (x0 = -4; x0 <= 4; x0 += .25) {
354		for (x1 = -4; x1 <= 4; x1 += 0.25) {
355			try {
356				string sy;
357				try {
358					try {
359						y = imprecise(f.f(imprecise(x0, prec), imprecise(x1, prec)), out_prec);
360						sy = sprintf("%.-eL", y);
361					} catch divide_by_zero(real x, real y) {
362						if (x == 0)
363							raise invalid_argument(f.name, 0, x);
364						raise infinity(x);
365					}
366				} catch infinity(real v) {
367					sy = "(long double) INFINITY";
368					if (v < 0)
369						sy = "-" + sy;
370				} catch nan() {
371					sy = "(long double) NAN";
372				}
373				printf("    { .line = __LINE__, .x0 = %.-eL, .x1 = %.-eL, .y = %s },\n", x0, x1, sy);
374			} catch invalid_argument(string s, int i, poly x) {
375			}
376		}
377	}
378	printf("};\n\n");
379	printf("int test_%s(void) {\n", f.name);
380	printf("    unsigned int i;\n");
381	printf("    int result = 0;\n");
382	printf("    for (i = 0; i < sizeof(%s)/sizeof(%s[0]); i++) {\n", vec, vec);
383	printf("        long double y = %s(%s[i].x0, %s[i].x1);\n", f.name, vec, vec);
384	printf("        result += check_long_double(\"%s\", %s[i].line, %s, %s[i].y, y);\n", f.name, vec,prec_name, vec);
385	printf("    }\n");
386	printf("    return result;\n");
387	printf("}\n");
388	if (is_full_func(f.name))
389		printf("#endif /* FULL_LONG_DOUBLE */\n");
390}
391
392real fmod(real x, real y) {
393	if (y == 0)
394		raise nan();
395	real n = x / y;
396	if (n < 0)
397		n = ceil(n);
398	else
399		n = floor(n);
400	return x - n * y;
401}
402real fdim(real x, real y) { return max(x-y, 0); }
403real fmax(real x, real y) { return max(x,y); }
404real fmin(real x, real y) { return min(x,y); }
405
406real hypot(real x, real y) { return sqrt(x*x + y*y); }
407
408/* Compute an IEEE remainder */
409real remainder(real x, real y) {
410	if (y == 0)
411		raise nan();
412	real q = x / y;
413	int n;
414	if (q < 0)
415		n = ceil(q - 0.5);
416	else
417		n = floor(q + 0.5);
418	if (abs(q-n) == 0.5) {
419		if (n % 2 != 0) {
420			if (n > 0)
421				n--;
422			else
423				n++;
424		}
425	}
426	return x - n * y;
427}
428
429real drem(real x, real y) {
430	return remainder (x, y);
431}
432
433real copysign(real x, real y) {
434	x = abs(x);
435	if (y < 0)
436		x = -x;
437	return x;
438}
439
440bool
441isoddint(real x) {
442	return x == floor(x) && (floor(x) & 1) == 1;
443}
444
445bool
446isevenint(real x) {
447	return x == floor(x) && (floor(x) & 1) == 0;
448}
449
450bool
451isint(real x) {
452	return x == floor(x);
453}
454
455/* Deal with the oddities of IEEE pow */
456real powl(real x, real y) {
457	if (x == 0 && isoddint(y) && y < 0)
458		raise infinity(1);
459	if (x == 0 && y < 0)
460		raise infinity(1);
461	if (x == 0 && y > 0)
462		return 0;
463	if (x == 1)
464		return 1;
465	if (y == 0)
466		return 1;
467	if (x < 0 && !isint(y))
468		raise nan();
469	return pow(x, y);
470}
471
472real scalb(real x, real y) {
473	if (!isint(y))
474		raise nan();
475	return x * 2 ** y;
476}
477
478/* Functions of the form f(x,y) */
479func_f_ff_t[] funcs_f_ff = {
480        { .f = atan2, .name = "atan2l" },
481        { .f = powl, .name = "powl" },
482        { .f = fmod, .name = "fmodl" },
483#        { .f = nextafter, .name = "nextafterl" },
484#        { .f = nexttoward, .name = "nexttowardl" },
485        { .f = fdim, .name = "fdiml" },
486        { .f = fmax, .name = "fmaxl" },
487        { .f = fmin, .name = "fminl" },
488	{ .f = hypot, .name = "hypotl" },
489	{ .f = scalb, .name = "scalbl" },
490        { .f = remainder, .name = "remainderl" },
491        { .f = drem, .name = "dreml" },
492        { .f = copysign, .name = "copysignl" },
493};
494
495void
496gen_real_f_fi(func_f_fi_t f)
497{
498	real x0, y;
499	int x1;
500	string vec = sprintf("%s_vec", f.name);
501
502	printf("\n");
503	if (is_full_func(f.name))
504		printf("#ifdef FULL_LONG_DOUBLE\n");
505	string prec_name = make_prec(f.name);
506	printf("long_double_test_f_fi_t %s[] = {\n", vec);
507	for (x0 = -4; x0 <= 4; x0 += .25) {
508		for (x1 = -16; x1 <= 16; x1 += 1) {
509			try {
510				string sy;
511				try {
512					try {
513						y = imprecise(f.f(imprecise(x0, prec), x1), out_prec);
514						sy = sprintf("%.-eL", y);
515					} catch divide_by_zero(real x, real y) {
516						if (x == 0)
517							raise invalid_argument(f.name, 0, x);
518						raise infinity(x);
519					}
520				} catch infinity(real v) {
521					sy = "(long double) INFINITY";
522					if (v < 0)
523						sy = "-" + sy;
524				} catch nan() {
525					sy = "(long double) NAN";
526				}
527				printf("    { .line = __LINE__, .x0 = %.-eL, .x1 = %d, .y = %s },\n", x0, x1, sy);
528			} catch invalid_argument(string s, int i, poly x) {
529			}
530		}
531	}
532	printf("};\n\n");
533	printf("int test_%s(void) {\n", f.name);
534	printf("    unsigned int i;\n");
535	printf("    int result = 0;\n");
536	printf("    for (i = 0; i < sizeof(%s)/sizeof(%s[0]); i++) {\n", vec, vec);
537	printf("        long double y = %s(%s[i].x0, %s[i].x1);\n", f.name, vec, vec);
538	printf("        result += check_long_double(\"%s\", %s[i].line, %s, %s[i].y, y);\n", f.name, vec,prec_name, vec);
539	printf("    }\n");
540	printf("    return result;\n");
541	printf("}\n");
542	if (is_full_func(f.name))
543		printf("#endif /* FULL_LONG_DOUBLE */\n");
544}
545
546/* Functions of the form f(x,y) */
547func_f_fi_t[] funcs_f_fi = {
548	{ .f = ldexpl, .name = "ldexpl" },
549	{ .f = scalbnl, .name = "scalbnl" },
550};
551
552exception invalid_int(string y);
553
554void
555gen_real_i_f(func_i_f_t f)
556{
557	real x;
558	int y;
559	string vec = sprintf("%s_vec", f.name);
560
561	printf("\n");
562	if (is_full_func(f.name))
563		printf("#ifdef FULL_LONG_DOUBLE\n");
564	printf("long_double_test_i_f_t %s[] = {\n", vec);
565	for (x = -10; x <= 10; x += .1) {
566		try {
567			string sy;
568			try {
569				y = f.f(imprecise(x, prec));
570				sy = sprintf("%d", y);
571			} catch invalid_int(string s) {
572				sy = s;
573			}
574			printf("    { .line = __LINE__, .x = %.-eL, .y = %s },\n", x, sy);
575		} catch invalid_argument(string s, int i, poly x) {
576		}
577	}
578	printf("};\n\n");
579	printf("int test_%s(void) {\n", f.name);
580	printf("    unsigned int i;\n");
581	printf("    int result = 0;\n");
582	printf("    for (i = 0; i < sizeof(%s)/sizeof(%s[0]); i++) {\n", vec, vec);
583	printf("        long long y = %s(%s[i].x);\n", f.name, vec);
584	printf("        result += check_long_long(\"%s\", %s[i].line, %s[i].y, y);\n", f.name, vec, vec);
585	printf("    }\n");
586	printf("    return result;\n");
587	printf("}\n");
588	if (is_full_func(f.name))
589		printf("#endif /* FULL_LONG_DOUBLE */\n");
590}
591
592int finite(real x) {
593	return 1;
594}
595
596int ilogb(real x) {
597	if (x == 0)
598		raise invalid_int("FP_ILOGB0");
599	return exponent(imprecise(x)) - 1;
600}
601
602int isinf(real x) {
603	return 0;
604}
605
606int isnan(real x) {
607	return 0;
608}
609
610int lrint(real x) {
611	return rintl(x);
612}
613
614int lround(real x) {
615	int ix = floor(x);
616	real diff = x - ix;
617	if ((diff == 0.5) && (x > 0) || (diff > 0.5))
618		ix++;
619	return ix;
620}
621
622/* Functions of the form i(x) */
623func_i_f_t[] funcs_i_f = {
624	{ .f = finite, .name = "finitel" },
625	{ .f = ilogb, .name = "ilogb" },
626	{ .f = isinf, .name = "isinfl" },
627	{ .f = isnan, .name = "isnanl" },
628	{ .f = lrint, .name = "lrintl" },
629	{ .f = lrint, .name = "llrintl" },
630	{ .f = lround, .name = "lroundl" },
631	{ .f = lround, .name = "llroundl" },
632};
633
634
635/*
636 * These functions aren't tested yet
637 *
638 *  long double modfl (long double, long double *);
639 *  float nexttowardf (float, long double);
640 *  double nexttoward (double, long double);
641 *  long double nextowardl (long double, long double);
642 *  long double remquol (long double, long double, int *);
643 *  long double fmal (long double, long double, long double);
644 *  long double lgammal_r (long double, int *);
645 *  void sincosl (long double, long double *, long double *);
646 */
647
648void
649main()
650{
651	for (int i = 0; i < dim(funcs_i_f); i++)
652		gen_real_i_f(funcs_i_f[i]);
653
654	for (int i = 0; i < dim(funcs_f_fi); i++)
655		gen_real_f_fi(funcs_f_fi[i]);
656
657	for (int i = 0; i < dim(funcs_f_ff); i++)
658		gen_real_f_ff(funcs_f_ff[i]);
659
660	for (int i = 0; i < dim(funcs_f_f); i++)
661		gen_real_f_f(funcs_f_f[i]);
662
663	printf("long_double_test_t long_double_tests[] = {\n");
664	for (int i = 0; i < dim(funcs_f_f); i++) {
665		if (is_full_func(funcs_f_f[i].name))
666			printf("#ifdef FULL_LONG_DOUBLE\n");
667		printf("    { .name = \"%s\", .test = test_%s },\n", funcs_f_f[i].name, funcs_f_f[i].name);
668		if (is_full_func(funcs_f_f[i].name))
669			printf("#endif /* FULL_LONG_DOUBLE */\n");
670	}
671	for (int i = 0; i < dim(funcs_f_ff); i++) {
672		if (is_full_func(funcs_f_ff[i].name))
673			printf("#ifdef FULL_LONG_DOUBLE\n");
674		printf("    { .name = \"%s\", .test = test_%s },\n", funcs_f_ff[i].name, funcs_f_ff[i].name);
675		if (is_full_func(funcs_f_ff[i].name))
676			printf("#endif /* FULL_LONG_DOUBLE */\n");
677	}
678	for (int i = 0; i < dim(funcs_f_fi); i++) {
679		if (is_full_func(funcs_f_fi[i].name))
680			printf("#ifdef FULL_LONG_DOUBLE\n");
681		printf("    { .name = \"%s\", .test = test_%s },\n", funcs_f_fi[i].name, funcs_f_fi[i].name);
682		if (is_full_func(funcs_f_fi[i].name))
683			printf("#endif /* FULL_LONG_DOUBLE */\n");
684	}
685	for (int i = 0; i < dim(funcs_i_f); i++) {
686		if (is_full_func(funcs_i_f[i].name))
687			printf("#ifdef FULL_LONG_DOUBLE\n");
688		printf("    { .name = \"%s\", .test = test_%s },\n", funcs_i_f[i].name, funcs_i_f[i].name);
689		if (is_full_func(funcs_i_f[i].name))
690			printf("#endif /* FULL_LONG_DOUBLE */\n");
691	}
692	printf("};\n");
693}
694
695main();
696