anchor.diff:         lw_max = log_weight(amax)
          integrand = exp(log_weight(alpha) - lw_max)  * alpha / amax
       else
  -       hb = h_dir(alpha)
  -
  -       integrand = hb**order
  +       if (order == 1) then
  +          integrand = h_dir(alpha)
  +       else if (order == 2) then
  +          integrand = h_var(alpha)
  +       end if
   
          lw_max = log_weight(amax)
   
  
positive.code:   real(real64)                :: alphabet_size
    real(real64), allocatable :: multi_z(:)  ! array of observed frequencies
    real(real64), allocatable :: multi(:)  ! multiplicities of frequency z
  
  contains
  
  
negative.code:     end do
      deallocate(multi0)
  
      allocate(phi(0:n_multi), stat=err)
  
    end subroutine compute_multiplicities
  
    subroutine finalize()
  