anchor.diff:      
       real(real64), intent(in) :: alpha
       integer(int32) :: i
  -    real(real64)   :: a(ngtz)
  +    real(real64)   :: a(nmulti)
   
       hdir = 0.0_real64
  -    do i = 1,ngtz 
  -       a(i) = - multi_gtz(i) * (fgtz(i) - one + alpha) * digamma(fgtz(i) + alpha) 
  +    do i = 1,nmulti 
  +       a(i) = - multi(i) * (multi_x(i) - one + alpha) * digamma(multi_x(i) + alpha) 
       end do
       hdir = sum(a)
       hdir = hdir / (ndata + alpha * alphabet_size)
  
positive.code:     
      real(real64), intent(in) :: alpha
      integer(int32) :: i
      real(real64)   :: a(ngtz)
  
      log_fpxa = log_gamma(ndata + one) + log_gamma(alpha * alphabet_size) & 
           - alphabet_size * log_gamma(alpha) - log_gamma(ndata + alpha * alphabet_size)
      do i = 1,ngtz 
         a(i) = multi_gtz(i) * (log_gamma(fgtz(i) - one + alpha) - log_gamma(fgtz(i) * one))
      end do
      log_fpxa = log_fpxa + sum(a)
      
  
negative.code: 
    integer(int32)              :: alphabet_size
    integer(int32)              :: ndata
    integer(int32)              :: nmulti
    integer(int32), allocatable :: multi_x(:)
    integer(int32), allocatable :: multi(:)
  
  contains
  
  