anchor.diff:  
       hdir = 0.0_real64
       do i = 1,n_multi 
  -       a(i) = - multi(i) * (multi_z(i) - one + alpha) * digamma(multi_z(i) + alpha) 
  +       a(i) = - multi(i) * (multi_z(i) + alpha) * digamma(multi_z(i) + one + alpha) 
       end do
       hdir = sum(a)
       hdir = hdir / (n_data + alpha * alphabet_size)
  
positive.code:     do i_ = 1, nmax + 1
         if (multi0(i_) > 0) then 
            k_ = k_ + 1
            multi_z(k_) = i_
            multi(k_) = multi0(i_)
         end if
      end do
  
negative.code:     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_) + alpha) - log_gamma(multi_z(i_) + one))
      end do
      log_fpxa = log_fpxa + sum(a)
      
  