anchor.diff:      real(float64) :: c, nu, ni, xi, jsum
   
       nu = n_data + alpha * alphabet_size
  -    phi = digamma(hn + alpha + one) - &
  +    phi = digamma(nk + alpha + one) - &
            digamma(nu + two)
       c = trigamma(nu + two)
   
       h_var = 0.0_float64
  -    do i_ = 1, size(hz)
  -       ni = hn(i_) + alpha
  +    do i_ = 1, size(zk)
  +       ni = nk(i_) + alpha
          xi = phi(i_)
  -       jsum = sum(hz * ni * (hn + alpha) * &
  +       jsum = sum(zk * ni * (nk + alpha) * &
               (xi * phi - c))
  -       h_var = h_var + hz(i_) * jsum
  -       h_var = h_var - hz(i_) * ni**2 * (xi**2 - c)
  +       h_var = h_var + zk(i_) * jsum
  +       h_var = h_var - zk(i_) * ni**2 * (xi**2 - c)
          xi = xi + 1 / (ni + one)
  -       h_var = h_var + hz(i_) * (ni + one) * ni * &
  +       h_var = h_var + zk(i_) * (ni + one) * ni * &
               (xi**2 + trigamma(ni + two) - c)
       end do
   
  
positive.code:   implicit none
  
    integer, intent(in)  :: n
    real(float64), intent(in)    :: hn1(n)
    real(float64), intent(in)    :: hz1(n)
    real(float64), intent(in)    :: nc
    real(float64),   intent(out) :: estimate
    real(float64),   intent(out) :: err_estimate
    integer :: err
  
    call initialize_from_multiplicities(hn1, hz1, nc)
  
    call compute_integration_range()
  
  
negative.code: 
  end module constants
  
  
  module counter
    use constants
    implicit none
    real(float64), allocatable :: nk(:)
    real(float64), allocatable :: zk(:)
    real(float64) :: n_data
    real(float64) :: n_bins
    real(float64) :: k1
  contains
  
    subroutine counts_reset()
      implicit none
      n_data = 0.0_float64
      n_bins = 0.0_float64
      k1 = 0.0_float64
      if (allocated(nk)) then
         deallocate(nk)
      end if
      if (allocated(zk)) then
         deallocate(zk)
      end if
    end subroutine counts_reset
  
    subroutine counts_fit(ar)
      implicit none
      integer, intent(in) :: ar(:)
      integer :: i, j, u
      integer :: x, err
      integer :: xmax
      integer, allocatable :: wrk(:)
  
      call counts_reset()
  
      xmax = maxval(ar)
      allocate(wrk(0:xmax), stat=err)
      wrk = 0
      u = 0
      do i = 1,ubound(ar, 1)
         x = ar(i)
         if (wrk(x) == 0) then
            u = u + 1
         end if
         wrk(x) = wrk(x) + 1
      end do
  
      allocate(nk(u), zk(u), stat=err)
  
      u = 0
      do i = 0, xmax
         x = wrk(i)
         if (x > 0) then
            u = u + 1
            nk(u) = i
            zk(u) = x
         end if
      end do
  
      n_data = sum(zk * nk)
      n_bins = sum(zk)
      k1 = sum(zk, nk>0)
  
      deallocate(wrk)
    end subroutine counts_fit
  
    subroutine add_empty_bins(alphabet_size)
      ! add empty bins to nk, zk
      implicit none
      real(float64), intent(in) :: alphabet_size
      integer :: k, err
      real(float64) :: unobserved
      integer, allocatable :: wrk(:)
  
      if (allocated(nk)) then
         if (alphabet_size > n_bins) then
            unobserved = alphabet_size - n_bins*1.0_float64
            if (minval(nk) < 1.e-3_float64) then
               k = minloc(nk, 1)
               zk(k) = zk(k) + unobserved
            else
               k = size(nk)
               allocate(wrk(k), stat=err)
  
               wrk(:k) = nk
               deallocate(nk)
               allocate(nk(k+1), stat=err)
               nk(:k) = wrk
               nk(k+1) = 0
  
               wrk(:k) = zk
               deallocate(zk)
               allocate(zk(k+1), stat=err)
               zk(:k) = wrk
               zk(k+1) = unobserved
  
               deallocate(wrk)
            end if
         end if
      end if
  
    end subroutine add_empty_bins
  end module counter
  
  
  module dirichlet_mod
    use constants
    use counter
    implicit none
  
    real(float64)              :: alphabet_size
    real(float64), allocatable :: phi(:)  ! wrk array for var
  
  contains
  
    subroutine initialize_from_counts(counts, nc)
      ! set n_multi, nk, multi
      integer, intent(in) :: counts(:)
      real(float64), intent(in) :: nc
      integer :: err
  
      alphabet_size = nc
  
      call counts_reset()
  
      call counts_fit(counts)
  
      call add_empty_bins(nc)
  
      allocate(phi(size(nk)), stat=err)
  
    end subroutine initialize_from_counts
  
    subroutine initialize_from_multiplicities(nk1, zk1, nc)
      ! set nk, zk, n_data, alphabet_size
      real(float64), intent(in) :: nk1(:)
      real(float64), intent(in) :: zk1(:)
      real(float64), intent(in) :: nc
      integer :: idx=-1, err, j, nm
  
      nm = size(nk1)
      ! check if zeros are included in mults arrays
      if (any(int(nk1) == 0)) then  ! we work with floats
         allocate(nk(nm), stat=err)
         allocate(zk(nm), stat=err)
         allocate(phi(nm), stat=err)
         nk = nk1
         zk = zk1
         idx = -1
         do j = 1, nm
            if (int(nk1(j)) == 0) then
               idx = j
               exit
            end if
         end do
         zk(idx) = zk1(idx) + nc - sum(zk1)
      else
         allocate(nk(nm + 1), stat=err)
         allocate(zk(nm + 1), stat=err)
         allocate(phi(nm + 1), stat=err)
         nk(1) = 0._float64
         zk(1) = nc - sum(zk1)
         nk(2:) = nk1
         zk(2:) = zk1
      end if
  
      alphabet_size = nc
      n_data = sum(zk * nk)
  
    end subroutine initialize_from_multiplicities
  
    subroutine finalize()
  
      if (allocated(zk)) then
         deallocate(zk)
      end if
  
      if (allocated(nk)) then
         deallocate(nk)
      end if
  
      if (allocated(phi)) then
  