anchor.diff:  
     end function h_dir
   
  -
  -  elemental real(real64) function integrand(alpha, amax, order)
  +  real(real64) function h_var(alpha)
  +    ! posterior average of the entropy given data and a specific alpha value
  +    ! computed from histogram multiplicities
  +    use gamma_funcs, only: digamma, trigamma
  +    use constants
  +
  +    real(real64), intent(in) :: alpha
  +    integer(int32) :: i_
  +    real(real64) :: c, nu, ni, xi, jsum
  +
  +    nu = n_data + alpha * alphabet_size
  +    phi = digamma(multi_z + alpha + one) - &
  +         digamma(nu + two)
  +    c = trigamma(nu + two)
  +
  +    h_var = 0.0
  +    do i_ = 0, size(multi)-1
  +       ni = multi_z(i_) + alpha
  +       xi = phi(i_)
  +       jsum = sum(multi * ni * (multi_z + alpha) * &
  +            (xi * phi - c))
  +       h_var = h_var + multi(i_) * jsum
  +       h_var = h_var - multi(i_) * ni**2 * (xi**2 - c)
  +       xi = xi + 1 / (ni + one)
  +       h_var = h_var + multi(i_) * (ni + one) * ni * &
  +            (xi**2 + trigamma(ni + two) - c)
  +    end do
  +
  +    h_var = h_var / (nu * (nu + one))
  +
  +  end function h_var
  +
  +
  +  real(real64) function integrand(alpha, amax, order)
       ! posterior average of the entropy given the data and alpha
       ! computed from histogram multiplicities
       use gamma_funcs, only: digamma
  
positive.code:     ! set n_multi, multi_z, multi
      integer(int32), intent(in) :: counts(:)
      real(real64), intent(in) :: nc
  
      alphabet_size = nc
      n_data = sum(counts)
  
negative.code:   end subroutine weight_std
  
    subroutine hnsb(estimate,err_estimate, err)
      use dirichlet_mod, only: h_dir, h_var
      real(real64), intent(out) :: estimate,err_estimate
      integer(int32), intent(out) :: err
      real(real64)              :: rslt,nrm
  