Quelle  tcexp.c

/*
C----------------------------------------------------------------------
C  Program to test CEXP
C
C  Accuracy tests are based on the identity
C
C            exp(z) = exp(z-w)*exp(w)
C
C  Data required:
C
C     None
C
C  Subprograms required from this package:
C
C     MACHAR - An environmental inquiry program providing
C              information on the floating point arithmetic
C              system.  Note that the call to MACHAR can
C              be deleted provided the following four
C              parameters are assigned the values indicated:
C
C              IBETA  - the radix of the floating point system;
C              IT     - the number of base-IBETA digits in the
C                       significand of a floating-point number;
C              XMIN   - the smallest non-vanishing floating-point
C                       power of the radix;
C              XMAX   - the largest finite floating-point no.
C
C     REN(K) - a function subprogram returning random real
C              numbers uniformly distributed over (0,1)
C
C     RESET, TABLAT, REPORT - programs to report results.
C
C
C  Generic Fortran subprograms required:
C
C         AIMAG (or DIMAG), AINT, CMPLX (or DCMPLX), EXP,
C
C         LOG, REAL (or DBLE), SQRT
C
C
C  Latest revision - October 31, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
C  Declare named common
C----------------------------------------------------------------------
*/

#include "celefunt.h"

int IFLAG, RFLAG, PFLAG;
int KX1, KX2, KX3, KY1, KY2, KY3;

#ifdef TDOUBLE
double AX, AY, BX, BY, RX6, RX7, RY6, RY7, R6, R7, XN;
double complex CF1, CF2, CF3, CV, CX, CX1, CX2, CX3, CY1, CY2, CY3, CZ,
  CZZ;
#endif

#ifdef TSINGLE
float AX, AY, BX, BY, RX6, RX7, RY6, RY7, R6, R7, XN;
float complex CF1, CF2, CF3, CV, CX, CX1, CX2, CX3, CY1, CY2, CY3, CZ,
  CZZ;
#endif

/*
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
*/


/*
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
*/
#ifdef TDOUBLE
double complex EXPDEL = 6.2416044877018563681e-2
+ 6.6487597751003112768e-2 * _Imaginary_I;
double complex CDEL = 0.0625 + 0.0625 * _Imaginary_I;
#endif

