* * $Id: b100m.F,v 1.1.1.1 1996/04/01 15:01:12 mclareni Exp $ * * $Log: b100m.F,v $ * Revision 1.1.1.1 1996/04/01 15:01:12 mclareni * Mathlib gen * * #include "gen/pilot.h" SUBROUTINE B100M C Routine to test MATHLIB routines KBINOM, BINOM and DBINOM (B100) C Note KBINOM is only called for integer values of the parameter #include "imp64r.inc" REAL BINOM C Set the total number of tests PARAMETER ( NT=54 ) C Specify the largest error allowed for a successful test PARAMETER ( TSTERR=1D-15 ) PARAMETER (RTSTERR=1D-6 ) INTEGER KBINOM DIMENSION KEX1(NT),EX2(NT),REX3(NT),EX4(NT),REX5(NT) PARAMETER (Z1 = 1, Z2 = 2, HALF = Z1/2) #include "iorc.inc" C Analytical values expected to be obtained DATA (KEX1(J),EX2(J),REX3(J),EX4(J),REX5(J),J=1,10) 1 / 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, 2 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, 3 -8, -8.0, -8.0,-0.750000000000000000D+01,-0.750000000E+01, 4 36, 36.0, 36.0, 0.318750000000000000D+02, 0.318750000E+02, 5 -120,-120.0,-120.0,-0.100937499999999993D+03,-0.100937500E+03, 6 330, 330.0, 330.0, 0.264960937499999943D+03, 0.264960937E+03, 7 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, 8 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, 9 -6, -6.0, -6.0,-0.550000000000000000D+01,-0.550000000E+01, + 21, 21.0, 21.0, 0.178750000000000000D+02, 0.178750000E+02/ DATA (KEX1(J),EX2(J),REX3(J),EX4(J),REX5(J),J=11,20) + /-56, -56.0, -56.0,-0.446875000000000000D+02,-0.446875000E+02, + 126, 126.0, 126.0, 0.949609375000000000D+02, 0.949609375E+02, + 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, + 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, + -4, -4.0, -4.0,-0.350000000000000000D+01,-0.350000000E+01, + 10, 10.0, 10.0, 0.787500000000000000D+01, 0.787500000E+01, + -20, -20.0, -20.0,-0.144374999999999993D+02,-0.144375000E+02, + 35, 35.0, 35.0, 0.234609374999999964D+02, 0.234609375E+02, + 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, + 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01/ DATA (KEX1(J),EX2(J),REX3(J),EX4(J),REX5(J),J=21,30) 1 /-2, -2.0, -2.0,-0.150000000000000000D+01,-0.150000000E+01, 2 3, 3.0, 3.0, 0.187500000000000000D+01, 0.187500000E+01, 3 -4, -4.0, -4.0,-0.218749999999999956D+01,-0.218750000E+01, 4 5, 5.0, 5.0, 0.246093749999999933D+01, 0.246093750E+01, 5 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, 6 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, 7 0, 0.0, 0.0, 0.500000000000000000D+00, 0.500000000E+00, 8 0, 0.0, 0.0,-0.125000000000000000D+00,-0.125000000E+00, 9 0, 0.0, 0.0, 0.625000000000000000D-01, 0.625000000E-01, + 0, 0.0, 0.0,-0.390625000000000000D-01,-0.390625000E-01/ DATA (KEX1(J),EX2(J),REX3(J),EX4(J),REX5(J),J=31,40) 1 / 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, 2 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, 3 2, 2.0, 2.0, 0.250000000000000000D+01, 0.250000000E+01, 4 1, 1.0, 1.0, 0.187500000000000000D+01, 0.187500000E+01, 5 0, 0.0, 0.0, 0.312499999999999972D+00, 0.312500000E+00, 6 0, 0.0, 0.0,-0.390624999999999965D-01,-0.390625000E-01, 7 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, 8 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, 9 4, 4.0, 4.0, 0.450000000000000000D+01, 0.450000000E+01, + 6, 6.0, 6.0, 0.787500000000000000D+01, 0.787500000E+01/ DATA (KEX1(J),EX2(J),REX3(J),EX4(J),REX5(J),J=41,54) + / 4, 4.0, 4.0, 0.656249999999999978D+01, 0.656250000E+01, + 1, 1.0, 1.0, 0.246093749999999978D+01, 0.246093750E+01, + 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, + 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, + 6, 6.0, 6.0, 0.650000000000000000D+01, 0.650000000E+01, + 15, 15.0, 15.0, 0.178750000000000000D+02, 0.178750000E+02, + 20, 20.0, 20.0, 0.268125000000000000D+02, 0.268125000E+02, + 15, 15.0, 15.0, 0.234609375000000000D+02, 0.234609375E+02, + 0, 0.0, 0.0, 0.000000000000000000D+00, 0.000000000E+00, + 1, 1.0, 1.0, 0.100000000000000000D+01, 0.100000000E+01, + 8, 8.0, 8.0, 0.850000000000000000D+01, 0.850000000E+01, + 28, 28.0, 28.0, 0.318750000000000000D+02, 0.318750000E+02, + 56, 56.0, 56.0, 0.690624999999999929D+02, 0.690625000E+02, + 70, 70.0, 70.0, 0.949609374999999893D+02, 0.949609375E+02/ CALL HEADER('B100',0) KERMAX=0 ERMAX=0D0 RERMAX=0E0 WRITE(LOUT,'(/15X,''ERRORS OBTAINED WHEN TESTING:'')') WRITE(LOUT,'(/7X,''KBINOM (X,K) = KB'')') WRITE(LOUT,'(/7X,''BINOM (X,K) = B BINOM (X+1/2,K) = BH'')') #if defined(CERNLIB_DOUBLE) WRITE(LOUT,'(/7X,''DBINOM (X,K) = DB DBINOM (X+1/2,K)= DBH'')') #endif WRITE(LOUT,'(//7X,'' for X = -8,-6,-4,-2,0,2,4,6,8'', + '' each with K from -1 to 4'')') WRITE(LOUT,'(/7X, #if defined(CERNLIB_DOUBLE) + ''KB'',6X,''DB'',8X,''B'',23X,''DBH'',13X,''BH'')') #endif #if !defined(CERNLIB_DOUBLE) + ''KB'',6X,''B'',13X,''BH'')') #endif J= 1 DO 1 N = -8,8,2 DO 1 K = -1,4 KB=KBINOM(N,K) IF (KEX1(J) .NE. 0) THEN KER1=ABS((KB-KEX1(J))/KEX1(J)) ELSE KER1=ABS(KB-KEX1(J)) ENDIF KERMAX=MAX(KERMAX,KER1) DX=N #if defined(CERNLIB_DOUBLE) D1=DBINOM(DX,K) IF (EX2(J) .NE. 0) THEN ER2=ABS((D1-EX2(J))/EX2(J)) ELSE ER2=ABS(D1-EX2(J)) ENDIF ERMAX=MAX(ERMAX,ER2) #endif RDX=DX RD1=BINOM(RDX,K) IF (REX3(J) .NE. 0) THEN RER3=ABS((RD1-REX3(J))/REX3(J)) ELSE RER3=ABS(RD1-REX3(J)) ENDIF RERMAX=MAX(RERMAX,RER3) DX=N+HALF #if defined(CERNLIB_DOUBLE) D2=DBINOM(DX,K) IF (EX2(J) .NE. 0) THEN ER4=ABS((D2-EX4(J))/EX4(J)) ELSE ER4=ABS(D2-EX4(J)) ENDIF ERMAX=MAX(ERMAX,ER4) #endif RDX=DX RD2=BINOM(RDX,K) IF (REX5(J) .NE. 0) THEN RER3=ABS((RD2-REX5(J))/REX5(J)) ELSE RER5=ABS(RD2-REX5(J)) ENDIF RERMAX=MAX(RERMAX,RER5) J=J+1 #if !defined(CERNLIB_DOUBLE) 1 WRITE(LOUT,'(5X,I4,D10.1,E16.1)') KER1,RER3,RER5 #endif #if defined(CERNLIB_DOUBLE) 1 WRITE(LOUT,'(5X,I4,D10.1,E10.1,D25.1,E16.1)') + KER1,ER2,RER3,ER4,RER5 #endif DX=Z2**32 RDX=DX RZ2=Z2 RD1=BINOM(RZ2**32,2) RTEST=RDX*(RDX-1)/2 IF (RTEST .NE. 0E0)RER7=ABS((RD1-RTEST)/RTEST) RERMAX=MAX(RERMAX,RER7) #if defined(CERNLIB_DOUBLE) D1=DBINOM(Z2**32,2) TEST=DX*(DX-1)/2 IF (TEST .NE. 0D0) ER6=ABS((D1-TEST)/TEST) ERMAX=MAX(ERMAX,ER6) WRITE(LOUT,'(/7X,''DBINOM(2**32,2)'',10X,''EXACT VALUE'')') WRITE(LOUT,'(3X,D25.15,D25.15)') D1,TEST #endif WRITE(LOUT,'(/7X,''BINOM(2**32,2)'',6X,''EXACT VALUE'')') WRITE(LOUT,'(4X,E18.9,E18.9)') RD1,RTEST WRITE(LOUT,'(/7X,''Largest Error for KBINOM was'',I4)')KERMAX WRITE(LOUT,'(/7X,''Largest Error for BINOM was'',E10.1)')RERMAX RKERMAX=KERMAX RERMAX= MAX(RERMAX,RKERMAX) #if defined(CERNLIB_DOUBLE) WRITE(LOUT,'(/7X,''Largest Error for DBINOM was'',D10.1)')ERMAX RERMAX= MAX(RERMAX,SNGL(ERMAX)) #endif WRITE(LOUT,'(/7X,''TESTING ERROR MESSAGES:''/)') N=KBINOM(2**30,4) IRC=ITEST('B100',RERMAX .LE. RTSTERR) CALL PAGEND('B100') RETURN END