anchor.diff:    real(real64)                :: alphabet_size
     real(real64), allocatable :: multi_z(:)  ! array of observed frequencies
     real(real64), allocatable :: multi(:)  ! multiplicities of frequency z
  -  real(real64), allocatable :: lgz1(:)  ! log_gamma(z + 1)
   
   contains
   
  
positive.code:     end do
      deallocate(multi0)
  
      ! compute log_gamma(z + 1) once for all
      allocate(lgz1(0:n_multi),stat=err)
      lgz1 = log_gamma(multi_z + one)
  
    end subroutine compute_multiplicities
  
    subroutine dirichlet_finalize()
  
      deallocate(multi_z,multi, lgz1)
  
    end subroutine dirichlet_finalize
  
  
negative.code:          - 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)))
  
      log_pna = log_pna + wsum
  
  