anchor.diff:      nbins = size(counts)
       nmax = maxval(counts)
       allocate(multi0(nmax+1),stat=err)
  +    ! multi0(n+1) is the number of states with frequency n
       multi0 = 0
  +    ! take into account the alphabet_size - nbins states with zero frequency
       multi0(1) = alphabet_size - nbins
       do i_ = 1,nbins
          ni_ = counts(i_)
  
positive.code:   end subroutine compute_multiplicities
  
    subroutine dirichlet_finalize()
  
      deallocate(multi_z,multi)
  
    end subroutine dirichlet_finalize
  
    pure real(real64) function log_fpxa(alpha) 
  
negative.code:     use constants
      
      real(real64), intent(in) :: alpha
      integer(int32) :: i_
      real(real64)   :: a(n_multi)
  
      log_fpxa = log_gamma(n_data + one) + log_gamma(alpha * alphabet_size) & 
           - alphabet_size * log_gamma(alpha) - log_gamma(n_data + alpha * alphabet_size)
      do i_ = 1, n_multi 
         a(i_) = multi(i_) * (log_gamma(multi_z(i_) - one + alpha) - log_gamma(multi_z(i_) * one))
      end do
      log_fpxa = log_fpxa + sum(a)
      
  