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