/* C---------------------------------------------------------------------- C Program to test CLOG C C Accuracy tests are based on the identity C C log(z) = log(z*z)/2 C and C log(z) = log(17z/16) - log(17/16) 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 CMPLX (or DCMPLX), LOG, REAL (or DBLE), SQRT C C The program as given here must be modified before compiling. C If a single (double) precision version is desired, change all C occurrences of CS (CD) in columns 1 and 2 to blanks. For a C calibration version, change all occurrences of CS and of CC C in columns 1 and 2 to blanks. C C Latest revision - November 20, 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; 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; float complex CF1, CF2, CF3, CV, CX, CX1, CX2, CX3, CY1, CY2, CY3, CZ, CZZ; #endif /* C---------------------------------------------------------------------- C Declare local variables C---------------------------------------------------------------------- */ int IEXP, IOUT, I1, J; #ifdef TDOUBLE double SCALE; #endif #ifdef TSINGLE float SCALE; #endif /* C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- */ /* 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); I1 = IT / 2 + 4; /* SCALE = BETA ** I1; */ #ifdef TDOUBLE SCALE = POW (BETA, (double) I1); #endif #ifdef TSINGLE SCALE = POW (BETA, (float) I1); #endif AX = ONE + ONE; BX = TEN; AY = ZERO; 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 Purify argument and evaluate identities C----------------------------------------------------------------- */ Z = X * SCALE; W = Z + X; X = W - Z; Z = Y * SCALE; W = Z + Y; Y = W - Z; CX = CONVC (X, Y); Z = X * X - Y * Y; W = X * Y; CY = CONVC (Z, W + W); CZ = FLOG (CX); CZZ = FLOG (CY); CZZ = CZZ * HALF; /* C----------------------------------------------------------------- C Accumulate results C----------------------------------------------------------------- */ TABLAT (); XL = XL + DELX; /* 200 CONTINUE */ } /* C----------------------------------------------------------------- C Process and output statistics C----------------------------------------------------------------- */ printf (" Test of LOG(Z) vs LOG(Z*Z)/2\n"); REPORT (IOUT); if (J == 1) { AX = THOUS; BX = AX + AX; AY = -AX; BY = -BX - BX; } else { AX = EPS; BX = HALF * HALF; AY = -AX; BY = -BX; } /* 300 CONTINUE*/ } /* C----------------------------------------------------------------- C Special tests C----------------------------------------------------------------- */ printf (" Special tests\n"); printf (" The identity LOG(Z) = -LOG(1/Z) will be tested\n"); printf (" Z LOG(Z) + LOG(1/Z)\n"); /* DO 320 I = 1, 5*/ for (I = 1; I <= 5; I++) { X = REN (I1) * TEN; Y = REN (I1) * TEN; CX = CONVC (X, Y); CY = ONE / CX; CZ = FLOG (CX) + FLOG (CY); 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 (ONE, ZERO); CZ = FLOG (CX); printf (" LOG(1.0,0.0) = (%15.7E,%15.7E)\n", creal (CZ), cimag (CZ)); printf (" Tests near extreme arguments. These should not\n"); printf (" but may\n trigger error messages\n"); CX = CONVC (XMIN, XMIN); printf (" LOG will be called with the argument "); printf (" (%15.4E,%15.4E)\n", creal (CX), cimag (CX)); CZ = FLOG (CX); printf (" LOG(XMIN,XMIN) = LOG(%15.7E,%15.7E) =\n", creal (CX), cimag (CX)); printf (" (%15.7E,%15.7E)\n", creal (CZ), cimag (CZ)); CX = CONVC (XMAX, XMAX); printf (" LOG will be called with the argument "); printf (" (%15.4E,%15.4E)\n", creal (CX), cimag (CX)); CZ = FLOG (CX); printf (" LOG(XMAX,XMAX) = LOG(%15.7E,%15.7E) =\n", creal (CX), cimag (CX)); printf (" (%15.7E,%15.7E)\n", creal (CZ), cimag (CZ)); /* C----------------------------------------------------------------- C Test of error returns C----------------------------------------------------------------- */ printf (" Test of error returns\n"); CX = CONVC (ZERO, ZERO); printf (" LOG will be called with the argument (%15.4E,%15.4E)\n", creal (CX), cimag (CX)); printf (" This should trigger an error message\n"); CZ = FLOG (CX); printf (" LOG 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 CLOG Test Program ---------- */ #include "tablat.c"