validate.c

任意精度的数学库

if (pass)
  fprintf(stdout,"... COS function passes\n");
  
/*
 *
 */

fprintf(stdout,"Spot check the SIN_COS function ... \n");

xx        = -LPI;
xdelta    = 0.475;
pass      = TRUE;
tolerance = 1.0E-14;

while (1)
  {
   yy_sin = sin(xx);
   yy_cos = cos(xx);

   m_apm_set_double(aa_mapm, xx);
   m_apm_sin_cos(bb_mapm, cc_mapm, DECIMAL_PLACES, aa_mapm);

   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);

   ydiff = yy_mapm_lib - yy_sin;

   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx,yy_sin);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   m_apm_to_string(buffer, DECIMAL_PLACES, cc_mapm);
   yy_mapm_lib = atof(buffer);

   ydiff = yy_mapm_lib - yy_cos;

   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx,yy_cos);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   xx += xdelta;
   if (xx > LPI)  
     break;
  }

if (pass)
  fprintf(stdout,"... SIN_COS function passes\n");

/*
 *
 */

fprintf(stdout,"Validating the TAN function ... (-PI -> +PI)\n");

xx        = -LPI;
xdelta    = 0.035;
pass      = TRUE;
tolerance = 1.0E-12;

while (1)
  {
   yy = tan(xx);

   m_apm_set_double(aa_mapm, xx);
   m_apm_tan(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);

   if (fabs(yy) > 1.0)
     ydiff = (yy_mapm_lib - yy) / yy;
   else
     ydiff = yy_mapm_lib - yy;

   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %23.16E   C-math-lib = %.15E \n",xx,yy);
      fprintf(stdout,"DIFF = %23.16E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   xx += xdelta;
   if (xx > LPI)  
     break;
  }

if (pass)
  fprintf(stdout,"... TAN function passes\n");
  
/*
 *
 */

fprintf(stdout,"Validating the ARC-SIN function ... (-1 -> +1)\n");

xx        = -0.9999;
xdelta    = 0.035;
pass      = TRUE;
tolerance = 1.0E-14;

while (1)
  {
   yy = asin(xx);

   m_apm_set_double(aa_mapm, xx);
   m_apm_arcsin(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);

   ydiff = yy_mapm_lib - yy;

   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx,yy);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   xx += xdelta;
   if (xx > 0.9999)  
     break;
  }

if (pass)
  fprintf(stdout,"... ARC-SIN function passes\n");
  
/*
 *
 */

fprintf(stdout,"Validating the ARC-COS function ... (-1 -> +1)\n");

xx        = -0.9999;
xdelta    = 0.035;
pass      = TRUE;
tolerance = 1.0E-14;

while (1)
  {
   yy = acos(xx);

   m_apm_set_double(aa_mapm, xx);
   m_apm_arccos(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);

   ydiff = yy_mapm_lib - yy;

   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx,yy);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   xx += xdelta;
   if (xx > 0.9999)  
     break;
  }

if (pass)
  fprintf(stdout,"... ARC-COS function passes\n");
  
/*
 *
 */

fprintf(stdout,
   "Validating the ARC-TAN function ... (-PI/2 < tan(x) < PI/2)\n");

xx        = -1.57;
xdelta    = 0.035;
pass      = TRUE;
tolerance = 1.0E-14;

while (1)
  {
   xx1 = tan(xx);
   sprintf(cbuf,"%.6E",xx1);
   xx1 = atof(cbuf);

   yy = atan(xx1);
   
   m_apm_set_string(aa_mapm, cbuf);
   m_apm_arctan(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);
   
   ydiff = yy_mapm_lib - yy;
   
   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx1,yy);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   xx += xdelta;
   if (xx > 1.57)
     break;
  }

if (pass)
  fprintf(stdout,"... ARC-TAN function passes\n");
  
/*
 *
 */

fprintf(stdout,
   "Spot check the ARC-TAN2 function (calls ARC-TAN) ...\n");

pass      = TRUE;
tolerance = 1.0E-14;

