anchor.diff:         err = err + ierr
          err_estimate = err_estimate / nrm
          err_estimate = sqrt(err_estimate - estimate**2)
  +       if (isnan(err_estimate)) then
  +          err_estimate = 0.0_real64
  +       end if
       end if
   
     end subroutine hnsb
  
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:        lw_max = log_weight(amax)
         integrand = exp(log_weight(alpha) - lw_max)  * alpha / amax
      else
         if (order == 1) then
            integrand = h_dir(alpha)
         else if (order == 2) then
            integrand = h_var(alpha)
         end if
  
         lw_max = log_weight(amax)
  
  