! Last change: DOS 3 Aug 2000 5:53 pm ! *** copyright 2000 *** ! *** filename chfniv.f95 *** John F. Monahan ** ! ********************** program chfniv ! demonstration chfniv -- straight inversion of characteristic fn ! inverse of characteristic function to get probabilities implicit none real, parameter :: omgam = .4227843 ! 1 - Euler's gamma real, parameter :: twopi = 6.2831853 real, parameter :: pio2 = 1.570796 integer, parameter :: lgnp = 12 ! log2(np) integer, parameter :: np = 4096 ! np = number of points complex, dimension(np) :: a,b real t,c,s,deltx,deltt,angle integer j,jj ! 20 format(' Demonstration of inverting characteristic function' & & /30x,'CDF from',5x,'CDF from' & & /10x,'x',8x,'Density',5x,'char fn',7x,'sum') 21 format(2x,4f12.6) ! open( unit=6, file='chfniv.out' ) ! ! delta(x) and delta(t), product is 2pi/np ! range of x is +/- np*deltx/2 = +/- 128 ! range of t is +/- np*deltt/2 = +/- 8*2pi deltx = 1/16. deltt = 6.283185/256. ! ! load vector with characteristic function values ! ! start at 2 do jj = 2,np j = jj - 1 if( jj .gt. np/2 ) j = jj - np - 1 t = deltt*real(j) ! characteristic function, omgam = 1 - Euler's gamma ! angle = t* ( omgam + log(abs(t)) ) c = cos(angle) s = sin(angle) ! characteristic function * deltt / 2*pi a(jj) = cmplx(c,-s) * deltt * exp(-abs(t)*pio2) / twopi ! difference with scaled Cauchy b(jj) = cmplx(1.-c,s) * deltt * exp(-abs(t)*pio2)/(twopi*t) end do ! loop on jj ! treat first one ( t = 0 ) special a(1) = deltt/twopi ! what's really the limit as t goes to zero? b(1) = cmplx( 0., 0.) ! ! compute transforms call fft2n(a,lgnp,1.) call fft2n(b,lgnp,1.) ! write header write(6,20) ! ! write them out nicely s = 0. do jj = 1,np j = np/2 + jj if( jj .gt. np/2 ) j = j - np t = deltx*real(jj-np/2-1) ! ! skip if they're too far away if( t .lt. -8. ) cycle if( t .gt. 32. ) cycle ! ! simple sum to approximate cdf s = s + deltx*real(a(j)) ! ! add in Cauchy df c = .5 + atan(t/pio2)/(twopi/2.) - aimag(b(j)) write(6,21) t,real(a(j)),c,s end do ! loop on jj stop end program chfniv SUBROUTINE FFT2N(A,N,SGN) ! ! 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, APRIL 2000 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 COMPLEX, DIMENSION(1:N), INTENT(IN OUT) :: A REAL, INTENT(IN) :: SGN REAL C,S COMPLEX, DIMENSION(25) :: WK COMPLEX WS,ALEFT,ARIGHT 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. ) WK(2) = CMPLX( 0., -SGN ) ! 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 ) 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., 0.) 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 ALEFT = A( L + J ) ARIGHT = WS*A( L + J + M2 ) A( L + J ) = ALEFT + ARIGHT A( L + J + M2 ) = ALEFT - ARIGHT END DO ! LOOP ON J ! 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 ALEFT = A(LL) ARIGHT = A(LM+1) A(LL) = ARIGHT A(LM+1) = ALEFT 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 END DO ! LOOP ON I RETURN END SUBROUTINE FFT2N ! *** end of filename chfniv.f95 *********************