anchor.diff:       end if
        n_data = sum(counts)*1.0d0
        nmax = maxval(counts)
  -     allocate(multi0(nmax),stat=err)
  -     multi0 = 0
  +     allocate(wrk(nmax),stat=err)
  +     wrk = 0
        do i = 1,nbins
           ni = counts(i)
           if (ni == 0) cycle
  -        multi0(ni) = multi0(ni) + 1
  +        wrk(ni) = wrk(ni) + 1
        end do
        estimate = 0.0d0
        do i = 1,nmax
  -        mi = multi0(i)
  +        mi = wrk(i)
           if (mi > 0) estimate = estimate - mi*i*log(i*1.0d0)
        end do
        estimate = estimate / n_data + log(n_data)
  -     deallocate(multi0)
  +     deallocate(wrk)
     else
        ! standard implementation
        nbins = size(counts)
  
positive.code:     integer              :: i_,k_,ni_
      integer              :: err
      integer              :: nmax
      integer, allocatable :: multi0(:)
      real(8)                :: n_empty_bins
      integer              :: n_multi
  
  
negative.code:     ! 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
  