xarray[1] = -1.638;
xarray[2] = -2.179;
xarray[3] = 0.0;
xarray[4] = 0.842;
xarray[5] = 1.976;

yarray[1] = -0.374;
yarray[2] = -4.127;
yarray[3] = 0.0;
yarray[4] = 0.732;
yarray[5] = 3.261;

for (i=1; i <= 5; i++)
  {
   for (j=1; j <= 5; j++)
     {
      if (i != 3 || j != 3)        /* domain error if both == 0.0 */
        {
         xx1 = xarray[i];
         xx2 = yarray[j];
   
         yy = atan2(xx2, xx1);
               
         m_apm_set_double(aa_mapm, xx1);
         m_apm_set_double(bb_mapm, xx2);
         m_apm_arctan2(cc_mapm, DECIMAL_PLACES, bb_mapm, aa_mapm);
         m_apm_to_string(buffer, DECIMAL_PLACES, cc_mapm);
         yy_mapm_lib = atof(buffer);
               
         ydiff = yy_mapm_lib - yy;
               
         if (fabs(ydiff) > tolerance)
           {
            pass = FALSE;
            total_pass = FALSE;
            
            fprintf(stdout,
            "x2 = %10.3E  x1 = %10.3E  C-math-lib = %.15E \n",xx2,xx1,yy);
            fprintf(stdout,
            "DIFF = %17.10E          M_APM_LIB  = %s \n",ydiff,buffer);
           }
        }
     }
  }

if (pass)
  fprintf(stdout,"... ARC-TAN2 function passes\n");


/*
 *
 */

fprintf(stdout,
   "Validating the SIN, COS, TAN functions ... (|input| > PI)\n");

xx        = 2.970;
xdelta    = 0.711;
pass      = TRUE;
tolerance = 1.0E-13;

while (1)
  {
   sprintf(cbuf,"%.6E",xx);
   xx1 = atof(cbuf);

   for (k=0; k <= 1; k++)        /* loop for pos angles, then neg */
     {
      if (k == 1) 
        xx1 = -xx1;

      yy = sin(xx1);
      
      m_apm_set_double(aa_mapm, xx1);
      m_apm_sin(bb_mapm, DECIMAL_PLACES, aa_mapm);
      m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
      yy_mapm_lib = atof(buffer);
      
      ydiff = yy_mapm_lib - yy;
      
      if (fabs(ydiff) > tolerance)
        {
         pass = FALSE;
         total_pass = FALSE;
   
         fprintf(stdout,
	   "SIN: x    = %12.6E   C-math-lib = %.15E \n",xx1,yy);

         fprintf(stdout,
	   "SIN: DIFF = %12.6E   M_APM_LIB  = %s \n",ydiff,buffer);
        }

      yy = cos(xx1);
      
      m_apm_set_double(aa_mapm, xx1);
      m_apm_cos(bb_mapm, DECIMAL_PLACES, aa_mapm);
      m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
      yy_mapm_lib = atof(buffer);
      
      ydiff = yy_mapm_lib - yy;
      
      if (fabs(ydiff) > tolerance)
        {
         pass = FALSE;
         total_pass = FALSE;
   
         fprintf(stdout,
	   "COS: x    = %12.6E   C-math-lib = %.15E \n",xx1,yy);

         fprintf(stdout,
	   "COS: DIFF = %12.6E   M_APM_LIB  = %s \n",ydiff,buffer);
        }

      yy = tan(xx1);
      
      m_apm_set_double(aa_mapm, xx1);
      m_apm_tan(bb_mapm, DECIMAL_PLACES, aa_mapm);
      m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
      yy_mapm_lib = atof(buffer);
      
      if (fabs(yy) > 1.0)
        ydiff = (yy_mapm_lib - yy) / yy;
      else	
        ydiff = yy_mapm_lib - yy;
      
      if (fabs(ydiff) > tolerance)
        {
         pass = FALSE;
         total_pass = FALSE;
   
         fprintf(stdout,
	   "TAN: x    = %12.6E   C-math-lib = %.15E \n",xx1,yy);

         fprintf(stdout,
	   "TAN: DIFF = %12.6E   M_APM_LIB  = %s \n",ydiff,buffer);
        }
     }

   if (xx > 40.0)
     break;

   xx += xdelta;
  }

