! Last change: DOS 3 Aug 2000 6:27 pm ! *** copyright 2000 *** ! *** filename tivmm1.f95 *** John F. Monahan ** ! ********************** program tivmm1 ! compute the distribution of the sum of reciprocal uniforms ! implicit none real(kind=8), dimension(3,65536) :: p real(kind=8), dimension(65536,3) :: pj real(kind=8), dimension(3) :: ppi real(kind=8) delta,xi,fln,pi,pold,pnew,psum,al2,rr,t integer i,j,k,n,np,lgnp,ms ! 21 format(' Case',i2,' Sum of',i2,' variables, using',i8,' points',& & /4x,'sum of probabilities',f9.6 & & /4x,'First value: real, imaginary, and exponent'/9x,3f16.9) 22 format(i8,5f12.8) ! ! open( unit=6, file='tivmm1.out' ) open( unit=8, file='tivmm1.dat' ) ! ! pound on it -- np = number of points lgnp = 16 np = 2**lgnp ! delta = 2**(-8), range=(0., 256=2**8) ms = 256 delta = 1.d0/real(ms,8) !kind al2 = log(2.d0) ! ! outer loop -- several values of j do j = 1,3 ! n = sample sizes -- 8, 16, 32 n = 2**(j+2) fln = real(n,8) !kind ! pold = delta pi = 0. psum = 0. ! now for many values of x do i = 1,np pi = 0.d0 xi = real(i+ms,8)*delta !kind ! delta*density of ( 1/U - 1 ) pnew = delta / ( xi * xi ) ! trapezoid rule pi = (pold + pnew) / 2.d0 psum = psum + pi pold = pnew ! ppi(1) = pi ppi(2) = 0.d0 ppi(3) = 0.d0 call adjst(ppi) p(1,i) = ppi(1) p(2,i) = ppi(2) p(3,i) = ppi(3) end do ! loop on i write(6,21) j,n,np,psum,p(1,1),p(2,1),p(3,1) ! ! transform call fft2ne(p,lgnp,1.d0) !w write(6,21) j,n,np,psum,p(1,1),p(2,1),p(3,1) ! ! raise transform to power m ! pwr: do i = 1,np rr = p(1,i)*p(1,i) + p(2,i)*p(2,i) if( rr .le. 0.d0 ) cycle pwr ! get log and power up rr = log(rr) * fln / 2.d0 k = rr / al2 rr = rr - k*al2 ! integer part in k, multiply exponent and add k ppi(3) = p(3,i) * fln + k ! fix two floating parts rr = exp(rr) ! get angle t = atan2( p(2,i), p(1,i) ) t = t * fln ppi(1) = rr * cos(t) ppi(2) = rr * sin(t) call adjst( ppi ) p(1,i) = ppi(1) p(2,i) = ppi(2) p(3,i) = ppi(3) end do pwr ! loop on i ! ! ! transform back call fft2ne(p,lgnp,-1.d0) ! ! compute cumulative and print psum = 0. do i = 1,np ! compute probability rr = p(1,i) * exp( p(3,i)*al2 ) psum = psum + rr pj(i,j) = rr !w write(6,21) j,n,i,psum,p(1,i),p(2,i),p(3,i) end do ! loop on i write(6,21) j,n,np,psum,p(1,1),p(2,1),p(3,1) ! done with this j end do ! loop on j ! ! ! write out all three -- but only a fraction ! to file tivmm1.dat for plotting do i = 1,np !w write(6,22) i,(pj(i,j),j=1,3) if( mod(i,16) .eq. 0 ) write(8,22) i,(pj(i,j),j=1,3) end do ! loop on i stop end program tivmm1 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 ! *** end of filename tivmm1.f95 *********************