anchor.diff:  
       ! 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
  +    allocate(hns(n_multi+1),stat=err)
  +    allocate(hzs(n_multi+1),stat=err)
  +    hns(1) = 0
  +    hzs(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_)
  +          hns(k_) = i_
  +          hzs(k_) = multi0(i_)
          end if
       end do
       deallocate(multi0)
   
       allocate(phi(n_multi+1), stat=err)
   
  -    n_data = sum(multi * multi_z)
  +    n_data = sum(hzs * hns)
   
     end subroutine initialize_from_counts
   
     subroutine finalize()
   
  -    if (allocated(multi)) then
  -       deallocate(multi)
  +    if (allocated(hzs)) then
  +       deallocate(hzs)
       end if
   
  -    if (allocated(multi_z)) then
  -       deallocate(multi_z)
  +    if (allocated(hns)) then
  +       deallocate(hns)
       end if
   
       if (allocated(phi)) then
  
positive.code:     real(real64) :: c, nu, ni, xi, jsum
  
      nu = n_data + alpha * alphabet_size
      phi = digamma(multi_z + alpha + one) - &
           digamma(nu + two)
      c = trigamma(nu + two)
  
      h_var = 0.0
      do i_ = 0, size(multi)-1
         ni = multi_z(i_) + alpha
         xi = phi(i_)
         jsum = sum(multi * ni * (multi_z + alpha) * &
              (xi * phi - c))
         h_var = h_var + multi(i_) * jsum
         h_var = h_var - multi(i_) * ni**2 * (xi**2 - c)
         xi = xi + 1 / (ni + one)
         h_var = h_var + multi(i_) * (ni + one) * ni * &
              (xi**2 + trigamma(ni + two) - c)
      end do
  
  
negative.code:     - alphabet_size * digamma(alpha) &
           - alphabet_size * digamma(n_data + alpha * alphabet_size)
  
      wsum = sum(hzs * (digamma(hns + alpha)))
  
      dlpna = dlpna + wsum
  
  