anchor.diff:      end do
       xs = exp(xs)
       
  -    fxs = log_fwa(xs)
  +    fxs = log_weight(xs)
       amax = xs(maxloc(fxs,1,fxs < largest))
  -    log_fwa_amax = log_fwa(amax)
  +    lw_max = log_weight(amax)
       
     end subroutine compute_integration_range
   
     real(real64) function m_func(x)
  -    ! integrating after change of variable x = log(alpha) in the integral(s)
  -    use dirichlet_mod, only: hdir
  +    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)
  -    m_func = exp(log_fwa(a) - log_fwa_amax) * hdir(a) * a
  -    
  -  end function m_func
  -
  -  real(real64) function m2_func(x)
  -    ! integrating after change of variable x = log(alpha) in the integral(s)
  -    use dirichlet_mod, only: hdir
  -
  -    real(real64), intent(in) :: x
  -    real(real64) :: a
  -
  -    a = exp(x)
  -    m2_func = exp(log_fwa(a) - log_fwa_amax) * hdir(a)**2 * a
  -    
  -  end function m2_func
  -
  -  real(real64) function nrm_func(x)
  -    ! integrating after change of variable x = log(alpha) in the integral(s)
  -
  -    real(real64), intent(in) :: x
  -    real(real64) :: a
  -
  -    a = exp(x)
  -    nrm_func = exp(log_fwa(a) - log_fwa_amax)  * a
  +    nrm_func = exp(log_weight(a) - lw_max)  * a
       
     end function nrm_func
   
  
positive.code:     integer(int32) :: i
      real(real64)   :: a(n_multi)
  
      hdir = 0.0_real64
      do i = 1,n_multi 
         a(i) = multi(i) * (multi_z(i) + alpha) * digamma(multi_z(i) + alpha + one) 
      end do
      hdir = - sum(a)
      hdir = hdir / (n_data + alpha * alphabet_size)
      hdir = hdir + digamma(n_data + alpha * alphabet_size + one)
  
    end function hdir
  
  end module dirichlet_mod
  
  module nsb_mod
    use iso_fortran_env
    implicit none
  
    real(real64) :: log_alpha1
    real(real64) :: log_alpha2
    real(real64) :: amax
    real(real64) :: log_fwa_amax
  
  contains
  
    elemental real(real64) function log_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
      
      real(real64), intent(in) :: a
      
  
negative.code:     end do
      xs = exp(xs)
      
      fxs = log_weight(xs)
      ! find amax such that the weight in the integral - fwa(alpha) - is maximal
      i = maxloc(fxs,1,fxs < largest)
      amax = xs(i)
      lw_max = log_weight(amax)
  
      ! recompute fxs
      fxs = exp(log_weight(xs) - lw_max)
  
      ! re-compute a reasonable integration range
      log_alpha1 = log(minval(xs, fxs > small_number))
  