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