Stochastically calculate staggered eigenvalue mode number nu(lambda)

Features:
  Arbitrary number of nHYP smearing steps with fixed coefficients
  Arbitrary number "npbp" of equally-spaced lambda at which to calculate nu,
    starting from a given startomega with spacing "spacing"
  Correct Omega --> Omega_* for given minmax polynomial
  Same stochastic sources used for all lambda

More details:
The step function is estimated with a minmax polynomial
Several different polynomial coefficients are hard-coded in coeffs.c
They are calculated by polynomial/approx_for_modenumber.c,
  as is the ratio (Omega --> Omega_*)^2
The function coefficients() in coeffs.c switches between the polynomials
  based on the input parameter "order"

The polynomial P(X) itself depends on the gauge field (obviously)
  P(X) psi = sum_{i = 0}^n c_i T_i(Z) psi, where
  X = 1 - 8Omega_*^2 / (Ddag D + Omega_*^2)
  Z = (2X - 1 - epsilon) / (1 - epsilon)
We apply it using the recursive Clenshaw algorithm
Each iterations requires that we invert Ddag D + Omega_*^2 to apply Z

NOTE: Due to MILC conventions for the normalization of the Dirac operator,
      X has a factor of 8 in the numerator instead of a factor of 2
      The input lambda is twice the Omega in the Dirac operator

The step function is then h(x) psi = (1 - X P(X^2) ) psi / 2
We calculate |h(X)^2 psi|^2 to get a more accurate result for the mode number

The coefficients c_i in P(X) are computed offline and copied into coeffs.c
The serial (GSL-based) code is in the polynomial directory
  https://www.gnu.org/software/gsl/
It finds the smallest Norder that satisfies the input epsilon and delta
It is based on code kindly shared by Agostino Patella
  http://inspirehep.net/author/profile/A.Patella.1

Sample infile:
+++
prompt 0
nx 4
ny 4
nz 4
nt 4
iseed 41

Nsmear 1
alpha_hyp0 0.5
alpha_hyp1 0.5
alpha_hyp2 0.4

npbp 5
order 32
npts 5
startomega 0
spacing 0.5

max_cg_iterations 5000
max_cg_restarts 1
error_per_site 1e-6

reload_serial config.4444