newlib/libm/test/math.c

   1 /*
   2   Test the library maths functions using trusted precomputed test
   3   vectors.
   4
   5   These vectors were originally generated on a sun3 with a 68881 using
   6   80 bit precision, but ...
   7
   8   Each function is called with a variety of interesting arguments.
   9   Note that many of the polynomials we use behave badly when the
  10   domain is stressed, so the numbers in the vectors depend on what is
  11   useful to test - eg sin(1e30) is pointless - the arg has to be
  12   reduced modulo pi, and after that there's no bits of significance
  13   left to evaluate with - any number would be just as precise as any
  14   other.
  15
  16
  17 */
  18
  19 #include "test.h"
  20 #include <math.h>
  21 #include <ieeefp.h>
  22 #include <float.h>
  23 #include <math.h>
  24 #include <errno.h>
  25 #include <stdio.h>
  26
  27 int merror;
  28 double mretval = 64;
  29 int traperror = 1;
  30 char *mname;
  31
  32 void translate_to (FILE *file,
  33             double r)
  34 {
  35   __ieee_double_shape_type bits;
  36   bits.value = r;
  37   fprintf(file, "0x%08x, 0x%08x", bits.parts.msw, bits.parts.lsw);
  38 }
  39
  40 int
  41 ffcheck (double is,
  42        one_line_type *p,
  43        char *name,
  44        int serrno,
  45        int merror)
  46 {
  47   /* Make sure the answer isn't to far wrong from the correct value */
  48   __ieee_double_shape_type correct, isbits;
  49   int mag;
  50   isbits.value = is;
  51
  52   correct.parts.msw = p->qs[0].msw;
  53   correct.parts.lsw = p->qs[0].lsw;
  54
  55   mag = mag_of_error(correct.value, is);
  56
  57   if (mag < p->error_bit)
  58   {
  59     inacc ++;
  60
  61     printf("%s:%d, inaccurate answer: bit %d (%08x%08x %08x%08x) (%g %g)\n",
  62            name,  p->line, mag,
  63            correct.parts.msw,
  64            correct.parts.lsw,
  65            isbits.parts.msw,
  66            isbits.parts.lsw,
  67            correct.value, is);
  68   }
  69
  70 #if 0
  71   if (p->qs[0].merror != merror)
  72   {
  73     /* Beware, matherr doesn't exist anymore.  */
  74     printf("testing %s_vec.c:%d, matherr wrong: %d %d\n",
  75            name, p->line, merror, p->qs[0].merror);
  76   }
  77
  78   if (p->qs[0].errno_val != errno)
  79   {
  80     printf("testing %s_vec.c:%d, errno wrong: %d %d\n",
  81            name, p->line, errno, p->qs[0].errno_val);
  82
  83   }
  84 #endif
  85   return mag;
  86 }
  87
  88 double
  89 thedouble (long msw,
  90        long lsw)
  91 {
  92   __ieee_double_shape_type x;
  93
  94   x.parts.msw = msw;
  95   x.parts.lsw = lsw;
  96   return x.value;
  97 }
  98
  99 int calc;
 100
 101
 102 frontline (FILE *f,
 103        int mag,
 104        one_line_type *p,
 105        double result,
 106        int merror,
 107        int errno,
 108        char *args,
 109        char *name)
 110 {
 111   if (reduce && p->error_bit < mag)
 112   {
 113     fprintf(f, "{%2d,", p->error_bit);
 114   }
 115   else
 116   {
 117     fprintf(f, "{%2d,",mag);
 118   }
 119
 120
 121   fprintf(f,"%2d,%3d,", merror,errno);
 122   fprintf(f, "__LINE__, ");
 123
 124   if (calc)
 125   {
 126     translate_to(f, result);
 127   }
 128   else
 129   {
 130     translate_to(f, thedouble(p->qs[0].msw, p->qs[0].lsw));
 131   }
 132
 133   fprintf(f, ", ");
 134
 135   fprintf(f,"0x%08x, 0x%08x", p->qs[1].msw, p->qs[1].lsw);
 136
 137
 138   if (args[2])
 139   {
 140     fprintf(f, ", ");
 141     fprintf(f,"0x%08x, 0x%08x", p->qs[2].msw, p->qs[2].lsw);
 142   }
 143
 144   fprintf(f,"}, /* %g=f(%g",result,
 145           thedouble(p->qs[1].msw, p->qs[1].lsw));
 146
 147   if (args[2])
 148   {
 149     fprintf(f,", %g", thedouble(p->qs[2].msw,p->qs[2].lsw));
 150   }
 151   fprintf(f, ")*/\n");
 152 }
 153
 154 finish (FILE *f,
 155        int vector,
 156        double result,
 157        one_line_type *p,
 158        char *args,
 159        char *name)
 160 {
 161   int mag;
 162
 163   mag = ffcheck(result, p,name,  merror, errno);
 164   if (vector)
 165   {
 166     frontline(f, mag, p, result, merror, errno, args , name);
 167   }
 168 }
 169
 170 run_vector_1 (int vector,
 171        one_line_type *p,
 172        char *func,
 173        char *name,
 174        char *args)
 175 {
 176   FILE *f;
 177   int mag;
 178   double result;
 179
 180   if (vector)
 181   {
 182
 183     VECOPEN(name, f);
 184
 185     if (redo)
 186     {
 187       double k;
 188
 189       for (k = -.2; k < .2; k+= 0.00132)
 190       {
 191
 192         fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
 193                 k,k+4);
 194
 195       }
 196
 197       for (k = -1.2; k < 1.2; k+= 0.01)
 198       {
 199
 200         fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
 201                 k,k+4);
 202
 203       }
 204       for (k = -M_PI *2; k < M_PI *2; k+= M_PI/2)
 205       {
 206
 207         fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
 208                 k,k+4);
 209
 210       }
 211
 212       for (k = -30; k < 30; k+= 1.7)
 213       {
 214
 215         fprintf(f,"{2,2, 1,1, 0,0, 0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
 216                 k,k+4);
 217
 218       }
 219       VECCLOSE(f, name, args);
 220       return;
 221     }
 222   }
 223
 224   newfunc(name);
 225   while (p->line)
 226   {
 227     double arg1 = thedouble(p->qs[1].msw, p->qs[1].lsw);
 228     double arg2 = thedouble(p->qs[2].msw, p->qs[2].lsw);
 229
 230     double r;
 231     double rf;
 232
 233     errno = 0;
 234     merror = 0;
 235     mname = 0;
 236
 237
 238     line(p->line);
 239
 240     merror = 0;
 241     errno = 123;
 242
 243     if (strcmp(args,"dd")==0)
 244     {
 245       typedef double (*pdblfunc) (double);
 246
 247       /* Double function returning a double */
 248
 249       result = ((pdblfunc)(func))(arg1);
 250
 251       finish(f,vector, result, p, args, name);
 252     }
 253     else  if (strcmp(args,"ff")==0)
 254     {
 255       float arga;
 256       double a;
 257
 258       typedef float (*pdblfunc) (float);
 259
 260       /* Double function returning a double */
 261
 262       if (arg1 < FLT_MAX )
 263       {
 264         arga = arg1;
 265         result = ((pdblfunc)(func))(arga);
 266         finish(f, vector, result, p,args, name);
 267       }
 268     }
 269     else if (strcmp(args,"ddd")==0)
 270      {
 271        typedef double (*pdblfunc) (double,double);
 272
 273        result = ((pdblfunc)(func))(arg1,arg2);
 274        finish(f, vector, result, p,args, name);
 275      }
 276      else  if (strcmp(args,"fff")==0)
 277      {
 278        double a,b;
 279
 280        float arga;
 281        float argb;
 282
 283        typedef float (*pdblfunc) (float,float);
 284
 285
 286        if (arg1 < FLT_MAX && arg2 < FLT_MAX)
 287        {
 288          arga = arg1;
 289          argb = arg2;
 290          result = ((pdblfunc)(func))(arga, argb);
 291          finish(f, vector, result, p,args, name);
 292        }
 293      }
 294      else if (strcmp(args,"did")==0)
 295      {
 296        typedef double (*pdblfunc) (int,double);
 297
 298        result = ((pdblfunc)(func))((int)arg1,arg2);
 299        finish(f, vector, result, p,args, name);
 300      }
 301      else  if (strcmp(args,"fif")==0)
 302      {
 303        double a,b;
 304
 305        float arga;
 306        float argb;
 307
 308        typedef float (*pdblfunc) (int,float);
 309
 310
 311        if (arg1 < FLT_MAX && arg2 < FLT_MAX)
 312        {
 313          arga = arg1;
 314          argb = arg2;
 315          result = ((pdblfunc)(func))((int)arga, argb);
 316          finish(f, vector, result, p,args, name);
 317        }
 318      }
 319
 320     p++;
 321   }
 322   if (vector)
 323   {
 324     VECCLOSE(f, name, args);
 325   }
 326 }
 327
 328 void
 329 test_math (void)
 330 {
 331   test_acos(0);
 332   test_acosf(0);
 333   test_acosh(0);
 334   test_acoshf(0);
 335   test_asin(0);
 336   test_asinf(0);
 337   test_asinh(0);
 338   test_asinhf(0);
 339   test_atan(0);
 340   test_atan2(0);
 341   test_atan2f(0);
 342   test_atanf(0);
 343   test_atanh(0);
 344   test_atanhf(0);
 345   test_ceil(0);
 346   test_ceilf(0);
 347   test_cos(0);
 348   test_cosf(0);
 349   test_cosh(0);
 350   test_coshf(0);
 351   test_erf(0);
 352   test_erfc(0);
 353   test_erfcf(0);
 354   test_erff(0);
 355   test_exp(0);
 356   test_expf(0);
 357   test_fabs(0);
 358   test_fabsf(0);
 359   test_floor(0);
 360   test_floorf(0);
 361   test_fmod(0);
 362   test_fmodf(0);
 363   test_gamma(0);
 364   test_gammaf(0);
 365   test_hypot(0);
 366   test_hypotf(0);
 367   test_j0(0);
 368   test_j0f(0);
 369   test_j1(0);
 370   test_j1f(0);
 371   test_jn(0);
 372   test_jnf(0);
 373   test_log(0);
 374   test_log10(0);
 375   test_log10f(0);
 376   test_log1p(0);
 377   test_log1pf(0);
 378   test_log2(0);
 379   test_log2f(0);
 380   test_logf(0);
 381   test_sin(0);
 382   test_sinf(0);
 383   test_sinh(0);
 384   test_sinhf(0);
 385   test_sqrt(0);
 386   test_sqrtf(0);
 387   test_tan(0);
 388   test_tanf(0);
 389   test_tanh(0);
 390   test_tanhf(0);
 391   test_y0(0);
 392   test_y0f(0);
 393   test_y1(0);
 394   test_y1f(0);
 395   test_y1f(0);
 396   test_ynf(0);
 397 }
 398
 399 /* These have to be played with to get to compile on machines which
 400    don't have the fancy <foo>f entry points
 401 */
 402
 403 #if 0
 404 float cosf (float a) { return cos((double)a); }
 405 float sinf (float  a) { return sin((double)a); }
 406 float log1pf (float a) { return log1p((double)a); }
 407 float tanf (float a) { return tan((double)a); }
 408 float ceilf (float a) { return ceil(a); }
 409 float floorf (float a) { return floor(a); }
 410 #endif
 411
 412 /*ndef HAVE_FLOAT*/
 413 #if 0
 414
 415 float fmodf(a,b) float a,b; { return fmod(a,b); }
 416 float hypotf(a,b) float a,b; { return hypot(a,b); }
 417
 418 float acosf(a) float a; { return acos(a); }
 419 float acoshf(a) float a; { return acosh(a); }
 420 float asinf(a) float a; { return asin(a); }
 421 float asinhf(a) float a; { return asinh(a); }
 422 float atanf(a) float a; { return atan(a); }
 423 float atanhf(a) float a; { return atanh(a); }
 424
 425 float coshf(a) float a; { return cosh(a); }
 426 float erff(a) float a; { return erf(a); }
 427 float erfcf(a) float a; { return erfc(a); }
 428 float expf(a) float a; { return exp(a); }
 429 float fabsf(a) float a; { return fabs(a); }
 430
 431 float gammaf(a) float a; { return gamma(a); }
 432 float j0f(a) float a; { return j0(a); }
 433 float j1f(a) float a; { return j1(a); }
 434 float log10f(a) float a; { return log10(a); }
 435
 436 float logf(a) float a; { return log(a); }
 437
 438 float sinhf(a) float a; { return sinh(a); }
 439 float sqrtf(a) float a; { return sqrt(a); }
 440
 441 float tanhf(a) float a; { return tanh(a); }
 442 float y0f(a) float a; { return y0(a); }
 443 float y1f(a) float a; { return y1(a); }
 444 #endif