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 <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 }