anchor.diff:      - alphabet_size * digamma(alpha) &
            - alphabet_size * digamma(n_data + alpha * alphabet_size)
   
  -    wsum = sum(multi * (digamma(multi_z + alpha)))
  +    wsum = sum(hzs * (digamma(hns + alpha)))
   
       dlpna = dlpna + wsum
   
  
positive.code: 
      ! further compress data into 'sparse' multiplicities
      n_multi = count(multi0 > 0)
      allocate(multi_z(n_multi+1),stat=err)
      allocate(multi(n_multi+1),stat=err)
      multi_z(1) = 0
      multi(1) = n_empty_bins
      k_ = 1
      do i_ = 1, nmax
         if (multi0(i_) > 0) then
            k_ = k_ + 1
            multi_z(k_) = i_
            multi(k_) = multi0(i_)
         end if
      end do
      deallocate(multi0)
  
      allocate(phi(n_multi+1), stat=err)
  
      n_data = sum(multi * multi_z)
  
    end subroutine initialize_from_counts
  
    subroutine finalize()
  
      if (allocated(multi)) then
         deallocate(multi)
      end if
  
      if (allocated(multi_z)) then
         deallocate(multi_z)
      end if
  
      if (allocated(phi)) then
  
negative.code: 
    integer(int32)              :: n_data
    real(real64)                :: alphabet_size
    real(real64), allocatable :: hns(:)  ! array of observed frequencies
    real(real64), allocatable :: hzs(:)  ! multiplicities of frequency z
    real(real64), allocatable :: phi(:)  ! wrk array for var
  
  contains
  
    subroutine initialize_from_counts(counts, nc)
      ! set n_multi, hns, multi
      use constants
      integer(int32), intent(in) :: counts(:)
      real(real64), intent(in) :: nc
  