anchor.diff:      integer(int32)              :: err
       integer(int32)              :: nmax
       integer(int32), allocatable :: multi0(:)
  +    integer(int32)              :: n_multi
   
       ! compute multiplicities
       ! nmax is the largest number of samples in a bin
  
positive.code:     real(real64), intent(in) :: alpha
      integer(int32) :: i_
  
      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)
  
  
negative.code:     - alphabet_size * digamma(alpha) &
           - alphabet_size * digamma(n_data + alpha * alphabet_size)
  
      wsum = sum(multi * (digamma(multi_z + alpha)))
  
      dlpna = dlpna + wsum
  
  