anchor.diff:  
     elemental real(real64) function log_weight(alpha)
       ! un-normalized weight for alpha in the integrals; prop. to p(alpha|x)
  -    use dirichlet_mod, only: log_pna, alpha_prior
  -
  -    real(real64), intent(in) :: alpha
  -
  -    log_weight = log(alpha_prior(alpha)) + log_pna(alpha)
  +    use dirichlet_mod, only: log_pna_u, alpha_prior
  +
  +    real(real64), intent(in) :: alpha
  +
  +    log_weight = log(alpha_prior(alpha)) + log_pna_u(alpha)
   
     end function log_weight
   
  
positive.code:     use gamma_funcs, only: digamma, trigamma, quadgamma
      use dirichlet_mod, only: alphabet_size, n_data, multi,&
           multi_z
      use dirichlet_mod, only: log_pna, alpha_prior
  
      real(real64), intent(in) :: alpha
      real(real64), intent(out) :: logw, dlogw
  
negative.code: 
      ! log weight
      prior = alpha_prior(alpha)
      logw = log(prior) + log_pna_u(alpha)
  
      ! log weight derivative
      dprior = alphabet_size**2 * quadgamma(alphabet_size * alpha + one) - &
  