anchor.diff:      log_fpxa = log_gamma(n_data + one) + log_gamma(alpha * alphabet_size) & 
            - alphabet_size * log_gamma(alpha) - log_gamma(n_data + alpha * alphabet_size)
       do i_ = 1, n_multi 
  -       a(i_) = multi(i_) * (log_gamma(multi_z(i_) - one + alpha) - log_gamma(multi_z(i_) * one))
  +       a(i_) = multi(i_) * (log_gamma(multi_z(i_) + alpha) - log_gamma(multi_z(i_) + one))
       end do
       log_fpxa = log_fpxa + sum(a)
       
  
positive.code: 
      hdir = 0.0_real64
      do i = 1,n_multi 
         a(i) = - multi(i) * (multi_z(i) - one + alpha) * digamma(multi_z(i) + alpha) 
      end do
      hdir = sum(a)
      hdir = hdir / (n_data + alpha * alphabet_size)
  
negative.code:   end if
  
    call initialize_dirichlet(counts, nc)
    call compute_multiplicities(counts)
  
    call compute_integration_range()
  
  