math.c

edb9315 wince4.2最新的bsp代码

/* 
  Test the library maths functions using trusted precomputed test
  vectors.

  These vectors were originally generated on a sun3 with a 68881 using
  80 bit precision, but ...

  Each function is called with a variety of interesting arguments.
  Note that many of the polynomials we use behave badly when the
  domain is stressed, so the numbers in the vectors depend on what is
  useful to test - eg sin(1e30) is pointless - the arg has to be
  reduced modulo pi, and after that there's no bits of significance
  left to evaluate with - any number would be just as precise as any
  other.


*/
#include "test.h"
#include <stdio.h>
#include <stdlib.h>

int inacc;

int merror;
double mretval = 64;
int traperror = 1;
char *mname;

int verbose;


void _DEFUN(translate_to,(file,r),
	    FILE *file _AND
	    double r)
{
  __ieee_double_shape_type bits;
  bits.value = r;
  fprintf(file, "0x%08x, 0x%08x", bits.parts.msw, bits.parts.lsw);
}


int 
_DEFUN(bigger,(a,b),
	   __ieee_double_shape_type *a  _AND
	   __ieee_double_shape_type *b)
{

  if (a->parts.msw > b->parts.msw) 
    {

      return 1;
    } 
  else if (a->parts.msw == b->parts.msw) 
    {
      if (a->parts.lsw > b->parts.lsw) 
	{
	  return 1;
	}
    }
  return 0;
}

/* Return the first bit different between two double numbers */
int 
_DEFUN(mag_of_error,(is, shouldbe),
       double is _AND
       double shouldbe)
{
  __ieee_double_shape_type a,b;
  int i;
  int a_big;
  unsigned  int mask;
  unsigned long int __x;
  unsigned long int msw, lsw;						  
  a.value = is;
  
  b.value = shouldbe;
  
  if (a.parts.msw == b.parts.msw 
      && a.parts.lsw== b.parts.lsw) return 64;


  /* Subtract the larger from the smaller number */

  a_big = bigger(&a, &b);

  if (!a_big) {
    int t;
    t = a.parts.msw;
    a.parts.msw = b.parts.msw;
    b.parts.msw = t;

    t = a.parts.lsw;
    a.parts.lsw = b.parts.lsw;
    b.parts.lsw = t;
  }



  __x = (a.parts.lsw) - (b.parts.lsw);							
  msw = (a.parts.msw) - (b.parts.msw) - (__x > (a.parts.lsw));
  lsw = __x;								

  


  /* Find out which bit the difference is in */
  mask = 0x80000000;
  for (i = 0; i < 32; i++)
  {
    if (((msw) & mask)!=0) return i;
    mask >>=1;
  }
  
  mask = 0x80000000;
  for (i = 0; i < 32; i++)
  {
    
    if (((lsw) & mask)!=0) return i+32;
    mask >>=1;
  }
  
  return 64;
  
}

int 
_DEFUN(ffcheck,( is, p, name, serrno, merror),
       double is _AND
       one_line_type *p _AND
       char *name _AND
       int serrno _AND
       int merror)
{
  /* Make sure the answer isn't to far wrong from the correct value */
  __ieee_double_shape_type correct, isbits;
  int mag;  
  isbits.value = is;
  
  correct.parts.msw = p->qs[0].msw;
  correct.parts.lsw = p->qs[0].lsw;
  
  mag = mag_of_error(correct.value, is);
  
  if (mag < p->error_bit)
  {
    inacc ++;
    
    printf("%s:%d, inaccurate answer: bit %d (%08x%08x %08x%08x) (%g %g)\n",
	   name,  p->line, mag,
	   correct.parts.msw,
	   correct.parts.lsw,
	   isbits.parts.msw,
	   isbits.parts.lsw,
	   correct.value, is);
  }      
  
  return mag;
}

double
_DEFUN(thedouble, (msw, lsw),
       long msw _AND
       long lsw)
{
  __ieee_double_shape_type x;
  
  x.parts.msw = msw;
  x.parts.lsw = lsw;
  return x.value;
}

int calc;
int reduce;


_DEFUN(frontline,(f, mag, p, result, merror, errno, args, name),
       FILE *f _AND
       int mag _AND
       one_line_type *p _AND
       double result _AND
       int merror _AND
       int errno _AND
       char *args _AND
       char *name)
{
  if (reduce && p->error_bit < mag) 
  {
    fprintf(f, "{%2d,", p->error_bit);
  }
  else 
  {
    fprintf(f, "{%2d,",mag);
  }


  fprintf(f,"%2d,%3d,", merror,errno);
  fprintf(f, "__LINE__, ");

  if (calc) 
  {
    translate_to(f, result);
  }
  else 
  {
    translate_to(f, thedouble(p->qs[0].msw, p->qs[0].lsw));
  }
  
  fprintf(f, ", ");      

  fprintf(f,"0x%08x, 0x%08x", p->qs[1].msw, p->qs[1].lsw);
  

  if (args[2]) 
  {
    fprintf(f, ", ");      
    fprintf(f,"0x%08x, 0x%08x", p->qs[2].msw, p->qs[2].lsw);
  }
	
  fprintf(f,"},	/* %g=f(%g",result,
  	  thedouble(p->qs[1].msw, p->qs[1].lsw));

  if (args[2])
  {
    fprintf(f,", %g", thedouble(p->qs[2].msw,p->qs[2].lsw));
  }
  fprintf(f, ")*/\n");      
}

