anchor.diff:      integer(int32) :: i_
       real(real64)   :: a(n_multi)
   
  -    log_fpxa = log_gamma(n_data + one) + log_gamma(alpha * alphabet_size) & 
  +    log_pna = log_gamma(n_data + one) + log_gamma(alpha * alphabet_size) & 
            - alphabet_size * log_gamma(alpha) - log_gamma(n_data + alpha * alphabet_size)
       do i_ = 1, n_multi 
          a(i_) = multi(i_) * (log_gamma(multi_z(i_) + alpha) - log_gamma(multi_z(i_) + one))
       end do
  -    log_fpxa = log_fpxa + sum(a)
  -    
  -  end function log_fpxa
  -
  -  real(real64) function hdir(alpha)
  +    log_pna = log_pna + sum(a)
  +    
  +  end function log_pna
  +
  +  real(real64) function h_bayes(alpha)
       ! posterior average of the entropy given the data and alpha
       ! computed from histogram multiplicities
       use gamma_funcs, only: digamma
  
positive.code:   call initialize_dirichlet(counts, nc)
    call compute_multiplicities(counts)
  
    estimate = hdir(alpha)
  
    call dirichlet_finalize()
  
  
negative.code:     end do
      xs = exp(xs)
      
      fxs = log_weight(xs)
      amax = xs(maxloc(fxs,1,fxs < largest))
      lw_max = log_weight(amax)
      
    end subroutine compute_integration_range
  
    real(real64) function m_func(x)
      use dirichlet_mod, only: h_bayes, integrand
  
      real(real64), intent(in) :: x
      real(real64) :: a, h_c,lw
  
      a = exp(x)
      call integrand(a, h_c, lw)
      m_func = exp(lw - lw_max) * h_c * a
      
    end function m_func
  
    real(real64) function m2_func(x)
      use dirichlet_mod, only: h_bayes, integrand
  
      real(real64), intent(in) :: x
      real(real64) :: a, h_c,lw
  
      a = exp(x)
      call integrand(a, h_c, lw)
      m2_func = exp(lw - lw_max) * h_c**2 * a
      
    end function m2_func
  
    real(real64) function nrm_func(x)
  
      real(real64), intent(in) :: x
      real(real64) :: a
  
      a = exp(x)
      nrm_func = exp(log_weight(a) - lw_max)  * a
      
    end function nrm_func
  
  