anchor.diff:           - alphabet_size * log_gamma(alpha) &
            - log_gamma(n_data + alpha * alphabet_size)
   
  -    wsum = sum(multi * (log_gamma(multi_z + alpha) - log_gamma(multi_z + one)))
  +    wsum = sum(hzs * (log_gamma(hns + alpha) - log_gamma(hns + one)))
   
       log_pna = log_pna + wsum
   
  
positive.code:     log_pna_u = log_gamma(alpha * alphabet_size) &
           - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size) &
           + sum(multi * (log_gamma(multi_z + alpha)))
  
    end function log_pna_u
  
  
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
  