anchor.diff:      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 = - n_empty_bins * alpha * digamma(alpha + one)
  +    h_bayes = 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)
   
  
positive.code:          - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size)
  
      wsum = sum(multi * (log_gamma(multi_z + alpha) - log_gamma(multi_z + one)))
  
      log_pna = log_pna + wsum
      
  
negative.code:   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(:)
  
  