anchor.diff:      allocate(multi0(nmax+1),stat=err)
       multi0 = 0
       multi0(1) = alphabet_size - nbins
  -    do i = 1,nbins
  -       ni = counts(i)
  -       multi0(ni+1) = multi0(ni+1) + 1
  -    end do
  -
  -     ! working arrays on bins visited at leat once 
  -    nmulti = count(multi0 > 0) 
  -    allocate(multi_x(nmulti),stat=err)
  -    allocate(multi(nmulti),stat=err)
  -
  -    k = 0
  -    do i = 1,nmax+1
  -       if (multi0(i) > 0) then 
  -          k = k + 1
  -          multi_x(k) = i
  -          multi(k) = multi0(i)
  +    do i_ = 1,nbins
  +       ni_ = counts(i_)
  +       multi0(ni_ + 1) = multi0(ni_ + 1) + 1
  +    end do
  +
  +    ! further compress data into 'sparse' multiplicities
  +    n_multi = count(multi0 > 0) 
  +    allocate(multi_z(n_multi),stat=err)
  +    allocate(multi(n_multi),stat=err)
  +
  +    k_ = 0
  +    do i_ = 1, nmax + 1
  +       if (multi0(i_) > 0) then 
  +          k_ = k_ + 1
  +          multi_z(k_) = i_
  +          multi(k_) = multi0(i_)
          end if
       end do
       deallocate(multi0)
  
positive.code:     
      real(real64), intent(in) :: alpha
      integer(int32) :: i
      real(real64)   :: a(nmulti)
  
      log_fpxa = log_gamma(ndata + one) + log_gamma(alpha * alphabet_size) & 
           - alphabet_size * log_gamma(alpha) - log_gamma(ndata + alpha * alphabet_size)
      do i = 1,nmulti 
         a(i) = multi(i) * (log_gamma(multi_x(i) - one + alpha) - log_gamma(multi_x(i) * one))
      end do
      log_fpxa = log_fpxa + sum(a)
      
  
negative.code:      stop
    end if
    if (nc > nbins) estimate = estimate - (nc - nbins)*alpha*log(alpha)
    estimate = estimate / (n_data + nc*alpha) + log(n_data + nc*alpha)
  
  end subroutine pseudo
  
  