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---------------------------------------------------------------------- COMMON /ELEFUN/ CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3, 1 CZ,CZZ,AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN, 2 IBETA,IT,KX1,KX2,KX3,KY1,KY2,KY3,K1,K2,K3,N, 3 IFLAG,RFLAG,PFLAG LOGICAL IFLAG,RFLAG,PFLAG INTEGER IBETA,IT,KX1,KX2,KX3,KY1,KY2,KY3,K1,K2,K3,N CS REAL DOUBLE PRECISION 1 AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN CS COMPLEX COMPLEX*16 1 CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,CZ,CZZ C---------------------------------------------------------------------- C Declare local variables C---------------------------------------------------------------------- INTEGER I,IEXP,IOUT,IRND,I1,J,MACHEP,MAXEXP,MINEXP,NDUM, 1 NEGEP,NGRD CS REAL DOUBLE PRECISION 1 BETA,CONV,DELX,DELY,EPS,EPSNEG,HALF,ONE,REN,SCALE,TEN, 2 THOUS,W,X,XL,XMAX,XMIN,Y,Z,ZERO CS COMPLEX COMPLEX*16 1 CONVC,CY CC 2 ,FLOG C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- DATA IOUT/6/ CS DATA ZERO,HALF,ONE,TEN,THOUS/0.0E0,0.5E0,1.0E0,10.0E0,1.0E3/ DATA ZERO,HALF,ONE,TEN,THOUS/0.0D0,0.5D0,1.0D0,10.0D0,1.0D3/ C---------------------------------------------------------------------- C Statement functions for conversion between integer or complex C and float C---------------------------------------------------------------------- CS CONV(NDUM) = REAL(NDUM) CS CONVC(X,Y) = CMPLX(X,Y) CONV(NDUM) = DBLE(NDUM) CONVC(X,Y) = DCMPLX(X,Y) C---------------------------------------------------------------------- CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) BETA = CONV(IBETA) ALBETA = LOG(BETA) AIT = CONV(IT) I1 = IT/2 + 4 SCALE = BETA ** I1 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 CALL RESET DELX = (BX - AX) / XN DELY = (BY - AY) XL = AX DO 200 I = 1, N 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 = LOG(CX) CZZ = LOG(CY) CC CZ = FLOG(CX) CC CZZ = FLOG(CY) CZZ = CZZ * HALF C----------------------------------------------------------------- C Accumulate results C----------------------------------------------------------------- CALL TABLAT XL = XL + DELX 200 CONTINUE C----------------------------------------------------------------- C Process and output statistics C----------------------------------------------------------------- WRITE (IOUT,1000) CALL REPORT(IOUT) IF (J .EQ. 1) THEN AX = THOUS BX = AX + AX AY = -AX BY = -BX - BX ELSE AX = EPS BX = HALF*HALF AY = -AX BY = -BX END IF 300 CONTINUE C----------------------------------------------------------------- C Special tests C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) DO 320 I = 1, 5 X = REN(I1) * TEN Y = REN(I1) * TEN CX = CONVC(X,Y) CY = ONE / CX CZ = LOG(CX) + LOG(CY) WRITE (IOUT,1031) CX,CZ 320 CONTINUE WRITE (IOUT,1040) CX = CONVC(ONE,ZERO) CZ = LOG(CX) WRITE (IOUT,1041) CZ WRITE (IOUT,1050) CX = CONVC(XMIN,XMIN) WRITE (IOUT,1061) CX CZ = LOG(CX) WRITE (IOUT,1051) CX,CZ CX = CONVC(XMAX,XMAX) WRITE (IOUT,1061) CX CZ = LOG(CX) WRITE (IOUT,1052) CX,CZ C----------------------------------------------------------------- C Test of error returns C----------------------------------------------------------------- WRITE (IOUT,1060) CX = CONVC(ZERO,ZERO) WRITE (IOUT,1062) CX CZ = LOG(CX) WRITE (IOUT,1065) CZ WRITE (IOUT,1100) STOP C--------------------------------------------------------------------- 1000 FORMAT(' Test of LOG(Z) vs LOG(Z*Z)/2'//) 1025 FORMAT(' Special tests'//) 1030 FORMAT(' The identity LOG(Z) = -LOG(1/Z) will be tested.'// 1 15X,'Z',22X,'LOG(Z) + LOG(1/Z)'/) 1031 FORMAT(2(E15.7,E15.7)/) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' LOG(1.0,0.0) = (',(E15.7,E15.7),')'//) 1050 FORMAT(//' Tests near extreme arguments. These should not,', 1 ' but may'/' trigger error messages'//) 1051 FORMAT(' LOG(XMIN,XMIN) = LOG(',(E15.7,E15.7),') ='/40X,'(', 1 (E15.7,E15.7),')'/) 1052 FORMAT(' LOG(XMAX,XMAX) = LOG(',(E15.7,E15.7),') ='/40X,'(', 1 (E15.7,E15.7),')'/) 1060 FORMAT(' Test of error returns'/) 1061 FORMAT(/' LOG will be called with the argument (', 1 (E15.4,E15.4),')'//) 1062 FORMAT(/' LOG will be called with the argument (', 1 (E15.4,E15.4),')'/ 2 ' This should trigger an error message'//) 1065 FORMAT(' LOG returned the value (',(E15.4,E15.4),')'//) 1100 FORMAT(' This concludes the tests') C---------- Last card of CLOG Test Program ---------- END CC COMPLEX FUNCTION FLOG(CX) C---------------------------------------------------------------------- C Calibration function for clog tests. C C Latest revision - October 31, 1990 C C Author - W. J. Cody C Argonne National Laboratory C C---------------------------------------------------------------------- CC COMPLEX CX CC COMPLEX*16 DCX C---------------------------------------------------------------------- CC DCX = CX CC FLOG = LOG(DCX) CC RETURN C---------- Last line of FLOG ---------- CC END