newlib/libm/test/math2.c

   1
   2 #include "test.h"
   3 #include <errno.h>
   4
   5
   6 int
   7 randi (void)
   8 {
   9   static int next;
  10   next = (next * 1103515245) + 12345;
  11   return ((next >> 16) & 0xffff);
  12 }
  13
  14 double randx (void)
  15 {
  16   double res;
  17
  18   do
  19   {
  20     union {
  21         short parts[4];
  22         double res;
  23       } u;
  24
  25     u.parts[0] = randi();
  26     u.parts[1] = randi();
  27     u.parts[2] = randi();
  28     u.parts[3] = randi();
  29     res = u.res;
  30
  31   } while (!finite(res));
  32
  33   return res ;
  34 }
  35
  36 /* Return a random double, but bias for numbers closer to 0 */
  37 double randy (void)
  38 {
  39   int pow;
  40   double r= randx();
  41   r = frexp(r, &pow);
  42   return ldexp(r, randi() & 0x1f);
  43 }
  44
  45 void
  46 test_frexp (void)
  47 {
  48   int i;
  49   double r;
  50   int t;
  51
  52   float xf;
  53   double gives;
  54
  55   int pow;
  56
  57
  58   /* Frexp of x return a and n, where a * 2**n == x, so test this with a
  59      set of random numbers */
  60   for (t = 0; t < 2; t++)
  61   {
  62     for (i = 0; i < 1000; i++)
  63     {
  64
  65       double x = randx();
  66       line(i);
  67       switch (t)
  68       {
  69       case 0:
  70         newfunc("frexp/ldexp");
  71         r = frexp(x, &pow);
  72         if (r > 1.0 || r < -1.0)
  73         {
  74           /* Answer can never be > 1 or < 1 */
  75           test_iok(0,1);
  76         }
  77
  78         gives = ldexp(r ,pow);
  79         test_mok(gives,x,62);
  80         break;
  81       case 1:
  82         newfunc("frexpf/ldexpf");
  83         if (x > FLT_MIN && x < FLT_MAX)
  84         {
  85           /* test floats too, but they have a smaller range so make sure x
  86              isn't too big. Also x can get smaller than a float can
  87              represent to make sure that doesn't happen too */
  88           xf = x;
  89           r = frexpf(xf, &pow);
  90           if (r > 1.0 || r < -1.0)
  91           {
  92             /* Answer can never be > 1 or < -1 */
  93             test_iok(0,1);
  94           }
  95
  96           gives = ldexpf(r ,pow);
  97           test_mok(gives,x, 32);
  98
  99         }
 100       }
 101
 102     }
 103
 104   }
 105
 106   /* test a few numbers manually to make sure frexp/ldexp are not
 107      testing as ok because both are broken */
 108
 109   r = frexp(64.0, &i);
 110
 111   test_mok(r, 0.5,64);
 112   test_iok(i, 7);
 113
 114   r = frexp(96.0, &i);
 115
 116   test_mok(r, 0.75, 64);
 117   test_iok(i, 7);
 118
 119 }
 120
 121 /* Test mod - this is given a real hammering by the strtod type
 122    routines, here are some more tests.
 123
 124    By definition
 125
 126    modf = func(value, &iptr)
 127
 128       (*iptr + modf) == value
 129
 130    we test this
 131
 132 */
 133 void
 134 test_mod (void)
 135 {
 136   int i;
 137
 138   newfunc("modf");
 139
 140
 141   for (i = 0; i < 1000; i++)
 142   {
 143     double intpart;
 144     double n;
 145     line(i);
 146     n  = randx();
 147     if (finite(n) && n != 0.0 )
 148     {
 149       double r = modf(n, &intpart);
 150       line(i);
 151       test_mok(intpart + r, n, 63);
 152     }
 153
 154   }
 155   newfunc("modff");
 156
 157   for (i = 0; i < 1000; i++)
 158   {
 159     float intpart;
 160     double nd;
 161     line(i);
 162     nd  = randx() ;
 163     if (nd < FLT_MAX && finitef(nd) && nd != 0.0)
 164     {
 165       float n = nd;
 166       double r = modff(n, &intpart);
 167       line(i);
 168       test_mok(intpart + r, n, 32);
 169     }
 170   }
 171
 172
 173 }
 174
 175 /*
 176 Test pow by multiplying logs
 177 */
 178 void
 179 test_pow (void)
 180 {
 181   unsigned int i;
 182   newfunc("pow");
 183
 184   for (i = 0; i < 1000; i++)
 185   {
 186     double n1;
 187     double n2;
 188     double res;
 189     double shouldbe;
 190
 191     line(i);
 192     n1 = fabs(randy());
 193     n2 = fabs(randy()/100.0);
 194     res = pow(n1, n2);
 195     shouldbe = exp(log(n1) * n2);
 196     test_mok(shouldbe, res,64);
 197   }
 198
 199   newfunc("powf");
 200
 201   for (i = 0; i < 1000; i++)
 202   {
 203     double n1;
 204     double n2;
 205     double res;
 206     double shouldbe;
 207
 208     errno = 0;
 209
 210     line(i);
 211     n1 = fabs(randy());
 212     n2 = fabs(randy()/100.0);
 213     res = powf(n1, n2);
 214     shouldbe = expf(logf(n1) * n2);
 215     if (!errno)
 216      test_mok(shouldbe, res,28);
 217   }
 218
 219
 220
 221
 222 }
 223
 224
 225
 226 void
 227 test_math2 (void)
 228 {
 229   test_mod();
 230   test_frexp();
 231   test_pow();
 232 }