anchor.diff:    end if
     ndata = sum(counts)*1.0_real64
     nmax = maxval(counts)
  -  allocate(multip(nmax),stat=err)
  -  multip = 0
  +  allocate(multiplicities(nmax),stat=err)
  +  multiplicities = 0
     do i = 1,nbins
        ni = counts(i)
        if (ni == 0) cycle
  -     multip(ni) = multip(ni) + 1
  +     multiplicities(ni) = multiplicities(ni) + 1
     end do
     estimate = 0.0_real64
     do i = 1,nmax
  -     mi = multip(i)
  +     mi = multiplicities(i)
        if (mi > 0) estimate = estimate - mi*i*log(i*1.0_real64)
     end do
     estimate = estimate / ndata + log(ndata)
  -  deallocate(multip)
  +  deallocate(multiplicities)
   
   end subroutine plugin
   
  
positive.code:   integer(int32), intent(in)  :: nc
    real(real64),   intent(in)  :: alpha
    real(real64),   intent(out) :: estimate
    integer(int32), allocatable :: sparse_multip_frq(:)
    integer(int32), allocatable :: sparse_multip(:)
  
    if (size(counts) == 1) then 
       estimate = 0.0_real64
       return
    end if
  
    call compute_multiplicities(counts, nc, sparse_multip_frq, sparse_multip)
  
    estimate = hdir(alpha)
  
    deallocate(sparse_multip_frq,sparse_multip)
  
  end subroutine dirichlet
  
  subroutine nsb(n,counts,nc,estimate,err_estimate)
    use iso_fortran_env
    use dirichlet_mod, only: compute_multiplicities
    use nsb_mod, only: hnsb
    use nsb_mod, only: compute_integration_range
    implicit none
  
negative.code:     
      real(real64), intent(in) :: alpha
      integer(int32) :: i
      real(real64)   :: a(ngtz)
  
      hdir = 0.0_real64
      do i = 1,ngtz 
         a(i) = - multi_gtz(i) * (fgtz(i) - one + alpha) * digamma(fgtz(i) + alpha) 
      end do
      hdir = sum(a)
      hdir = hdir / (ndata + alpha * alphabet_size)
  