anchor.diff:      
       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))
  +    real(real64)   :: a(n_multi)
  +
  +    log_fpxa = log_gamma(n_data + one) + log_gamma(alpha * alphabet_size) & 
  +         - alphabet_size * log_gamma(alpha) - log_gamma(n_data + alpha * alphabet_size)
  +    do i = 1,n_multi 
  +       a(i) = multi(i) * (log_gamma(multi_z(i) - one + alpha) - log_gamma(multi_z(i) * one))
       end do
       log_fpxa = log_fpxa + sum(a)
       
  
positive.code:   real(real64),   intent(in)  :: alpha
    real(real64),   intent(out) :: estimate
  
    integer(int32) :: nbins,ndata
    integer(int32) :: i
    real(real64)   :: ni
  
  
negative.code:     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
  
      ! 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)
  