anchor.diff:    integer(int32)              :: alphabet_size
     integer(int32)              :: n_data
     integer(int32)              :: n_multi
  +  integer(int32)              :: n_empty_bins
     integer(int32), allocatable :: multi_z(:)
     integer(int32), allocatable :: multi(:)
   
  
positive.code:     real(real64), intent(in) :: alpha
      integer(int32) :: i_
  
      h_bayes = 0.0_real64
      h_bayes = - sum(multi * (multi_z + alpha) * digamma(multi_z + alpha + one))
      h_bayes = h_bayes / (n_data + alpha * alphabet_size)
      h_bayes = h_bayes + digamma(n_data + alpha * alphabet_size + one)
  
  
negative.code:          - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size)
  
      asum = n_empty_bins * alpha * digamma(alpha + one)
      asum = asum + sum(multi * (multi_z + alpha) &
           * digamma(multi_z + alpha + one))
  
      bsum = n_empty_bins * (log_gamma(alpha) - log_gamma(one))
      bsum = bsum + sum(multi * (log_gamma(multi_z + alpha) &
           - log_gamma(multi_z + one)))
  
      lpna = lpna + bsum
  