1 #include <assert.h> 2 #include <fenv.h> 3 #include <math.h> 4 #include <signal.h> 5 #include <stdlib.h> 6 #include <stdio.h> 7 #include <string.h> 8 #include <machine/frame.h> 9 10 int max_error = 4; 11 #include "common.h" 12 13 14 /* maximum allowed FP difference for our tests */ 15 #define EPSILON 0.00000000023283064365386962890625 /* 2^(-32) */ 16 17 #define ERR(x, y) e(__LINE__, (x), (y)) 18 19 static void signal_handler(int signum) 20 { 21 struct sigframe_sigcontext *sigframe; 22 23 /* report signal */ 24 sigframe = (struct sigframe_sigcontext *) ((char *) &signum - 25 (char *) &((struct sigframe_sigcontext *) NULL)->sf_signum); 26 printf("Signal %d at 0x%x\n", signum, sigframe->sf_scp->sc_eip); 27 28 /* count as error */ 29 e(0); 30 fflush(stdout); 31 32 /* handle signa again next time */ 33 signal(signum, signal_handler); 34 } 35 36 static void test_fpclassify(double value, int class, int test_sign) 37 { 38 /* test fpclassify */ 39 if (fpclassify(value) != class) e(101); 40 if (test_sign) 41 { 42 if (fpclassify(-value) != class) e(102); 43 44 /* test signbit */ 45 if (signbit(value)) e(103); 46 if (!signbit(-value)) e(104); 47 } 48 } 49 50 /* Maximum normal double: (2 - 2^(-53)) * 2^1023 */ 51 #define NORMAL_DOUBLE_MAX 1.797693134862315708145274237317e+308 52 53 /* Minimum normal double: 2^(-1022) */ 54 #define NORMAL_DOUBLE_MIN 2.2250738585072013830902327173324e-308 55 56 /* Maximum subnormal double: (1 - 2^(-53)) * 2^(-1022) */ 57 #define SUBNORMAL_DOUBLE_MAX 2.2250738585072008890245868760859e-308 58 59 /* Minimum subnormal double: 2^(-52) * 2^(-1023) */ 60 #define SUBNORMAL_DOUBLE_MIN 4.9406564584124654417656879286822e-324 61 62 static void test_fpclassify_values(void) 63 { 64 double d; 65 char negzero[] = { 0, 0, 0, 0, 0, 0, 0, 0x80 }; 66 67 /* test some corner cases for fpclassify and signbit */ 68 test_fpclassify(INFINITY, FP_INFINITE, 1); 69 test_fpclassify(NAN, FP_NAN, 0); 70 test_fpclassify(0, FP_ZERO, 0); 71 test_fpclassify(1, FP_NORMAL, 1); 72 test_fpclassify(NORMAL_DOUBLE_MAX, FP_NORMAL, 1); 73 test_fpclassify(NORMAL_DOUBLE_MIN, FP_NORMAL, 1); 74 test_fpclassify(SUBNORMAL_DOUBLE_MAX, FP_SUBNORMAL, 1); 75 test_fpclassify(SUBNORMAL_DOUBLE_MIN, FP_SUBNORMAL, 1); 76 77 /* 78 * unfortunately the minus operator does not change the sign of zero, 79 * so a special case is needed to test it 80 */ 81 assert(sizeof(negzero) == sizeof(double)); 82 test_fpclassify(*(double *) negzero, FP_ZERO, 0); 83 if (!signbit(*(double *) negzero)) e(4); 84 85 /* test other small numbers for fpclassify and signbit */ 86 d = 1; 87 while (d >= NORMAL_DOUBLE_MIN) 88 { 89 test_fpclassify(d, FP_NORMAL, 1); 90 d /= 10; 91 } 92 while (d >= SUBNORMAL_DOUBLE_MIN) 93 { 94 test_fpclassify(d, FP_SUBNORMAL, 1); 95 d /= 10; 96 } 97 test_fpclassify(d, FP_ZERO, 0); 98 99 /* test other large numbers for fpclassify and signbit */ 100 d = 1; 101 while (d <= NORMAL_DOUBLE_MAX / 10) 102 { 103 test_fpclassify(d, FP_NORMAL, 1); 104 d *= 10; 105 } 106 } 107 108 /* expected rounding: up, down or identical */ 109 #define ROUND_EQ 1 110 #define ROUND_DN 2 111 #define ROUND_UP 3 112 113 static void test_round_value_mode_func(double value, int mode, double (*func)(double), int exp) 114 { 115 int mode_old; 116 double rounded; 117 118 /* update and check rounding mode */ 119 mode_old = fegetround(); 120 fesetround(mode); 121 if (fegetround() != mode) e(5); 122 123 /* perform rounding */ 124 rounded = func(value); 125 126 /* check direction of rounding */ 127 switch (exp) 128 { 129 case ROUND_EQ: if (rounded != value) e(6); break; 130 case ROUND_DN: if (rounded >= value) e(7); break; 131 case ROUND_UP: if (rounded <= value) e(8); break; 132 default: assert(0); 133 } 134 135 /* check whether the number is sufficiently close */ 136 if (fabs(value - rounded) >= 1) e(9); 137 138 /* check whether the number is integer */ 139 if (remainder(rounded, 1)) e(10); 140 141 /* re-check and restore rounding mode */ 142 if (fegetround() != mode) e(11); 143 fesetround(mode_old); 144 } 145 146 static void test_round_value_mode(double value, int mode, int exp_nearbyint, 147 int exp_ceil, int exp_floor, int exp_trunc) 148 { 149 /* test both nearbyint and trunc */ 150 #if 0 151 test_round_value_mode_func(value, mode, nearbyint, exp_nearbyint); 152 #endif 153 test_round_value_mode_func(value, mode, ceil, exp_ceil); 154 test_round_value_mode_func(value, mode, floor, exp_floor); 155 test_round_value_mode_func(value, mode, trunc, exp_trunc); 156 } 157 158 static void test_round_value(double value, int exp_down, int exp_near, int exp_up, int exp_trunc) 159 { 160 /* test each rounding mode */ 161 test_round_value_mode(value, FE_DOWNWARD, exp_down, exp_up, exp_down, exp_trunc); 162 test_round_value_mode(value, FE_TONEAREST, exp_near, exp_up, exp_down, exp_trunc); 163 test_round_value_mode(value, FE_UPWARD, exp_up, exp_up, exp_down, exp_trunc); 164 test_round_value_mode(value, FE_TOWARDZERO, exp_trunc, exp_up, exp_down, exp_trunc); 165 } 166 167 static void test_round_values(void) 168 { 169 int i; 170 171 /* test various values */ 172 for (i = -100000; i < 100000; i++) 173 { 174 test_round_value(i + 0.00, ROUND_EQ, ROUND_EQ, ROUND_EQ, ROUND_EQ); 175 test_round_value(i + 0.25, ROUND_DN, ROUND_DN, ROUND_UP, (i >= 0) ? ROUND_DN : ROUND_UP); 176 test_round_value(i + 0.50, ROUND_DN, (i % 2) ? ROUND_UP : ROUND_DN, ROUND_UP, (i >= 0) ? ROUND_DN : ROUND_UP); 177 test_round_value(i + 0.75, ROUND_DN, ROUND_UP, ROUND_UP, (i >= 0) ? ROUND_DN : ROUND_UP); 178 } 179 } 180 181 static void test_remainder_value(double x, double y) 182 { 183 int mode_old; 184 double r1, r2; 185 186 assert(y != 0); 187 188 /* compute remainder using the function */ 189 r1 = remainder(x, y); 190 191 /* compute remainder using alternative approach */ 192 mode_old = fegetround(); 193 fesetround(FE_TONEAREST); 194 r2 = x - rint(x / y) * y; 195 fesetround(mode_old); 196 197 /* Compare results */ 198 if (fabs(r1 - r2) > EPSILON && fabs(r1 + r2) > EPSILON) 199 { 200 printf("%.20g != %.20g\n", r1, r2); 201 e(13); 202 } 203 } 204 205 static void test_remainder_values_y(double x) 206 { 207 int i, j; 208 209 /* try various rational and transcendental values for y (except zero) */ 210 for (i = -50; i <= 50; i++) 211 if (i != 0) 212 for (j = 1; j < 10; j++) 213 { 214 test_remainder_value(x, (double) i / (double) j); 215 test_remainder_value(x, i * M_E + j * M_PI); 216 } 217 } 218 219 static void test_remainder_values_xy(void) 220 { 221 int i, j; 222 223 /* try various rational and transcendental values for x */ 224 for (i = -50; i <= 50; i++) 225 for (j = 1; j < 10; j++) 226 { 227 test_remainder_values_y((double) i / (double) j); 228 test_remainder_values_y(i * M_E + j * M_PI); 229 } 230 } 231 232 int main(int argc, char **argv) 233 { 234 fenv_t fenv; 235 int i; 236 237 start(47); 238 239 /* no FPU errors, please */ 240 if (feholdexcept(&fenv) < 0) e(14); 241 242 /* some signals count as errors */ 243 for (i = 0; i < _NSIG; i++) 244 if (i != SIGINT && i != SIGTERM && i != SIGKILL) 245 signal(i, signal_handler); 246 247 /* test various floating point support functions */ 248 test_fpclassify_values(); 249 test_remainder_values_xy(); 250 test_round_values(); 251 quit(); 252 return -1; 253 } 254