anchor.diff:  
       real(real64), intent(in) :: alpha
   
  +    ! p(alpha | n) propto p(alpha) * p(n | alpha)
       fwa = fpa(alpha) * exp(log_fpxa(alpha) - log_fpxa_amax) / fpa_amax
   
     end function fwa
  
positive.code: 
      hdir = 0.0_real64
      do i = 1,n_multi 
         a(i) = - multi(i) * (multi_z(i) + alpha) * digamma(multi_z(i) + one + alpha) 
      end do
      hdir = sum(a)
      hdir = hdir / (n_data + alpha * alphabet_size)
      hdir = hdir + digamma(n_data + alpha * alphabet_size + 1.0_real64)
  
    end function hdir
  
  
negative.code: contains
  
    elemental real(real64) function fpa(a) 
      ! prop. to p(alpha) - the prior for alpha in NSB estimator
      use constants
      use gamma_funcs, only: trigamma
      use dirichlet_mod, only: alphabet_size
  