test/test47.c

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