!     Last change:  DOS  24 Jul 2000    8:29 pm
!                                       *** copyright 2000 ***
! *** filename sum24.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM SUM24
!  TEST OF ADDITION ALGORITHM THAT ADDS FROM 1 TO N BY PAIRING          
!                                                                       
      IMPLICIT NONE
      INTEGER, PARAMETER :: N = 32768   ! HOW MANY TO ADD
      INTEGER, PARAMETER :: BMAX = 30   ! N < 2**BMAX (DIMENSION OF B)
      REAL XXA,YYA,ZZA,XX4,YY4,ZZ4,XI,YI,ZI
      REAL, PARAMETER :: G = .5772157   ! EULER'S CONSTANT
      REAL, DIMENSION(N) :: X,Y,Z
      REAL, DIMENSION(BMAX) :: BX,BY,BZ
      INTEGER i
      INTEGER IXXA,IX1,IXX4,IYYA,IY1,IYY4,IZZA,IZ1,IZZ4
      EQUIVALENCE (XXA,IXXA),(X(1),IX1),(XX4,IXX4)
      EQUIVALENCE (YYA,IYYA),(Y(1),IY1),(YY4,IYY4)
      EQUIVALENCE (ZZA,IZZA),(Z(1),IZ1),(ZZ4,IZZ4)
!
      INTERFACE
      SUBROUTINE SUM4(X,B)
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, DIMENSION(:), INTENT(IN OUT) :: B
      END SUBROUTINE SUM4
      END INTERFACE
!
  21  FORMAT(' N =',I12,4X,'ANALYTIC SUM',E16.9,Z12/        &
     &          'SUM2',E16.9,Z12,4X,'SUM4',E16.9,Z12)
  25  FORMAT(/' FIRST EXAMPLE - ADD NUMBERS FROM 1 TO N')
  26  FORMAT(/' SECOND EXAMPLE - ADD SQUARES FROM 1 TO N')
  27  FORMAT(/' THIRD EXAMPLE - HARMONIC SUM')
!                   WRITE RESULTS HERE
      OPEN( UNIT=6, file='sum24.out' )
!                   DON'T FORGET TO INITIALIZE SUMMING ARRAYS
!                    BX, BY, BZ TO ZERO FOR SUM4
      BX = 0.
      BY = 0.
      BZ = 0.
!                                                                       
!  BOTH SUM2 AND SUM4 ADD THE SAME WAY, BUT DIFFERENT STORAGE
!   SUM2 NEEDS ALL NUMBERS IN ONE LONG LIST
!   SUM4 CAN WORK IN A LOOP AS NUMBERS ARE CREATED
!         BUT HAS A LIMIT OF 2**M NUMBERS (LIKE M=30) 
! 
!            ADD USING SUM4 ON THE FLY AND FILL UP X, Y AND Z
!             FOR SUM2 TO USE
      DO I = 1,N                                                      
!            FIRST ADD NUMBERS FROM 1 TO N
      X(I) = I                                                          
      XI = X(I)                                                         
      CALL SUM4(XI,BX)                                               
!            SECOND EXAMPLE - ADD SQUARES FROM 1 TO N
      Y(I) = X(I)*X(I)                                                  
      YI = Y(I)                                                         
      CALL SUM4(YI,BY)                                               
!            THIRD EXAMPLE - HARMONIC SERIES
      Z(I) = 1.0/REAL(N+1-I)                                           
      ZI = Z(I)                                                         
      CALL SUM4(ZI,BZ)                                               
      ENDDO                   ! END LOOP ON I                                                          
! 
!  GET RESULTS FROM SUM4  
!
      XX4 = 0.                                                           
      YY4 = 0.                                                           
      ZZ4 = 0.                                                           
      DO I = 1,BMAX
      XX4 = XX4 + BX(I)
      YY4 = YY4 + BY(I)
      ZZ4 = ZZ4 + BZ(I)
      ENDDO                   ! END LOOP ON I                                        
!
!  COMPARE RESULTS WITH ANALYTIC SUMMATION
!
      WRITE(6,25)
      CALL SUM2(X,N)                                                    
      XXA = REAL(N)*REAL(N+1)/2.
      WRITE(6,21) N,XXA,IXXA,X(1),IX1,XX4,IXX4                            
!
      WRITE(6,26) 
      CALL SUM2(Y,N)                                                    
      YYA = REAL(N)*REAL(N+1)*REAL(2*N+1)/6.
      WRITE(6,21) N,YYA,IYYA,Y(1),IY1,YY4,IYY4                   
!
      WRITE(6,27) 
      CALL SUM2(Z,N)                                                    
      ZZA = G + LOG(REAL(N))
      WRITE(6,21) N,ZZA,IZZA,Z(1),IZ1,ZZ4,IZZ4                            
      STOP
      END PROGRAM SUM24
      SUBROUTINE SUM2(X,N)                                              
!  RETURNS SUM OF X(I) FROM 1 TO N IN X(1)                              
!  ADDS PAIRS TOGETHER AT EACH STEP
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION( N ), INTENT(IN OUT) :: X
      INTEGER INCR,SKIP,J,NN
!                  BRANCH FOR SPECIAL CASE                              
      IF( N .EQ. 1) RETURN                                              
!                  INITIALIZATION                                       
      INCR = 1                                                          
      SKIP = 2*INCR                                                     
      NN = N-1                                                          
!                  RETURN HERE AFTER PASS THROUGH DATA                  
      DO                                                          
        DO J = 1,NN,SKIP                                                
        X(J) = X(J) + X(J+INCR)                                           
        ENDDO                   ! END LOOP ON J                                                          
      INCR = 2*INCR                                                     
      SKIP = 2*INCR                                                     
      IF( INCR .GT. NN ) EXIT
      ENDDO
!                  PASS THROUGH IF AFTER LAST ADDITIONS                 
      RETURN                                                            
      END SUBROUTINE SUM2
      SUBROUTINE SUM4(X,B)                                            
!  THIS ROUTINE ADD NUMBERS IN PAIRS AND KEEPS TRACK IN ARRAY B         
!  GREAT FOR PUTTING IN A LOOP                                          
!  WHEN DONE, JUST ADD UP THE VECTOR B IN INCREASING INDEX ORDER        
!  *** WARNINGS ***
!      1) B MUST BE INITIALIZED TO ZERO
!      2) NUMBER OF ITEMS ADDED CAN'T EXCEED 2**(LENGTH OF B) - 1
!
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, DIMENSION(:), INTENT(IN OUT) :: B
      INTEGER I
      REAL S
!                   INITIALIZE INDEX AND SUM                                                                      
      I = 1
      S = X
      DO                                                          
      IF( B(I) .EQ. 0. ) EXIT
!             IF NOT EMPTY, ADD, ZERO OUT AND MOVE UP IN B              
      S = S + B(I)
      B(I) = 0.
      I = I + 1
      ENDDO                                            
!                   ERROR IF INDEX CLIMBS HIGHER THAN SIZE OF B
!                    *** NOT CHECKED ***
      B(I) = S
      RETURN                                                            
      END SUBROUTINE SUM4
! *** end of filename sum24.f95         *********************