anchor.diff:    call initialize_dirichlet(counts, nc)
     call compute_multiplicities(counts)
   
  -  estimate = h_bayes(alpha)
  +  estimate = h_dir(alpha)
   
     call dirichlet_finalize()
   
  
positive.code: 
      err = 0
      if (ascale < 1.e-20) then
         estimate = h_bayes(amax)
         err_estimate = 0.0
      else
         call quad(nrm_func,log_alpha1,log_alpha2, nrm, ierr)
  
negative.code:     real(real64), intent(in) :: alpha
      integer(int32) :: i_
  
      h_dir = - sum(multi * (multi_z + alpha) * digamma(multi_z + alpha + one))
      h_dir = h_dir / (n_data + alpha * alphabet_size)
      h_dir = h_dir + digamma(n_data + alpha * alphabet_size + one)
  
    end function h_dir
  
    elemental real(real64) function integrand(alpha, amax, order)
      ! posterior average of the entropy given the data and alpha
  