!*****************
! File Description:  
!   
! Creator: Hoang-Trong Minh Tuan 
! 
! Rev: 0.1 
! Update: 
!    0.2: early 2010 - Tuan
!         rename combinator_me to count_ball2bin
!         add count_ball2bin_T (transpose)
!    0.3: (plan) add date-time function to determine expiration use of the program
!*******************/
MODULE math_utility
  USE matrix_utility
  USE PRECISIONINFO
  IMPLICIT NONE

  SAVE

CONTAINS
  !! =================
  !! factorial of N = N!
  !! =================
  FUNCTION fact(N)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: N
    REAL(KIND=REAL8) :: fact
    INTEGER :: I

    IF (N == 0) THEN
       fact = 1
    ELSE
       fact = 1.
       DO I = 1, N, 1
          fact = fact * I
       END DO
    END IF
  END FUNCTION fact

  !! ===========================
  !! return (N+M-1)!/[(N-1)!*M!]
  !! ===========================
  FUNCTION combination(M,N)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: M,N
    INTEGER :: combination
    REAL(KIND=REAL8) :: tmp

    tmp = fact(N+M-1)/(fact(N-1)*fact(M))
    combination = FLOOR(tmp)
  END FUNCTION combination


  !! =======================================
  !! arrange M balls into N bins
  !! Output: matrix of size C(N+M-1,N) x N  
  !! =======================================
!  SUBROUTINE combinator_mem(M, N, matrix, rows, cols)
  SUBROUTINE count_ball2bin(M, N, matrix, rows, cols)
    USE matrix_utility
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: M, N
    INTEGER, ALLOCATABLE ,DIMENSION(:,:) :: matrix
    INTEGER, INTENT(OUT) :: rows, cols
    !local
    INTEGER :: count, test
    INTEGER :: balls, bins

    balls = M
    bins = N

    IF (ALLOCATED(matrix)) DEALLOCATE(matrix)

    rows = combination(M,N)
    cols = N
    !     print *, "Number of combinations: ", rows
    ALLOCATE(matrix(rows,N), STAT = test)
    IF (test /= 0) THEN
       PRINT *, "cannot allocate memory"
       STOP
    END IF
    CALL mat_reset(matrix, rows, cols)
    
    count = 1
    CALL combinator_mem_recursive(M,N,matrix)

  CONTAINS
    !! =====================
    !! part of count_ball2bin
    !! =====================
    RECURSIVE SUBROUTINE combinator_mem_recursive(M,N,matrix)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: M, N
      INTEGER, DIMENSION(:,:), INTENT(INOUT) :: matrix
      !local
      INTEGER :: I, K

      DO I = M, 0, -1
         IF (N == 1) THEN
            matrix(count, bins-N+1) = I
            count = count + 1
            RETURN
         ELSE
            IF (M == 0 .AND. N > 1) THEN
               DO K = bins-N+1, bins, 1
                  matrix(count,K) = 0
               END DO
               count = count + 1
            ELSE
               matrix(count, bins-N+1) = I

               CALL combinator_mem_recursive(M-I, N-1, matrix)             
               IF (I > 0) THEN
                  matrix(count,1:bins-N) = matrix(count-1,1:bins-N)
               END IF
            END IF
         END IF
      END DO
    END SUBROUTINE combinator_mem_recursive
  END SUBROUTINE count_ball2bin


  !! =======================================
  !! Arrange M balls into N bins.
  !! Now RETURN a transpose matrix of count_ball2bin
  !!
  !! Output: matrix of size N x C(N+M-1,N)
  !! =======================================
  SUBROUTINE count_ball2bin_T(M, N, matrix, rows, cols)
    USE matrix_utility
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: M, N
    INTEGER, ALLOCATABLE ,DIMENSION(:,:) :: matrix
    INTEGER, INTENT(OUT) :: rows, cols
    !local
    INTEGER :: count, test
    INTEGER :: balls, bins
    balls = M
    bins = N

    IF (.NOT. ALLOCATED(matrix)) THEN
       cols = combination(M,N)
       rows = N
       !     print *, "Number of combinations: ", rows
       ALLOCATE(matrix(N, cols), STAT = test)
       IF (test /= 0) THEN
          PRINT *, "cannot allocate memory"
          STOP
       END IF
       CALL mat_reset(matrix, rows, cols)
    END IF

    count = 1
    CALL new_combinator_mem_recursive(M,N,matrix)

  CONTAINS
    ! =====================
    ! part of count_ball2bin_T
    ! =====================
    RECURSIVE SUBROUTINE new_combinator_mem_recursive(M,N,matrix)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: M, N
      INTEGER, DIMENSION(:,:), INTENT(INOUT) :: matrix
      !local
      INTEGER :: I, K

      DO I = M, 0, -1
         IF (N == 1) THEN
            matrix(bins-N+1, count) = I
            count = count + 1
            RETURN
         ELSE
            IF (M == 0 .AND. N > 1) THEN
               DO K = bins-N+1, bins, 1
                  matrix(K, count) = 0
               END DO
               count = count + 1
            ELSE
               matrix(bins-N+1, count) = I
               CALL new_combinator_mem_recursive(M-I, N-1, matrix)
               IF (I > 0) THEN
                  matrix(1:bins-N,count) = matrix(1:bins-N,count-1)
               END IF
            END IF
         END IF
      END DO
    END SUBROUTINE new_combinator_mem_recursive
  END SUBROUTINE count_ball2bin_T

  !! =================================
  !! Kronecker product
  !! =================================
  SUBROUTINE kron(a,b,z,m1,m2,n1,n2)
    IMPLICIT NONE
    REAL(KIND=dp):: a(m1,m2),b(n1,n2),z(n1*m1,n2*m2)
    INTEGER :: m1, m2, n1, n2
    INTEGER :: i, j
    !$acc region
    DO j =1,m2*n2,n2
       DO i=1,m1*n1, n1
          z(i:i+n1-1,j:j+n2-1) = a((i-1)/n1+1,(j-1)/n2+1) * b(:,:)
       ENDDO
    ENDDO
    !$acc end region
  END SUBROUTINE kron

  LOGICAL FUNCTION isnan(a)
    REAL(KIND=dp) ::  a
    if (a.ne.a) then
       isnan = .true.
    else
       isnan = .false.
    end if
    return
  END FUNCTION isnan

  SUBROUTINE check_expiration()
    INTEGER*4, DIMENSION(3) :: expired_date
    INTEGER*4, DIMENSION(3) :: today, now
    LOGICAL :: is_expired 
    
    CALL idate(today)   ! today(1)=day, (2)=month, (3)=year
    CALL itime(now)     ! now(1)=hour, (2)=minute, (3)=second
    !! 04/01/1995; time 14:07:40
    expired_date = (/ 10, 6,  2012 /)  !month,day,year
    
    is_expired = .TRUE.

    !IF (today(1) .LE. expired_date(1) .AND. today(2) .LE. expired_date(2) .AND. today(3) .LE. expired_date(3)) THEN
    IF ((today(3) .LT. expired_date(3)) .OR. &
         (today(3) .EQ. expired_date(3) .AND. today(2) .LT. expired_date(2)) .OR. &
         (today(3) .EQ. expired_date(3) .AND. today(2) .EQ. expired_date(2) .AND. today(1) .LE. expired_date(1)) &
         ) THEN 
       is_expired = .FALSE.
    ENDIF

    IF (is_expired) THEN
       PRINT *, "The program is expired. Please contact the author"
       STOP
    END IF
  END SUBROUTINE  check_expiration
END MODULE math_utility


