! Last change: DOS 3 Aug 2000 5:51 pm ! *** copyright 2000 *** ! *** filename chex144.f95 *** John F. Monahan ** ! ********************** PROGRAM CHEX144 ! TEST PROGRAM ON RAMAKRISHNA'S CONVOLUTION PROBLEM ! ! USES EXTENDED RANGE COMPLEX ARITHMETIC WHERE NUMBERS ARE ! REPRESENTED BY TRIPLE (X(1), X(2), X(3) ) AND THE VALUE ! IS GIVEN BY ( X(1) + I X(2) )*2**X(3) ! IMPLICIT NONE REAL(KIND=8), DIMENSION(3,16384) :: A REAL(KIND=8), DIMENSION(3) :: AI REAL(KIND=8) R,R2,T,ALR,AL2,ALM REAL(KIND=8) DLGAMA INTEGER I,J,K,M,B,BP1,N,NN ! ! 21 FORMAT(3F20.8) 22 FORMAT(I6,F14.6,4X,3F16.10) 25 FORMAT(' RAMAKRISHNA HASH TABLE/BALL & URN PROBLEM' & & /' POISSON-LIKE COEFFICIENTS, REAL, IMAGINARY AND EXPONENT') 26 FORMAT(//' OVERFLOW PROBABILITES'/' N PR(N,',I2,',',I2, & & ')',11X,'REAL',10X,'IMAGINARY',8X,'EXPONENT') ! OPEN( UNIT=6, FILE='chex144.out' ) ! CASE REPRESENTED IN RAMAKRISHNAN PAPER M = 30 B = 25 N = 10 ! N = 10 GIVES M*B = 750 < (2**10) = 1024 BP1 = B+1 NN = 2 ** N ! AL2 = LOG(2.D0) ! WANT PROBABILITIES FOR N AROUND 400 R = 400.D0/30.D0 ALR = LOG(R) ALM = LOG( REAL(M,8) ) !KIND ! ! INITIALIZE TAIL PART TO ZERO DO I = BP1,NN A(1,I) = 0. A(2,I) = 0. A(3,I) = 0. END DO ! LOOP ON I ! FILL IN THE FRONT PART WITH POISSON PROBS ! WITH RESCALAR R AI(1) = 1. AI(2) = 0. AI(3) = 0. A(1,1) = AI(1) A(2,1) = AI(2) A(3,1) = AI(3) DO I = 2,BP1 AI(1) = AI(1) * R / REAL(I-1,8) !KIND CALL ADJST( AI ) A(1,I) = AI(1) A(2,I) = AI(2) A(3,I) = AI(3) END DO ! LOOP ON I ! WRITE(6,25) DO I = 1,BP1 WRITE(6,21) (A(J,I),J=1,3) END DO ! LOOP ON I ! FORWARD TRANSFORM CALL FFT2NE(A,N,1.D0) ! !WW DO I = 1,BP1 !WW WRITE(6,21) (A(J,I),J=1,3) !WW END DO ! LOOP ON I ! POWER UP BY M DO I = 1,NN R2 = A(1,I)*A(1,I) + A(2,I)*A(2,I) IF( R2 .GT. 0.D0 ) THEN R2 = LOG(R2) * (M/2.D0) ! R2 IS A SQUARE K = R2 / AL2 R2 = R2 - K*AL2 A(3,I) = A(3,I)*M + K R2 = EXP(R2) T = ATAN2( A(2,I), A(1,I) ) T = T * M A(1,I) = R2 * COS(T) A(2,I) = R2 * SIN(T) END IF ! ( R2 .GT. 0.D0 ) END DO ! LOOP ON I ! TRANSFORM BACK CALL FFT2NE(A,N,-1.D0) ! WRITE RESULTS WRITE(6,26) M,B DO I = 1,NN !WW WRITE(6,21) (A(J,I),J=1,3) ! ! COMPUTE LOG OF NORM OF ADJUSTED SERIES ! A(I) * I! / (M*R)**I A(3,I) = A(3,I)*AL2 + (DLGAMA(REAL(I,8)) - (I-1.)*(ALR+ALM)) !KIND ! IF NORM TOO BIG OR TOO SMALL -- OUT OF RANGE -- SKIP IF( (A(3,I).LE.2.D0) .AND. (A(3,I).GE.-20.D0) ) THEN ! ! ATTACH SCALE TO REAL PART TO GET PROBABILITY R = A(1,I)*EXP( A(3,I) ) ! SHIFTED BY ONE LIKE USUAL WITH FFT WRITE(6,22) I-1,R,(A(J,I),J=1,3) END IF ! (BOTH) END DO ! LOOP ON I ! STOP END PROGRAM CHEX144 SUBROUTINE FFT2NE(A,N,SGN) ! *** --- *** MODIFIED FOR EXTENDED RANGE ARITHMETIC *** --- ** ! ! USES EXTENDED RANGE COMPLEX ARITHMETIC WHERE NUMBERS ARE ! REPRESENTED BY TRIPLE (X(1), X(2), X(3) ) AND THE VALUE ! IS GIVEN BY ( X(1) + I X(2) )*2**X(3) ! ! SUBROUTINE FFT2N -- DISCRETE FOURIER TRANSFORM ! ! A COMPLEX VECTOR INPUT: VECTOR TO BE TRANSFORMED ! OUTPUT: TRANSFORMED VECTOR ! N INTEGER INPUT: LOG(BASE 2) OF LENGTH OF VECTOR A ! LENGTH OF A IS 2 ** N ! SGN REAL INPUT: INDICATES FORWARD OR INVERSE TRANSFO ! SGN = 1. THEN FORWARD TRANSFORM ! SGN = -1. THEN INVERSE TRANSFORM ! ! DISCRETE FOURIER TRANSFORM OF A VECTOR A(J) IS GIVEN BY ! NN-1 ! T( A ) = SUM A(J) EXP( - I SGN 2 PI J K / NN ) ! J K=0 ! ! WHERE I = SQRT( -1 ) ! NN = 2 ** N ! PI = 3.1415923565... ! ! NOTE THAT THE INDEXING OF A IS SHIFTED BY ONE -- ! A(0) STORED IN A(1), ... , A( NN-1 ) STORED IN A( NN ) ! ! J F MONAHAN (NOV,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695-8203 ! RECODED DECEMBER, 1999 FOR FORTRAN 95 ! ! THIS IS A 'POWER OF 2' ALGORITHM BASED ON THE PRESENTATION IN ! ! A. V. AHO, J. E. HOPCROFT, AND J. D. ULLMAN (1974) THE DESIGN AND ! ANALYSIS OF COMPUTER ALGORITHMS, ADDISON-WESLEY, READING. ! ! ! BEGIN WITH DECLARATIONS IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION(3,N), INTENT(IN OUT) :: A REAL(KIND=8), INTENT(IN) :: SGN REAL(KIND=8) C,S,ALEFT,ARIGHT COMPLEX(KIND=16), DIMENSION(25) :: WK COMPLEX(KIND=16) WS LOGICAL CARRY,NEW LOGICAL, DIMENSION(25) :: RBITS INTEGER I,J,L,LM,LL,M,MM,M2,NN ! ! COMPUTE WK(I) = EXP( - 2 PI IMAG SGN / 2**I ) ! ! SPECIAL VALUES FIRST WK(1) = CMPLX( -1., 0. ,16) !KIND WK(2) = CMPLX( 0., -SGN ,16) !KIND ! SKIP IF N IS ONLY 1 OR 2 IF( N .GT. 2 ) THEN ! AVOID SINES AND COSINES -- JUST USE SQUARE ROOTS C = 0. S = 1. DO I = 3,N C = SQRT( ( 1. + C ) / 2. ) S = S / ( 2.*C ) WK(I) = CMPLX( C, -SGN*S ,16) !KIND END DO ! LOOP ON I END IF ! ( N .GT. 2 ) ! ! THIS IS THE MAIN LOOP OF THE ALGORITHM ! DO MM = 1,N ! CONVERT TO REVERSE LOOP M = N - MM ! M STEPS DOWN FROM N-1 TO 0 ! ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I ! M2 M2 = 2**M ! LM CONTROLS THE NUMBER OF SUBLOOPS LM = 2**(MM-1) ! DO LL = 1,LM ! L IS THE MAJOR INDEX VARIABLE L = 2*M2*(LL-1) ! INITIALIZE TO GET W ** REV( L / M2 ) WS = CMPLX(1.D0, 0.D0, 16) !KIND CARRY = .TRUE. ! START AT 2 (SHIFT) SINCE L/M2 IS ALWAYS EVEN DO I = 2,N NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) WS = WS *WK(I) END DO ! LOOP ON I ! WS IS NOW W ** REV( L / M2 ) ! NOW DO THE REAL COMPUTATION DO J = 1,M2 CALL COMBO(L+J,L+J+M2,WS) ! SPECIAL ROUTINE FOR EXTENDED END DO ! LOOP ON J ARITHMETIC ! DONE WITH LOOP ON L END DO ! LOOP ON LL ! DONE WITH MAIN LOOP ON M END DO ! LOOP ON MM ! ! SCRAMBLE THE NUMBERS -- SWITCH A(J) AND A( REV(J) ) ! BUT DON'T DO IT TWICE ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I NN = 2 ** N ! LOOP THROUGH ALL THE NUMBERS DO LL = 1,NN ! DO THE BIT REVERSAL LM = 0 M2 = NN CARRY = .TRUE. ! START AT 1 HERE -- USE ALL THE NUMBERS DO I = 1,N M2 = M2 / 2 NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) LM = LM + M2 END DO ! LOOP ON I ! DON'T SWITCH TWICE IF( LL-1 .LT. LM ) THEN DO J = 1,3 ! A IS (3,N) MATRIX ALEFT = A(J,LL) ARIGHT = A(J,LM+1) A(J,LL) = ARIGHT A(J,LM+1) = ALEFT END DO ! LOOP ON J END IF ! ( LL-1 .LT. LM ) END DO ! LOOP ON LL ! IF INVERSE TRANSFORM, REMEMBER TO DIVIDE BY NN IF( SGN .GT. 0. ) RETURN ! C = REAL(NN) DO I = 1,NN ! A(I) = A(I) / C A(3,I) = A(3,I) - REAL(N) END DO ! LOOP ON I RETURN CONTAINS SUBROUTINE COMBO(I,J,WS) ! MIMICS IN EXTENDED ARITHMETIC THE BASIC STEP IN THE FFT ! A(I) = A(I) + WS*A(J) ! A(J) = A(I) - WS*A(J) ! ! USES EXTENDED RANGE COMPLEX ARITHMETIC WHERE NUMBERS ARE ! REPRESENTED BY TRIPLE (X(1), X(2), X(3) ) AND THE VALUE ! IS GIVEN BY ( X(1) + I X(2) )*2**X(3) ! IMPLICIT NONE INTEGER, INTENT(IN) :: I,J COMPLEX(KIND=16), INTENT(IN) :: WS REAL(KIND=8), DIMENSION(3) :: AI REAL(KIND=8) AMEX,DI,DJ INTEGER K ! COMPLEX(KIND=16) AIC,AJC,APC,AMC ! AMEX = MAX( A(3,I), A(3,J) ) AIC = CMPLX( A(1,I), A(2,I) ,16) !KIND AJC = CMPLX( A(1,J), A(2,J) ,16) !KIND ! DI = 0. IF( (A(3,I)-AMEX) .GT. -60 ) DI = 2. ** (A(3,I)-AMEX) DJ = 0. IF( (A(3,J)-AMEX) .GT. -60 ) DJ = 2. ** (A(3,J)-AMEX) ! APC = DI*AIC + WS*DJ*AJC AMC = DI*AIC - WS*DJ*AJC ! AI(1) = REAL( APC ,8) !KIND AI(2) = QIMAG( APC ) !KIND AI(3) = AMEX CALL ADJST(AI) DO K = 1,3 A(K,I) = AI(K) END DO AI(1) = REAL( AMC ,8) !KIND AI(2) = QIMAG( AMC ) !KIND AI(3) = AMEX CALL ADJST(AI) DO K = 1,3 A(K,J) = AI(K) END DO RETURN END SUBROUTINE COMBO END SUBROUTINE FFT2NE SUBROUTINE ADJST(AI) ! ! ROUTINE TO NORMALALIZE EXTENDED RANGE COMPLEX ARITHMETIC ! WHERE THE NUMBER IS REPRESENTED BY ( X(1) + I X(2) )*2**X(3) ! ! KEEPS THE COMPLEX NORM OF THE NUMBER LESS THAN TWO ! NOTE THAT ZERO IS REPRESENTED BY X(3) = -1600000 ! IMPLICIT NONE REAL(KIND=8), DIMENSION(3), INTENT(IN OUT) :: AI REAL(KIND=8) AM ! AM = MAX( ABS(AI(1)), ABS(AI(2)) ) IF( AM .EQ. 0. ) THEN AI(3) = -1600000. AI(1) = 0. AI(2) = 0. RETURN END IF ! ( AM .EQ. 0. ) DO WHILE( AM .LT. 1.0 ) AI(3) = AI(3) - 1. AM = AM * 2. AI(1) = AI(1) * 2. AI(2) = AI(2) * 2. END DO ! WHILE ! DO WHILE( AM .GT. 2.0 ) AI(3) = AI(3) + 1. AM = AM / 2. AI(1) = AI(1) / 2. AI(2) = AI(2) / 2. END DO ! WHILE RETURN END SUBROUTINE ADJST REAL(KIND=8) FUNCTION DLGAMA(X) ! COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: X REAL(KIND=8), DIMENSION(7) :: P ! COEFS FOR RATIONAL FUNCTION REAL(KIND=8), DIMENSION(6) :: Q REAL(KIND=8), DIMENSION(3) :: PT ! COEFS FOR STIRLING CORRECTION REAL(KIND=8), DIMENSION(2) :: QT REAL(KIND=8) XN,A,V REAL(KIND=8), PARAMETER :: LSQR2PI = .918938533204672742D0 DATA P/ .378601050348257245d4, .207745979389418732d4, & & .893581804523749814d3, .222112396168011795d3, & & .489543462279099381d2, .61260674503360843d1, & & .778079585613300576d0 / DATA Q/ .378601050348257187d4, .476793860503687915d3, & & -.867230987531102994d3, .835500586679197696d2, & & .507884753288954097d2, -.134004147857813483d2 / DATA PT/ .288119283935546015d0, .498030766924499634d0, & & .691561607375687d-1 / DATA QT/ .345743140722674507d1, .609161691641660296d1 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0.D0 ) THEN DLGAMA = -1.D0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8.D0 ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = (PT(1)+(PT(2)+PT(3)/XN)/XN) / (QT(1)+(QT(2)+1.D0/XN)/XN) DLGAMA = (X-0.5D0)*LOG(X) - X + V/X + LSQR2PI ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.D0 XN = X DO WHILE( XN .LT. 2.D0 ) V = V/XN XN = XN + 1.D0 END DO ! WHILE( XN .LT. 2.D0 ) DO WHILE( XN .GT. 3.D0 ) XN = XN - 1.D0 V = V * XN END DO ! WHILE( XN .GT. 3.D0 ) A = XN-2.D0 ! RATIONAL FUNCTION APPROXIMATION V =(P(1)+A*(P(2)+A*(P(3)+A*(P(4)+A*(P(5)+A*(P(6)+A*P(7))))))) & & * V / (Q(1)+A*(Q(2)+A*(Q(3)+A*(Q(4)+A*(Q(5)+A*(Q(6)+A)))))) DLGAMA = LOG(V) ENDIF RETURN END FUNCTION DLGAMA ! *** end of filename chex144.f95 *********************