#ifdef TSINGLE
float complex EXPDEL = 6.2416044877018563681e-2f
+ 6.6487597751003112768e-2f * _Imaginary_I;
float complex CDEL = 0.0625f + 0.0625f * _Imaginary_I;
#endif
/*
C----------------------------------------------------------------------
C  Statement functions for conversion between integer or complex
C  and float
C----------------------------------------------------------------------
*/
/*
C----------------------------------------------------------------------
*/
int
main ()
{
  /*
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
     */
  IBETA = 2;
  IT = IT_VAL;
  XMIN = XMIN_VAL;
  XMAX = XMAX_VAL;
  BETA = CONV (IBETA);
  ALBETA = LOG (BETA);
  AIT = CONV (IT);
  AX = ZERO;
  BX = ONE;
  AY = ONE16;
  BY = BX;
  N = 2000;
  XN = CONV (N);
  I1 = 0;
  RFLAG = TRUE;
  IFLAG = TRUE;
  PFLAG = FALSE;
  /*
C---------------------------------------------------------------------
C  Random argument accuracy tests
C---------------------------------------------------------------------
*/
  /* DO 300 J = 1, 3 */
  for (J = 1; J <= 3; J++)
    {
      RESET ();
      DELX = (BX - AX) / XN;
      DELY = (BY - AY);
      XL = AX;
      /* DO 200 I = 1, N */
      for (I = 1; I <= N; I++)
 {
   X = DELX * REN (I1) + XL;
   Y = DELY * REN (I1) + AY;
   /*
C-----------------------------------------------------------------
C  Carefully purify arguments and evaluate identities
C-----------------------------------------------------------------
*/
   CX = CONVC (X, Y);
   CY = CX - CDEL;
   CZ = FEXP (CX);
   CZZ = FEXP (CY);
   CZZ = CZZ + CZZ * EXPDEL;
   /*
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
*/
   TABLAT ();
   XL = XL + DELX;
 }
      /*  200    CONTINUE */
      /*
C-----------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------
*/
      printf (" Test of EXP(Z) vs EXP(Z-DEL)EXP(DEL), DEL = \n");
      printf ("(1/16,1/16)\n");

      REPORT (IOUT);
      if (J == 1)
 {
   AX = ONEP6;
   BX = THREE;
 }
      else
 {
   AX = SIXTEN;
   BX = AX + ONE;
 }
      AY = AX;
      BY = BX;
      /*  300 CONTINUE */
    }
  /*
C---------------------------------------------------------------------
C  Special tests
C---------------------------------------------------------------------
*/
  printf (" Special tests\n");
  printf (" The identity  EXP(Z)*EXP(-Z) = 1.0  will be tested.\n");
  printf ("        Z                        F(Z)*F(-Z) - 1\n");

  /* DO 320 I = 1, 5 */
  for (I = 1; I <= 5; I++)
    {
      X = REN (I1) * BETA;
      Y = REN (I1) * BETA;
      CX = CONVC (X, Y);
      CY = -CX;
      CZ = FEXP (CX) * FEXP (CY) - ONE;
      printf ("(%15.7E,%15.7E) (%15.7E,%15.7E)\n",
       creal (CX), cimag (CX), creal (CZ), cimag (CZ));
      /*  320 CONTINUE*/
    }
  printf (" Test of special arguments\n");
  CX = CONVC (ZERO, ZERO);
  CZ = FEXP (CX) - ONE;
  printf (" EXP(0.0,0.0) - 1.0E0 = (%15.7E,%15.7E)\n",
   creal (CZ), cimag (CZ));
  printf (" Test cos(y), sin(y), and exp(x)\n");
  CX = CONVC (ZERO, TEN);
  CZ = FEXP (CX);
  X = COS (TEN);
  Y = SIN (TEN);
  W = (CONVR (CZ) - X) / X;
  printf (" [REAL(EXP(0,10)) - COS(10)]/COS(10) = %15.7E\n", W);

  W = (CONVI (CZ) - Y) / Y;
  printf (" [IMAG(EXP(0,10)) - SIN(10)]/SIN(10) = %15.7E\n", W);
  CX = CONVC (TEN, ZERO);
  CZ = FEXP (CX);
  X = EXP (TEN);
  W = (CONVR (CZ) - X) / X;
  printf (" [REAL(EXP(10,0)) - EXP(10)]/EXP(10) = %15.7E\n", W);
  printf (" Tests near extreme arguments\n");
  printf (" The real part of the first result below may underflow.\n");
  X = AINT (LOG (XMIN));
  CX = CONVC (X, ONE);
  CZ = FEXP (CX);
  printf (" EXP(%13.6E,%13.6E) = (%13.6E,%13.6E)\n",
   creal (CX), cimag (CX), creal (CZ), cimag (CZ));
  Y = SQRT (XMAX);
  X = AINT (LOG (XMAX));
  CX = CONVC (X, ONE);
  CZ = FEXP (CX);
  printf (" EXP(%13.6E,%13.6E) = (%13.6E,%13.6E)\n",
   creal (CX), cimag (CX), creal (CZ), cimag (CZ));
  CX = CONVC (X * HALF, TEN);
  CY = CX * HALF;
  CZ = FEXP (CX);
  CZZ = FEXP (CY);
  CZZ = CZZ * CZZ;
  printf (" If EXP(%13.6E,%13.6E) = (%13.6E,%13.6E)\n",
   creal (CX), cimag (CX), creal (CZ), cimag (CZ));
  printf (" is not about\n EXP(%13.6E,%13.6E)**2 = \n",
   creal (CY), cimag (CY));
  printf (" (%13.6E,%13.6E) there is an\n argument reduction error\n",
   creal (CZZ), cimag (CZZ));

  /*
C---------------------------------------------------------------------
C  Test of error returns
C---------------------------------------------------------------------
*/
  printf (" Test of error returns\n");
  X = -ONE / SQRT (XMIN);
  CX = CONVC (X, ONE);
  printf (" EXP will be called with the argument (%15.4E,%15.4E)\n",
   creal (CX), cimag (CX));
  printf (" The real component of the argument is so negative\n");
  printf (" that the real exponential should underflow to zero.\n");
  printf (" Thus, the result should be (0,0).\n");
  CZ = FEXP (CX);
  printf (" EXP returned the value (%15.4E,%15.4E)\n",
   creal (CZ), cimag (CZ));
  Y = SQRT (XMAX);
  CX = CONVC (TEN, Y);
  printf (" EXP will be called with the argument (%15.4E,%15.4E)\n",
   creal (CX), cimag (CX));
  printf (" The imaginary component of the argument is too large\n");
  printf (" for the sin and cos computations to make sense.\n");
  CZ = FEXP (CX);
  printf (" EXP returned the value (%15.4E,%15.4E)\n",
   creal (CZ), cimag (CZ));
  CX = CONVC (-X, -ONE);
  printf (" EXP will be called with the argument (%15.4E,%15.4E)\n",
   creal (CX), cimag (CX));
  printf (" The real component of the argument is so large\n");
  printf (" that the real exponential should overflow.  The\n");
  printf (" result should be (inf,inf) with an error return.\n");
  CZ = FEXP (CX);
  printf (" EXP returned the value (%15.4E,%15.4E)\n",
   creal (CZ), cimag (CZ));
  printf (" This concludes the tests\n");
  exit (0);
}

/*
C---------------------------------------------------------------------
*/
/*
C---------- Last card of CEXP Test Program ----------
*/

#include "tablat.c"