1/*-
2 * Copyright (c) 2012 David Schultz <das@FreeBSD.org>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 *    notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 *    notice, this list of conditions and the following disclaimer in the
12 *    documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 */
26
27/*
28 * Test that floating-point arithmetic works as specified by the C standard.
29 */
30
31#include <sys/cdefs.h>
32__FBSDID("$FreeBSD: releng/11.0/tools/regression/usr.bin/cc/float.c 230368 2012-01-20 06:57:21Z das $");
33
34#include <fenv.h>
35#include <float.h>
36#include <math.h>
37#include <stdio.h>
38
39#ifdef  __i386__
40#include <ieeefp.h>
41#endif
42
43#define	ALL_STD_EXCEPT	(FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \
44			 FE_OVERFLOW | FE_UNDERFLOW)
45
46#define	TWICE(x)		((x) + (x))
47#define	test(desc, pass)	test1((desc), (pass), 0)
48#define	skiptest(desc, pass)	test1((desc), (pass), 1)
49
50#pragma STDC FENV_ACCESS ON
51
52static const float one_f = 1.0 + FLT_EPSILON / 2;
53static const double one_d = 1.0 + DBL_EPSILON / 2;
54static const long double one_ld = 1.0L + LDBL_EPSILON / 2;
55
56static int testnum, failures;
57
58static void
59test1(const char *testdesc, int pass, int skip)
60{
61
62	testnum++;
63	printf("%sok %d - %s%s\n", pass || skip ? "" : "not ", testnum,
64	    skip ? "(SKIPPED) " : "", testdesc);
65	if (!pass && !skip)
66		failures++;
67}
68
69/*
70 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
71 */
72static int
73fpequal(long double d1, long double d2)
74{
75
76	if (d1 != d2)
77		return (isnan(d1) && isnan(d2));
78	return (copysignl(1.0, d1) == copysignl(1.0, d2));
79}
80
81void
82run_zero_opt_test(double d1, double d2)
83{
84
85	test("optimizations don't break the sign of 0",
86	     fpequal(d1 - d2, 0.0)
87	     && fpequal(-d1 + 0.0, 0.0)
88	     && fpequal(-d1 - d2, -0.0)
89	     && fpequal(-(d1 - d2), -0.0)
90	     && fpequal(-d1 - (-d2), 0.0));
91}
92
93void
94run_inf_opt_test(double d)
95{
96
97	test("optimizations don't break infinities",
98	     fpequal(d / d, NAN) && fpequal(0.0 * d, NAN));
99}
100
101static inline double
102todouble(long double ld)
103{
104
105	return (ld);
106}
107
108static inline float
109tofloat(double d)
110{
111
112	return (d);
113}
114
115void
116run_tests(void)
117{
118	volatile long double vld;
119	long double ld;
120	volatile double vd;
121	double d;
122	volatile float vf;
123	float f;
124	int x;
125
126	test("sign bits", fpequal(-0.0, -0.0) && !fpequal(0.0, -0.0));
127
128	vd = NAN;
129	test("NaN equality", fpequal(NAN, NAN) && NAN != NAN && vd != vd);
130
131	feclearexcept(ALL_STD_EXCEPT);
132	test("NaN comparison returns false", !(vd <= vd));
133	/*
134	 * XXX disabled; gcc/amd64 botches this IEEE 754 requirement by
135	 * emitting ucomisd instead of comisd.
136	 */
137	skiptest("FENV_ACCESS: NaN comparison raises invalid exception",
138	    fetestexcept(ALL_STD_EXCEPT) == FE_INVALID);
139
140	vd = 0.0;
141	run_zero_opt_test(vd, vd);
142
143	vd = INFINITY;
144	run_inf_opt_test(vd);
145
146	feclearexcept(ALL_STD_EXCEPT);
147	vd = INFINITY;
148	x = (int)vd;
149	/* XXX disabled (works with -O0); gcc doesn't support FENV_ACCESS */
150	skiptest("FENV_ACCESS: Inf->int conversion raises invalid exception",
151	    fetestexcept(ALL_STD_EXCEPT) == FE_INVALID);
152
153	/* Raising an inexact exception here is an IEEE-854 requirement. */
154	feclearexcept(ALL_STD_EXCEPT);
155	vd = 0.75;
156	x = (int)vd;
157	test("0.75->int conversion rounds toward 0, raises inexact exception",
158	     x == 0 && fetestexcept(ALL_STD_EXCEPT) == FE_INEXACT);
159
160	feclearexcept(ALL_STD_EXCEPT);
161	vd = -42.0;
162	x = (int)vd;
163	test("-42.0->int conversion is exact, raises no exception",
164	     x == -42 && fetestexcept(ALL_STD_EXCEPT) == 0);
165
166	feclearexcept(ALL_STD_EXCEPT);
167	x = (int)INFINITY;
168	/* XXX disabled; gcc doesn't support FENV_ACCESS */
169	skiptest("FENV_ACCESS: const Inf->int conversion raises invalid",
170	    fetestexcept(ALL_STD_EXCEPT) == FE_INVALID);
171
172	feclearexcept(ALL_STD_EXCEPT);
173	x = (int)0.5;
174	/* XXX disabled; gcc doesn't support FENV_ACCESS */
175	skiptest("FENV_ACCESS: const double->int conversion raises inexact",
176	     x == 0 && fetestexcept(ALL_STD_EXCEPT) == FE_INEXACT);
177
178	test("compile-time constants don't have too much precision",
179	     one_f == 1.0L && one_d == 1.0L && one_ld == 1.0L);
180
181	test("const minimum rounding precision",
182	     1.0F + FLT_EPSILON != 1.0F &&
183	     1.0 + DBL_EPSILON != 1.0 &&
184	     1.0L + LDBL_EPSILON != 1.0L);
185
186	/* It isn't the compiler's fault if this fails on FreeBSD/i386. */
187	vf = FLT_EPSILON;
188	vd = DBL_EPSILON;
189	vld = LDBL_EPSILON;
190	test("runtime minimum rounding precision",
191	     1.0F + vf != 1.0F && 1.0 + vd != 1.0 && 1.0L + vld != 1.0L);
192
193	test("explicit float to float conversion discards extra precision",
194	     (float)(1.0F + FLT_EPSILON * 0.5F) == 1.0F &&
195	     (float)(1.0F + vf * 0.5F) == 1.0F);
196	test("explicit double to float conversion discards extra precision",
197	     (float)(1.0 + FLT_EPSILON * 0.5) == 1.0F &&
198	     (float)(1.0 + vf * 0.5) == 1.0F);
199	test("explicit ldouble to float conversion discards extra precision",
200	     (float)(1.0L + FLT_EPSILON * 0.5L) == 1.0F &&
201	     (float)(1.0L + vf * 0.5L) == 1.0F);
202
203	test("explicit double to double conversion discards extra precision",
204	     (double)(1.0 + DBL_EPSILON * 0.5) == 1.0 &&
205	     (double)(1.0 + vd * 0.5) == 1.0);
206	test("explicit ldouble to double conversion discards extra precision",
207	     (double)(1.0L + DBL_EPSILON * 0.5L) == 1.0 &&
208	     (double)(1.0L + vd * 0.5L) == 1.0);
209
210	/*
211	 * FLT_EVAL_METHOD > 1 implies that float expressions are always
212	 * evaluated in double precision or higher, but some compilers get
213	 * this wrong when registers spill to memory.  The following expression
214	 * forces a spill when there are at most 8 FP registers.
215	 */
216	test("implicit promption to double or higher precision is consistent",
217#if FLT_EVAL_METHOD == 1 || FLT_EVAL_METHOD == 2 || defined(__i386__)
218	       TWICE(TWICE(TWICE(TWICE(TWICE(
219	           TWICE(TWICE(TWICE(TWICE(1.0F + vf * 0.5F)))))))))
220	     == (1.0 + FLT_EPSILON * 0.5) * 512.0
221#else
222	     1
223#endif
224	    );
225
226	f = 1.0 + FLT_EPSILON * 0.5;
227	d = 1.0L + DBL_EPSILON * 0.5L;
228	test("const assignment discards extra precision", f == 1.0F && d == 1.0);
229
230	f = 1.0 + vf * 0.5;
231	d = 1.0L + vd * 0.5L;
232	test("variable assignment discards explicit extra precision",
233	     f == 1.0F && d == 1.0);
234	f = 1.0F + vf * 0.5F;
235	d = 1.0 + vd * 0.5;
236	test("variable assignment discards implicit extra precision",
237	     f == 1.0F && d == 1.0);
238
239	test("return discards extra precision",
240	     tofloat(1.0 + vf * 0.5) == 1.0F &&
241	     todouble(1.0L + vd * 0.5L) == 1.0);
242
243	fesetround(FE_UPWARD);
244	/* XXX disabled (works with -frounding-math) */
245	skiptest("FENV_ACCESS: constant arithmetic respects rounding mode",
246	    1.0F + FLT_MIN == 1.0F + FLT_EPSILON &&
247	    1.0 + DBL_MIN == 1.0 + DBL_EPSILON &&
248	    1.0L + LDBL_MIN == 1.0L + LDBL_EPSILON);
249	fesetround(FE_TONEAREST);
250
251	ld = vld * 0.5;
252	test("associativity is respected",
253	     1.0L + ld + (LDBL_EPSILON * 0.5) == 1.0L &&
254	     1.0L + (LDBL_EPSILON * 0.5) + ld == 1.0L &&
255	     ld + 1.0 + (LDBL_EPSILON * 0.5) == 1.0L &&
256	     ld + (LDBL_EPSILON * 0.5) + 1.0 == 1.0L + LDBL_EPSILON);
257}
258
259int
260main(int argc, char *argv[])
261{
262
263	printf("1..26\n");
264
265#ifdef  __i386__
266	fpsetprec(FP_PE);
267#endif
268	run_tests();
269
270	return (failures);
271}
272