_DEFUN(finish,(f, vector,  result , p, args, name),
       FILE *f _AND
       int vector _AND
       double result _AND
       one_line_type *p _AND
       char *args _AND
       char *name)
{
  int mag;

  mag = ffcheck(result, p,name,  merror, 0);    
  if (vector) 
  {    
    frontline(f, mag, p, result, merror, 0, args , name);
  }
} 

static _CONST char *iname = "foo";
void 
_DEFUN(newfunc,(string),
       _CONST char *string)
{
  if (strcmp(iname, string)) 
  {
    printf("testing %s\n", string);
    fflush(stdout);
    iname = string;
  }
  
}


static int theline;
static int count;

void line(li)
int li;
{
  if (verbose)  
  {
    printf("  %d\n", li);
  }
  theline = li;
  
  count++;
}

int redo;  

void _DEFUN(run_vector_1,(vector, p, func, name, args),
       int vector _AND
       one_line_type *p _AND
       char *func _AND
       char *name _AND
       char *args)
{
  FILE *f;
//  int mag;
  double result;  
  
  if (vector)
  {

    VECOPEN(name, f);

    if (redo)
    {
      double k;

      for (k = -.2; k < .2; k+= 0.00132) 
      {

	fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
		k,k+4);

      }

      for (k = -1.2; k < 1.2; k+= 0.01) 
      {

	fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
		k,k+4);

      }
      for (k = -M_PI *2; k < M_PI *2; k+= M_PI/2) 
      {

	fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
		k,k+4);

      }

      for (k = -30; k < 30; k+= 1.7) 
      {

	fprintf(f,"{2,2, 1,1, 0,0, 0x%08x,0x%08x, 0x%08x, 0x%08x},\n",
		k,k+4);

      }
      VECCLOSE(f, name, args);
      return;
    }
  }
 
  newfunc(name);
  while (p->line) 
  {
    double arg1 = thedouble(p->qs[1].msw, p->qs[1].lsw);
    double arg2 = thedouble(p->qs[2].msw, p->qs[2].lsw);

//    double r;
//    double rf;
    
    merror = 0;
    mname = 0;

    
    line(p->line);          

    merror = 0;

    if (strcmp(args,"dd")==0)
    {
      typedef double _EXFUN((*pdblfunc),(double));
      
      /* Double function returning a double */
      
      result = ((pdblfunc)(func))(arg1);
      
      finish(f,vector, result, p, args, name);       
    }  
    else  if (strcmp(args,"ff")==0)
    {
      float arga;
//      double a;
      
      typedef float (*pdblfuncf)(float);		/*bug fixed*/
      pdblfuncf funcf=(pdblfuncf)(func);
      
      /* Double function returning a double */
      
      if (arg1 < FLT_MAX )
      {
	arga = (float)arg1;      
//	result = ((pdblfunc)(func))(arga);
	result = (funcf)(arga);
	finish(f, vector, result, p,args, name);       
      }
    }      
    else if (strcmp(args,"ddd")==0)
     {
       typedef double _EXFUN((*pdblfunc),(double,double));
      
       result = ((pdblfunc)(func))(arg1,arg2);
       finish(f, vector, result, p,args, name);       
     }  
     else  if (strcmp(args,"fff")==0)
     {
//       double a,b;
       
       float arga;
       float argb;
      
      typedef float (*pdblfuncff)(float,float);		/*bug fixed*/
      pdblfuncff funcf=(pdblfuncff)(func);
      

       if (arg1 < FLT_MAX && arg2 < FLT_MAX) 
       {
	 arga = (float)arg1;      
	 argb = (float)arg2;
	 result = funcf(arga, argb);
	 finish(f, vector, result, p,args, name);       
       }
     }      
     else if (strcmp(args,"did")==0)
     {
       typedef double _EXFUN((*pdblfunc),(int,double));
      
       result = ((pdblfunc)(func))((int)arg1,arg2);
       finish(f, vector, result, p,args, name);       
     }  
     else  if (strcmp(args,"fif")==0)
     {
//       double a,b;
       
       float arga;
       float argb;
      
       typedef float _EXFUN((*pdblfunc),(int,float));
      

       if (arg1 < FLT_MAX && arg2 < FLT_MAX) 
       {
	 arga = (float)arg1;      
	 argb = (float)arg2;
	 result = ((pdblfunc)(func))((int)arga, argb);
	 finish(f, vector, result, p,args, name);       
       }
     }      

    p++;
  }
  if (vector)
  {
    VECCLOSE(f, name, args);
  }
  printf("finished\n");
}

