anchor.diff:       return
     end if
   
  +  call initialize_dirichlet(counts, nc)
     call compute_multiplicities(counts, nc)
   
     estimate = hdir(alpha)
  
positive.code: subroutine dirichlet(n,counts,nc,alpha,estimate)
    ! posterior mean entropy (averaged over Dirichlet distribution) given alpha 
    use iso_fortran_env
    use dirichlet_mod, only: compute_multiplicities,dirichlet_finalize
    use dirichlet_mod, only: hdir
    implicit none
  
  
negative.code:      return
    end if
  
    call initialize_dirichlet(counts, nc)
    call compute_multiplicities(counts, nc)
  
    call compute_integration_range()
  