anchor.diff:  
   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
  -  implicit none
  -
  -  integer              :: n_data
  -  real(float64)                :: alphabet_size
  -  real(float64), allocatable :: hn(:)  ! array of observed frequencies
  -  real(float64), allocatable :: hz(:)  ! multiplicities of frequency z
  +  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, hn, multi
  +    ! set n_multi, nk, multi
       integer, intent(in) :: counts(:)
       real(float64), intent(in) :: nc
  -    integer              :: nbins
  -    integer              :: i_,k_,ni_
  -    integer              :: err
  -    integer              :: nmax
  -    integer, allocatable :: wrk(:)
  -    real(float64)                :: n_empty_bins
  -    integer              :: n_multi
  +    integer :: err
   
       alphabet_size = nc
   
  -    ! compute multiplicities
  -    ! nmax is the largest number of samples in a bin
  -    nbins = size(counts)
  -    nmax = maxval(counts)
  -    allocate(wrk(nmax),stat=err)
  -    ! wrk(n) is the number of states with frequency n
  -    wrk = 0
  -    ! take into account the alphabet_size - nbins states with zero frequency
  -    n_empty_bins = alphabet_size - nbins
  -    do i_ = 1,nbins
  -       ni_ = counts(i_)
  -       if (ni_ == 0) then
  -          n_empty_bins = n_empty_bins + 1.0_float64
  -       else
  -          wrk(ni_) = wrk(ni_) + 1
  -       end if
  -    end do
  -
  -    ! further compress data into 'sparse' multiplicities
  -    n_multi = count(wrk > 0)
  -    allocate(hn(n_multi+1),stat=err)
  -    allocate(hz(n_multi+1),stat=err)
  -    hn(1) = 0
  -    hz(1) = n_empty_bins
  -    k_ = 1
  -    do i_ = 1, nmax
  -       if (wrk(i_) > 0) then
  -          k_ = k_ + 1
  -          hn(k_) = i_
  -          hz(k_) = wrk(i_)
  -       end if
  -    end do
  -    deallocate(wrk)
  -
  -    allocate(phi(n_multi+1), stat=err)
  -
  -    n_data = sum(hz * hn)
  +    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(hn1, hz1, nc)
  -    ! set hn, hz, n_data, alphabet_size
  -    real(float64), intent(in) :: hn1(:)
  -    real(float64), intent(in) :: hz1(:)
  +  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(hn1)
  +    nm = size(nk1)
       ! check if zeros are included in mults arrays
  -    if (any(int(hn1) == 0)) then  ! we work with floats
  -       allocate(hn(nm), stat=err)
  -       allocate(hz(nm), stat=err)
  +    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)
  -       hn = hn1
  -       hz = hz1
  +       nk = nk1
  +       zk = zk1
          idx = -1
          do j = 1, nm
  -          if (int(hn1(j)) == 0) then
  +          if (int(nk1(j)) == 0) then
                idx = j
                exit
             end if
          end do
  -       hz(idx) = hz1(idx) + nc - sum(hz1)
  +       zk(idx) = zk1(idx) + nc - sum(zk1)
       else
  -       allocate(hn(nm + 1), stat=err)
  -       allocate(hz(nm + 1), stat=err)
  +       allocate(nk(nm + 1), stat=err)
  +       allocate(zk(nm + 1), stat=err)
          allocate(phi(nm + 1), stat=err)
  -       hn(1) = 0._float64
  -       hz(1) = nc - sum(hz1)
  -       hn(2:) = hn1
  -       hz(2:) = hz1
  +       nk(1) = 0._float64
  +       zk(1) = nc - sum(zk1)
  +       nk(2:) = nk1
  +       zk(2:) = zk1
       end if
   
       alphabet_size = nc
  -    n_data = sum(hz * hn)
  +    n_data = sum(zk * nk)
   
     end subroutine initialize_from_multiplicities
   
     subroutine finalize()
   
  -    if (allocated(hz)) then
  -       deallocate(hz)
  -    end if
  -
  -    if (allocated(hn)) then
  -       deallocate(hn)
  +    if (allocated(zk)) then
  +       deallocate(zk)
  +    end if
  +
  +    if (allocated(nk)) then
  +       deallocate(nk)
       end if
   
       if (allocated(phi)) then
  
positive.code:     - alphabet_size * digamma(alpha) &
           - alphabet_size * digamma(n_data + alpha * alphabet_size)
  
      wsum = sum(hz * (digamma(hn + alpha)))
  
      dlpna = dlpna + wsum
  
  
negative.code:          - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size)
  
      wsum = sum(zk * (log_gamma(nk + alpha) - log_gamma(nk + one)))
  
      log_pna = log_pna + wsum
  
  