int
_DEFUN_VOID(randi)
{
  static int next;
  next = (next * 1103515245) + 12345;
  return ((next >> 16) & 0xffff);
}

double _DEFUN_VOID(randx)
{
  double res;
  
  do 
  {
    union {
	short parts[4];
	double res;
      } u;
    
    u.parts[0] = randi();
    u.parts[1] = randi();
    u.parts[2] = randi();
    u.parts[3] = randi();
    res = u.res;
    
  } while (!_finite(res));
  
  return res ;
}

/* Return a random double, but bias for numbers closer to 0 */
double _DEFUN_VOID(randy)
{
  int pow;
  double r= randx();
  r = frexp(r, &pow);
  return ldexp(r, randi() & 0x1f);
}


void
_DEFUN(test_mok,(value, shouldbe, okmag),
       double value _AND
       double shouldbe _AND
       int okmag)
{
  __ieee_double_shape_type a,b;
  int mag = mag_of_error(value, shouldbe);
  if (mag == 0) 
  {
    /* error in the first bit is ok if the numbers are both 0 */
    if (value == 0.0 && shouldbe == 0.0)
     return;
    
  }
  a.value = shouldbe;
  b.value = value;
  
  if (mag < okmag) 
  {
    printf("%s:%d, wrong answer: bit %d ",
	   iname, 
	   theline,
	   mag);
     printf("%08x%08x %08x%08x) ",
	    a.parts.msw,	     a.parts.lsw,
	    b.parts.msw,	     b.parts.lsw);
    printf("(%g %g)\n",   a.value, b.value);
    inacc++;
  }
}


/*
Test pow by multiplying logs  
*/
void
_DEFUN_VOID(test_pow)
{
  unsigned int i;  
  newfunc("pow");

  for (i = 0; i < 1000; i++) 
  {
    double n1;
    double n2;
    double res;
    double shouldbe;

    line(i);  
    n1 = fabs(randy());
    n2 = fabs(randy()/100.0);
    res = pow(n1, n2);
    shouldbe = exp(log(n1) * n2);
    test_mok(shouldbe, res,64);
  }
}

extern void test_acos(int);
extern void test_asin(int);
extern void test_atan(int);
extern void test_atan2(int);
extern void test_ceil(int);
extern void test_ceilf(int);
extern void test_cos(int);
extern void test_cosh(int);
extern void test_exp(int);
extern void test_fabs(int);
extern void test_fabsf(int);
extern void test_floor(int);
extern void test_floorf(int);
extern void test_fmod(int);
extern void test_fmodf(int);
extern void test_log(int);
extern void test_log10(int);
extern void test_sin(int);
extern void test_sinh(int);
extern void test_sqrt(int);
extern void test_sqrtf(int);
extern void test_tan(int);
extern void test_tanh(int);
extern void test_y0(int);
extern void test_y1(int);
extern void test_yn(int);
extern void test_pow(int);


void
_DEFUN_VOID(test_math_ce)
{
  test_acos(0);
  test_asin(0);
  test_atan(0);
  test_atan2(0);
  test_ceil(0);
  test_ceilf(0);
  test_cos(0);
  test_cosh(0);
  test_exp(0);
  test_fabs(0);
  test_fabsf(0);
  test_floor(0);
  test_floorf(0);
  test_fmod(0);
  test_fmodf(0);
  test_log(0);
  test_log10(0);
  test_sin(0);
  test_sinh(0);
  test_sqrt(0);
  test_sqrtf(0);
  test_tan(0);
  test_tanh(0);
  test_y0(0);
  test_y1(0);
  test_yn(0);
  test_pow(0);
	
}

#ifdef USE_CRUNCH_LIBRARY
extern void  test_acosf(int);
extern void  test_asinf(int);
extern void  test_atanf(int);
extern void  test_atan2f(int);
extern void  test_cosf(int);
extern void  test_coshf(int);
extern void  test_expf(int);
extern void  test_logf(int);
extern void  test_log10f(int);
extern void  test_sinf(int);
extern void  test_sinhf(int);
extern void  test_tanf(int);
extern void  test_tanhf(int);

void
_DEFUN_VOID(test_math_new)
{
  test_acosf(0);
  test_asinf(0);
  test_atanf(0);
  test_atan2f(0);
  test_cosf(0);
  test_coshf(0);
  test_expf(0);
  test_logf(0);
  test_log10f(0);
  test_sinf(0);
  test_sinhf(0);
  test_tanf(0);
  test_tanhf(0);
}

#endif

void
_DEFUN_VOID(test_math)
{
	test_math_ce();
#ifdef USE_CRUNCH_LIBRARY
	//test_math_new();
#endif
}