anchor.diff:      integer              :: i_,k_,ni_
       integer              :: err
       integer              :: nmax
  -    integer, allocatable :: multi0(:)
  +    integer, allocatable :: wrk(:)
       real(8)                :: n_empty_bins
       integer              :: n_multi
   
  
positive.code:      end if
       n_data = sum(counts)*1.0d0
       nmax = maxval(counts)
       allocate(multi0(nmax),stat=err)
       multi0 = 0
       do i = 1,nbins
          ni = counts(i)
          if (ni == 0) cycle
          multi0(ni) = multi0(ni) + 1
       end do
       estimate = 0.0d0
       do i = 1,nmax
          mi = multi0(i)
          if (mi > 0) estimate = estimate - mi*i*log(i*1.0d0)
       end do
       estimate = estimate / n_data + log(n_data)
       deallocate(multi0)
    else
       ! standard implementation
       nbins = size(counts)
  
negative.code:   real(8), intent(out) :: y
    y = trigamma(x)
  end subroutine gamma1
  
  module counts
    implicit none
    integer, allocatable :: nk(:)
    integer, allocatable :: zk(:)
    integer :: ndata
    integer :: nbins
    integer :: k1
  contains
    subroutine fit(n, counts)
      implicit none
      integer, intent(in) :: n
      integer, intent(in) :: counts(n)
      integer :: i, j, u
      integer :: x, err
      integer :: xmax
      integer, allocatable :: wrk(:)
  
      xmax = maxval(counts)
      allocate(wrk(0:xmax), stat=err)
      wrk = 0
      u = 0
      do i = 1,n
         x = counts(i)
         if (wrk(x) == 0) then
            u = u + 1
         end if
         wrk(x) = wrk(x) + 1
      end do
  
      if (allocated(nk)) then
         deallocate(nk)
      end if
      if (allocated(zk)) then
         deallocate(zk)
      end if
      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
  
      ndata = sum(zk * nk)
      nbins = sum(zk)
      k1 = sum(zk, nk>0)
  
      deallocate(wrk)
    end subroutine fit
  end module counts
  
  