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 #include <fenv.h> 33 #include <float.h> 34 #include <math.h> 35 #include <stdio.h> 36 37 #ifdef __i386__ 38 #include <ieeefp.h> 39 #endif 40 41 #define ALL_STD_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \ 42 FE_OVERFLOW | FE_UNDERFLOW) 43 44 #define TWICE(x) ((x) + (x)) 45 #define test(desc, pass) test1((desc), (pass), 0) 46 #define skiptest(desc, pass) test1((desc), (pass), 1) 47 48 #pragma STDC FENV_ACCESS ON 49 50 static const float one_f = 1.0 + FLT_EPSILON / 2; 51 static const double one_d = 1.0 + DBL_EPSILON / 2; 52 static const long double one_ld = 1.0L + LDBL_EPSILON / 2; 53 54 static int testnum, failures; 55 56 static void 57 test1(const char *testdesc, int pass, int skip) 58 { 59 60 testnum++; 61 printf("%sok %d - %s%s\n", pass || skip ? "" : "not ", testnum, 62 skip ? "(SKIPPED) " : "", testdesc); 63 if (!pass && !skip) 64 failures++; 65 } 66 67 /* 68 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. 69 */ 70 static int 71 fpequal(long double d1, long double d2) 72 { 73 74 if (d1 != d2) 75 return (isnan(d1) && isnan(d2)); 76 return (copysignl(1.0, d1) == copysignl(1.0, d2)); 77 } 78 79 void 80 run_zero_opt_test(double d1, double d2) 81 { 82 83 test("optimizations don't break the sign of 0", 84 fpequal(d1 - d2, 0.0) 85 && fpequal(-d1 + 0.0, 0.0) 86 && fpequal(-d1 - d2, -0.0) 87 && fpequal(-(d1 - d2), -0.0) 88 && fpequal(-d1 - (-d2), 0.0)); 89 } 90 91 void 92 run_inf_opt_test(double d) 93 { 94 95 test("optimizations don't break infinities", 96 fpequal(d / d, NAN) && fpequal(0.0 * d, NAN)); 97 } 98 99 static inline double 100 todouble(long double ld) 101 { 102 103 return (ld); 104 } 105 106 static inline float 107 tofloat(double d) 108 { 109 110 return (d); 111 } 112 113 void 114 run_tests(void) 115 { 116 volatile long double vld; 117 long double ld; 118 volatile double vd; 119 double d; 120 volatile float vf; 121 float f; 122 int x; 123 124 test("sign bits", fpequal(-0.0, -0.0) && !fpequal(0.0, -0.0)); 125 126 vd = NAN; 127 test("NaN equality", fpequal(NAN, NAN) && NAN != NAN && vd != vd); 128 129 feclearexcept(ALL_STD_EXCEPT); 130 test("NaN comparison returns false", !(vd <= vd)); 131 /* 132 * XXX disabled; gcc/amd64 botches this IEEE 754 requirement by 133 * emitting ucomisd instead of comisd. 134 */ 135 skiptest("FENV_ACCESS: NaN comparison raises invalid exception", 136 fetestexcept(ALL_STD_EXCEPT) == FE_INVALID); 137 138 vd = 0.0; 139 run_zero_opt_test(vd, vd); 140 141 vd = INFINITY; 142 run_inf_opt_test(vd); 143 144 feclearexcept(ALL_STD_EXCEPT); 145 vd = INFINITY; 146 x = (int)vd; 147 /* XXX disabled (works with -O0); gcc doesn't support FENV_ACCESS */ 148 skiptest("FENV_ACCESS: Inf->int conversion raises invalid exception", 149 fetestexcept(ALL_STD_EXCEPT) == FE_INVALID); 150 151 /* Raising an inexact exception here is an IEEE-854 requirement. */ 152 feclearexcept(ALL_STD_EXCEPT); 153 vd = 0.75; 154 x = (int)vd; 155 test("0.75->int conversion rounds toward 0, raises inexact exception", 156 x == 0 && fetestexcept(ALL_STD_EXCEPT) == FE_INEXACT); 157 158 feclearexcept(ALL_STD_EXCEPT); 159 vd = -42.0; 160 x = (int)vd; 161 test("-42.0->int conversion is exact, raises no exception", 162 x == -42 && fetestexcept(ALL_STD_EXCEPT) == 0); 163 164 feclearexcept(ALL_STD_EXCEPT); 165 x = (int)INFINITY; 166 /* XXX disabled; gcc doesn't support FENV_ACCESS */ 167 skiptest("FENV_ACCESS: const Inf->int conversion raises invalid", 168 fetestexcept(ALL_STD_EXCEPT) == FE_INVALID); 169 170 feclearexcept(ALL_STD_EXCEPT); 171 x = (int)0.5; 172 /* XXX disabled; gcc doesn't support FENV_ACCESS */ 173 skiptest("FENV_ACCESS: const double->int conversion raises inexact", 174 x == 0 && fetestexcept(ALL_STD_EXCEPT) == FE_INEXACT); 175 176 test("compile-time constants don't have too much precision", 177 one_f == 1.0L && one_d == 1.0L && one_ld == 1.0L); 178 179 test("const minimum rounding precision", 180 1.0F + FLT_EPSILON != 1.0F && 181 1.0 + DBL_EPSILON != 1.0 && 182 1.0L + LDBL_EPSILON != 1.0L); 183 184 /* It isn't the compiler's fault if this fails on FreeBSD/i386. */ 185 vf = FLT_EPSILON; 186 vd = DBL_EPSILON; 187 vld = LDBL_EPSILON; 188 test("runtime minimum rounding precision", 189 1.0F + vf != 1.0F && 1.0 + vd != 1.0 && 1.0L + vld != 1.0L); 190 191 test("explicit float to float conversion discards extra precision", 192 (float)(1.0F + FLT_EPSILON * 0.5F) == 1.0F && 193 (float)(1.0F + vf * 0.5F) == 1.0F); 194 test("explicit double to float conversion discards extra precision", 195 (float)(1.0 + FLT_EPSILON * 0.5) == 1.0F && 196 (float)(1.0 + vf * 0.5) == 1.0F); 197 test("explicit ldouble to float conversion discards extra precision", 198 (float)(1.0L + FLT_EPSILON * 0.5L) == 1.0F && 199 (float)(1.0L + vf * 0.5L) == 1.0F); 200 201 test("explicit double to double conversion discards extra precision", 202 (double)(1.0 + DBL_EPSILON * 0.5) == 1.0 && 203 (double)(1.0 + vd * 0.5) == 1.0); 204 test("explicit ldouble to double conversion discards extra precision", 205 (double)(1.0L + DBL_EPSILON * 0.5L) == 1.0 && 206 (double)(1.0L + vd * 0.5L) == 1.0); 207 208 /* 209 * FLT_EVAL_METHOD > 1 implies that float expressions are always 210 * evaluated in double precision or higher, but some compilers get 211 * this wrong when registers spill to memory. The following expression 212 * forces a spill when there are at most 8 FP registers. 213 */ 214 test("implicit promption to double or higher precision is consistent", 215 #if FLT_EVAL_METHOD == 1 || FLT_EVAL_METHOD == 2 || defined(__i386__) 216 TWICE(TWICE(TWICE(TWICE(TWICE( 217 TWICE(TWICE(TWICE(TWICE(1.0F + vf * 0.5F))))))))) 218 == (1.0 + FLT_EPSILON * 0.5) * 512.0 219 #else 220 1 221 #endif 222 ); 223 224 f = 1.0 + FLT_EPSILON * 0.5; 225 d = 1.0L + DBL_EPSILON * 0.5L; 226 test("const assignment discards extra precision", f == 1.0F && d == 1.0); 227 228 f = 1.0 + vf * 0.5; 229 d = 1.0L + vd * 0.5L; 230 test("variable assignment discards explicit extra precision", 231 f == 1.0F && d == 1.0); 232 f = 1.0F + vf * 0.5F; 233 d = 1.0 + vd * 0.5; 234 test("variable assignment discards implicit extra precision", 235 f == 1.0F && d == 1.0); 236 237 test("return discards extra precision", 238 tofloat(1.0 + vf * 0.5) == 1.0F && 239 todouble(1.0L + vd * 0.5L) == 1.0); 240 241 fesetround(FE_UPWARD); 242 /* XXX disabled (works with -frounding-math) */ 243 skiptest("FENV_ACCESS: constant arithmetic respects rounding mode", 244 1.0F + FLT_MIN == 1.0F + FLT_EPSILON && 245 1.0 + DBL_MIN == 1.0 + DBL_EPSILON && 246 1.0L + LDBL_MIN == 1.0L + LDBL_EPSILON); 247 fesetround(FE_TONEAREST); 248 249 ld = vld * 0.5; 250 test("associativity is respected", 251 1.0L + ld + (LDBL_EPSILON * 0.5) == 1.0L && 252 1.0L + (LDBL_EPSILON * 0.5) + ld == 1.0L && 253 ld + 1.0 + (LDBL_EPSILON * 0.5) == 1.0L && 254 ld + (LDBL_EPSILON * 0.5) + 1.0 == 1.0L + LDBL_EPSILON); 255 } 256 257 int 258 main(int argc, char *argv[]) 259 { 260 261 printf("1..26\n"); 262 263 #ifdef __i386__ 264 fpsetprec(FP_PE); 265 #endif 266 run_tests(); 267 268 return (failures); 269 } 270