anchor.diff:      end do
       deallocate(multiplicities)
       
  -  end subroutine dirichlet_initialize
  +  end subroutine compute_multiplicities
   
     subroutine dirichlet_finalize()
   
  
positive.code: 
    !! this is the old 'standard' implementation.
    !! the code using multiplicities is faster
    !  nbins = size(hist)
    !  if (nbins == 1) then 
    !     estimate = 0.0_real64
    !     return
    !  end if
    !  ndata = sum(hist)*1.0_real64
    !  estimate = 0.0_real64
    !  do i = 1,nbins
    !     ni = hist(i)*1.0_real64
    !     if (ni > 0) estimate = estimate - ni*log(ni)
    !  end do
    !  estimate = estimate / ndata + log(ndata)
  
    ! compute multiplicities 
    ! nmax is the largest number of samples in a bin
    nbins = size(hist)
    if (nbins == 1) then 
       estimate = 0.0_real64
       return
    end if
    ndata = sum(hist)*1.0_real64
    nmax = maxval(hist)
    allocate(multiplicities(nmax),stat=err)
    multiplicities = 0
    do i = 1,nbins
       ni = hist(i)
       if (ni == 0) cycle
       multiplicities(ni) = multiplicities(ni) + 1
    end do
  
negative.code:     integer(int32), allocatable :: multiplicities(:)
  
      alphabet_size = nc 
      ndata = sum(counts)
  
      ! compute multiplicities 
      ! nmax is the largest number of samples in a bin
      nbins = size(counts)
      nmax = maxval(counts)
      allocate(multiplicities(nmax+1),stat=err)
      multiplicities = 0
      multiplicities(1) = alphabet_size - nbins
      do i = 1,nbins
         ni = counts(i)
         multiplicities(ni+1) = multiplicities(ni+1) + 1
      end do
  
  