anchor.diff:      real(real64), intent(in) :: alpha
       integer(int32) :: i_
   
  -    h_dir = - sum(multi * (multi_z + alpha) * digamma(multi_z + alpha + one))
  +    h_dir = - sum(hzs * (hns + alpha) * digamma(hns + alpha + one))
       h_dir = h_dir / (n_data + alpha * alphabet_size)
       h_dir = h_dir + digamma(n_data + alpha * alphabet_size + one)
   
  
positive.code:     ! compute value and derivative of log p(a | x)
      use constants
      use gamma_funcs, only: digamma, trigamma, quadgamma
      use dirichlet_mod, only: alphabet_size, n_data, multi,&
           multi_z
      use dirichlet_mod, only: log_pna_u, alpha_prior
  
      real(real64), intent(in) :: alpha
  
negative.code:     - alphabet_size * digamma(alpha) &
           - alphabet_size * digamma(n_data + alpha * alphabet_size)
  
      wsum = sum(hzs * (digamma(hns + alpha)))
  
      dlpna = dlpna + wsum
  
  