anchor.diff:  contains
   
     elemental real(real64) function fpa(a) 
  -    ! p(alpha) prior for alpha in NSB estimator
  +    ! prop. to p(alpha) - the prior for alpha in NSB estimator
       use constants
       use gamma_funcs, only: trigamma
       use dirichlet_mod, only: alphabet_size
  
positive.code:     ! recompute fxs
      fxs = exp(lfwa(xs))
  
      ! recompute integration range
      a1 = alpha2
      a2 = alpha1
      do i = 1,nx
         x = log(xs(i))
         f = fxs(i)
         if (x > alpha1 .and. x < a1 .and. f > small_number) a1 = log(xs(i-1))
         if (x > a2 .and. x < alpha2 .and. f > small_number) a2 = log(xs(i+1))
      end do
      
      alpha1 = a1
      alpha2 = a2 
  
      ! recompute the maximum
      dx = (alpha2 - alpha1) / (nx * 1.0_real64)
  
negative.code:     alpha1 = -20.0_real64
      alpha2 = 10.0_real64
  
      ! initialize amax
      amax = 1.0_real64
      log_fpxa_amax = log_fpxa(amax) ! log p(n|a_max)
      fpa_amax = fpa(amax) ! propto p(a)
  
      ! set intervals equally spaced on log scale
      dx = (alpha2 - alpha1) / real(nx -1.0_real64)
      do i = 1,nx
         xs(i) = alpha1 + (i - 1) * dx
      end do
      xs = exp(xs)
      
      fxs = lfwa(xs)
      ! find amax such that the weight in the integral - fwa(alpha) - is maximal
      i = maxloc(fxs,1,fxs < largest)
      amax = xs(i)
      log_fpxa_amax = log_fpxa(amax)
  