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$"); 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 52 static const float one_f = 1.0 + FLT_EPSILON / 2; 53 static const double one_d = 1.0 + DBL_EPSILON / 2; 54 static const long double one_ld = 1.0L + LDBL_EPSILON / 2; 55 56 static int testnum, failures; 57 58 static void 59 test1(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 */ 72 static int 73 fpequal(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 81 void 82 run_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 93 void 94 run_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 101 static inline double 102 todouble(long double ld) 103 { 104 105 return (ld); 106 } 107 108 static inline float 109 tofloat(double d) 110 { 111 112 return (d); 113 } 114 115 void 116 run_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 259 int 260 main(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