anchor.diff:      log_pna_u = log_gamma(alpha * alphabet_size) &
            - alphabet_size * log_gamma(alpha) &
            - log_gamma(n_data + alpha * alphabet_size) &
  -         + sum(hzs * (log_gamma(hns + alpha)))
  +         + sum(hz * (log_gamma(hn + alpha)))
   
     end function log_pna_u
   
  
positive.code:     real(real64) :: c, nu, ni, xi, jsum
  
      nu = n_data + alpha * alphabet_size
      phi = digamma(hns + alpha + one) - &
           digamma(nu + two)
      c = trigamma(nu + two)
  
      h_var = 0.0
      do i_ = 0, size(hzs)-1
         ni = hns(i_) + alpha
         xi = phi(i_)
         jsum = sum(hzs * ni * (hns + alpha) * &
              (xi * phi - c))
         h_var = h_var + hzs(i_) * jsum
         h_var = h_var - hzs(i_) * ni**2 * (xi**2 - c)
         xi = xi + 1 / (ni + one)
         h_var = h_var + hzs(i_) * (ni + one) * ni * &
              (xi**2 + trigamma(ni + two) - c)
      end do
  
  
negative.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, hz,&
           hn
      use dirichlet_mod, only: log_pna_u, alpha_prior
  
      real(real64), intent(in) :: alpha
  