anchor.diff:  
       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)
  +       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 + 1.0_real64)
  +    hdir = hdir + digamma(n_data + alpha * alphabet_size + one)
   
     end function hdir
   
  
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:   end function log_fpxa
  
    real(real64) function hdir(alpha)
      ! posterior average of the entropy given the data and alpha
      ! computed from histogram multiplicities
      use gamma_funcs, only: digamma
      use constants
      
  