anchor.diff:      real(real64), intent(in) :: alpha
       integer(int32) :: i_
   
  -    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)
  -
  -  end function h_bayes
  +    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
  
positive.code:        lw_max = log_weight(amax)
         integrand = exp(log_weight(alpha) - lw_max)  * alpha / amax
      else
         hb = h_bayes(alpha)
  
         integrand = hb**order
  
  
negative.code: 
      err = 0
      if (ascale < 1.e-20) then
         estimate = h_dir(amax)
         err_estimate = 0.0
      else
         call quad(nrm_func,log_alpha1,log_alpha2, nrm, ierr)
  