if (pass)
  fprintf(stdout,"... SIN, COS, TAN (> PI) functions pass\n");


/*
 *
 */

fprintf(stdout,
   "Validating the EXP function ... (-120 -> 120)\n");

xx        = -120.4;
xdelta    = 0.8307;
pass      = TRUE;
tolerance = 1.0E-14;

while (1)
  {
   sprintf(cbuf,"%.6E",xx);
   xx1 = atof(cbuf);

   yy = exp(xx1);
   
   m_apm_set_string(aa_mapm, cbuf);
   m_apm_exp(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);
   
   ydiff = (yy_mapm_lib - yy) / yy;
   
   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx1,yy);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   if (xx > 120.4)
     break;

   xx += xdelta;
  }

if (pass)
  fprintf(stdout,"... EXP function passes\n");
  
  
/*
 *
 */

fprintf(stdout,
   "Validating the LOG function ... (3.1E-7 -> 3.3E+6)\n");

xx        = -15.01;
xdelta    = 0.351;
pass      = TRUE;
tolerance = 1.0E-14;

while (1)
  {
   xx1 = exp(xx);
   sprintf(cbuf,"%.6E",xx1);
   xx1 = atof(cbuf);

   yy = log(xx1);
   
   m_apm_set_string(aa_mapm, cbuf);
   m_apm_log(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);
   
   ydiff = yy_mapm_lib - yy;
   
   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx1,yy);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }

   if (xx > 15.0)
     break;

   xx += xdelta;
  }

if (pass)
  fprintf(stdout,"... LOG function passes\n");
  
/*
 *
 */

fprintf(stdout,
   "Spot check the LOG10 function (calls LOG) ...\n");

pass      = TRUE;
tolerance = 5.0E-14;

xarray[1]  = 1.638;
xarray[2]  = 8.321e-14;
xarray[3]  = 1.0;
xarray[4]  = 7.3231e17;
xarray[5]  = 3.735e-8;
xarray[6]  = 9.997e4;
xarray[7]  = 1.0528e3;
xarray[8]  = 6.319e-6;
xarray[9]  = 2.108e9;
xarray[10] = 8.417e-3;
xarray[11] = 5.347e7;
xarray[12] = 4.291e-23;
xarray[13] = 3.874e28;
xarray[14] = 9.0e2;
xarray[15] = 4.935e-4;

for (i=1; i <= 15; i++)
  {
   xx1 = xarray[i];
   yy  = log10(xx1);
               
   m_apm_set_double(aa_mapm, xx1);
   m_apm_log10(bb_mapm, DECIMAL_PLACES, aa_mapm);
   m_apm_to_string(buffer, DECIMAL_PLACES, bb_mapm);
   yy_mapm_lib = atof(buffer);
               
   ydiff = yy_mapm_lib - yy;
               
   if (fabs(ydiff) > tolerance)
     {
      pass = FALSE;
      total_pass = FALSE;

      fprintf(stdout,"x    = %17.10E   C-math-lib = %.15E \n",xx1,yy);
      fprintf(stdout,"DIFF = %17.10E   M_APM_LIB  = %s \n",ydiff,buffer);
     }
  }

if (pass)
  fprintf(stdout,"... LOG10 function passes\n");
  
/*
 *
 */

fprintf(stdout,
   "Spot check the POW function (calls EXP) ...\n");

pass      = TRUE;
tolerance = 1.0E-14;

xarray[1]  = 1.638;         yarray[1]  = 7.14;
xarray[2]  = 8.321e-12;     yarray[2]  = -5.2;
xarray[3]  = 1.0;           yarray[3]  = 4.19;
xarray[4]  = 7.3231e17;     yarray[4]  = 2.11;
xarray[5]  = 2.0;           yarray[5]  = 32.0;
xarray[6]  = 9.997e4;       yarray[6]  = -1.0;
xarray[7]  = 1.0528e3;      yarray[7]  = -3.83e-2;