%\documentclass[11pt]{cernrep}
%\usepackage{cite}
%\usepackage{graphicx}
%
%\begin{document}
%\bibliographystyle{unsrt}
%
\section{PARTON DISTRIBUTION FUNCTIONS\protect\footnote{Section
coordinators: A.~Vogt, W.~Giele}$^{,}$~\protect\footnote{Contributing authors:
S.~Alekhin, M.~Botje, W.~Giele, J.~Pumplin, F.~Olness, G.~Salam, A.~Vogt}}
\label{sec:pdf,qcdsm}

The experimental uncertainties in current and future hadronic colliders are 
decreasing to a level where more careful consideration has to be given to the 
uncertainties in the theoretical predictions. One important source of these 
uncertainties has its origin in Parton Distribution Functions (PDFs). The PDF 
uncertainties in turn are a reflection of those in the experimental data used 
as an input to the PDF fits and in the uncertainties of the theoretical 
calculations used to fit those data. As a consequence, sophisticated 
statistical methods and approximations have to be developed to handle the 
propagation of uncertainties. 
We will give a summary of the current status of several methods being pursued. 
To fully exploit these new types of PDF fits a uniform computer interface has 
been defined and developed. This code provides easy access to all the present 
and future fits. The code is available from the website {\bf\tt pdf.fnal.gov}. 
Such an interface is most efficient if the evolution codes of the various 
groups fiiting the PDFs are standardized to a sufficient degree. For this 
purpose detailed and accurate reference tables for the PDF evolution are 
provided at the end of this section. 

\subsection{Methods for estimating parton distribution function uncertainties
\protect\footnote{Contributing authors:
S.~Alekhin, M.~Botje, W.~Giele, J.~Pumplin, F.~Olness}}

\subsubsection{A mathematical framework}

At first sight PDF fitting is rather straightforward. However, a more
detailed look reveals many difficult issues. As the PDF uncertainties
will affect all areas of phenomenology at hadron colliders, 
a clear mathematical framework of a PDF fit is essential \cite{Giele:2001mr}. 
From this
formulation, all explicit methods can be derived. Also, the mathematical
description will make explicit all assumptions needed before one can 
make a fit. These assumptions are crucial and do not find their origin
in experimental results,  but rather in theoretical prejudice. Such 
assumptions are unavoidable as we fit a system with an infinite 
number of degrees of freedom (the PDF functional) 
to a finite set of data points.

We want to calculate the probability density function $P^{\cal O}_{pdf}$
which reflects the uncertainty in predicting the observable $\cal O$
due to the PDF uncertainties. The function $P^{\cal O}_{pdf}(x_e)$ gives 
the probability density to measure a value $x_e$ for observable $\cal O$.

To calculate the PDF probability density function for observable $\cal O$
we have to integrate over the functional space of all possible PDFs $V({\cal F})$.
The integration is weighted by three probability density functions:
the prior probability density function, $P_{prior}$, 
the experimental response function of the observable, $P_{exp}^{\cal O}$
and the probability density function of the fitted experiments, 
$P_{exp}^{input}$. The resulting formula is given by
\begin{equation}\label{eqn:pdf;qcdsm;Ifunctional}
P_{pdf}^{\cal O}(x_e) = \int_{V({\cal F})} d\,{\cal F}\ P_{prior}({\cal F})\times 
P_{exp}^{input}({\cal F})\times P_{exp}^{\cal O}\left(x_e|x_t({\cal F})\right)\ .
\end{equation}

The prior probability density function contains theoretical constraints
on the PDF functionals such as sum rules and other potential
fit constraints (e.g. an explicit $\alpha_S(M_Z)$ value).
The most crucial property of the prior function is that it  
defines the functional integral by imposing smoothness constraints
to make the number of degrees of freedom become finite. The simplest example
is an explicit finite parametrization of the PDF functionals
\begin{equation}
P_{pdf}^{\cal O}(x_e) = \int_{V(\{\lambda\})}\!\!\!\!\! 
d\,\lambda_1\,d\,\lambda_2\ldots\,d\,\lambda_n\ 
P_{prior}(\{\lambda\})\times P_{exp}^{input}(\{\lambda\})\times 
P_{exp}^{\cal O}\left(x_e|x_t(\{\lambda\})\right)\ ,
\end{equation}
where the PDF parameters are given by the list $\{\lambda\}$. Note that
through the functional parametrization choice we have restricted
the integration to a specific subset of potential PDF functionals \mbox{$\cal F$}.
Such a choice is founded on physics assumptions with no {\it a priori}
justification. The resulting phenomenology depends on this assumption.

The experimental response function forms the interface with the experimental
measurement. It gives the probability density to
measure a value $x_e$ given a prediction $x_t$.
The experimental response functions contain all information about
the experiment needed to analyze the measurement.
The prediction $x_t$ is an 
approximation using a perturbative estimate of the observable
amended with relevant nonperturbative corrections.
A simple example would be
\[ P_{exp}^{\cal O}(x_e|x_t ) =\frac{1}{\delta\sqrt{2\pi}}
\exp{\left(-\frac{1}{2}(x_e-x_t)^2/\delta^2\right)}\ ,\]
where we have a Gaussian response function 
with a one sigma uncertainty $\delta$.
Note that the form of the response function depends on the actual
experiment under consideration. It is sometimes convenient to
get a result that is independent of an experiment and its particular
detector: To obtain the
theoretical prediction for the probability density function of the
observable, one can simply replace
the experimental response function by a delta function
(i.e. assume a perfect detector)
\[P_{exp}^{theory}(x_e|x_t) = \delta(x_e-x_t)\ .\]

Finally the probability function for the fitted experiments is simply
a product of all the experimental response functions
\begin{equation}
P_{exp}^{input}({\cal F}) = 
\prod_i^{N_{exp}} P_{exp}^{(i)}(x_m^{(i)}|x_t^{(i)}({\cal F})) \ ,
\end{equation} 
where $x_m^{(i)}$ denotes the set of measurements provided by
experiment $(i)$.
This function measures the probability density that
the theory prediction based on PDF \mbox{$\cal F$}\ 
describes the combined experimental measurements.

Often the experimental uncertainties can be approximated by a Gaussian
form of the experimental reponse function, i.e. a $\chi^2$ 
description of the uncertainties:
\begin{equation}
P_{pdf}^{\cal O}(x_e)\propto \int_{V(\{\lambda\})}\!\!\!\!\! 
d\,\lambda_1\,d\,\lambda_2\ldots\,d\,\lambda_n\ 
e^{-\frac{1}{2}\sum_i \chi_i^2(x_m^{(i)}-x_t^{(i)}(\{\lambda\}))}\times 
e^{-\frac{1}{2}\chi_{\cal O}^2(x_e^{(i)}-x_t(\{\lambda\}))}\ ,
\end{equation}
where we have chosen a specific parametrization for the PDF functionals.
This approximation leads to a more traditional approach for determining
PDFs with uncertainties. These methods are outlined in sections
\ref{sec:pdf;qcdsm;covariance} and \ref{sec:pdf;qcdsm;offset}.
To go beyond the Gaussian approximation more elaborate methods are needed.
Sections \ref{sec:pdf;qcdsm;ransam} and \ref{sec:pdf;qcdsm;lagrange} 
describe techniques to accomplish this.

\subsubsection{Random sampling method}\label{sec:pdf;qcdsm;ransam}

This method attempts to calculate Eq. (\ref{eqn:pdf;qcdsm;Ifunctional}) without
any approximations by using random sampling, i.e. a Monte Carlo evaluation of 
the integral \cite{Giele:2001mr}. By generating a sufficiently 
large sample of PDFs Eq. (\ref{eqn:pdf;qcdsm;Ifunctional})
can be approximated by
\begin{equation}\label{eqn:pdf;qcdsm;Irandomsampling}
P_{pdf}^{\cal O}(x_e)\approx\frac{1}{N}
\sum_{i=1}^N P_{prior}({\cal F}_i)\times P_{exp}^{input}({\cal F}_i)\times  
P_{exp}^{\cal O}\left(x_e|x_t({\cal F}_i)\right)\ .
\end{equation}

The simple implementation of Eq. (\ref{eqn:pdf;qcdsm;Irandomsampling})
would lead to a highly inefficient
random sampling and an unreasonably large number of PDFs would be
required. By using an optimization procedure, this problem can be
solved. The optimization procedure on the functional level of
Eq. (\ref{eqn:pdf;qcdsm;Ifunctional}) is simply redefining the PDF \mbox{$\cal F$}\ to ${\cal F}^u$
such that the Jacobian of the transformation equals
\begin{equation}
\label{eqn:pdf;qcdsm;transform}
\frac{d\,{\cal F}^u}{d\,{\cal F}} = P_{prior}({\cal F})\times P_{exp}^{input}({\cal F})\ ,
\end{equation}
so that
\begin{equation}
P_{pdf}^{\cal O}(x_e) = 
\int_{V({\cal F}^u)}d\,{\cal F}^u\ P_{exp}^{\cal O}\left(x_e|x_t({\cal F}^u)\right)\ .
\end{equation}
Applying this to the random sampling evaluation gives
\begin{equation}\label{eqn:pdf;qcdsm;density}
P_{pdf}^{\cal O}(x_e)\approx\frac{1}{N}\sum_{i=1}^N 
P_{exp}^{\cal O}\left(x_e|x_t({\cal F}_i^u)\right)\ .
\end{equation}
In the random sampling approximation the redefined PDFs ${\cal F}_i^u$ 
are easily identified as the unweighted PDFs with respect to
the combined probability density 
$P_{prior}({\cal F})\times P_{exp}^{input}({\cal F})$. That is, the density of
${\cal F}_i^u$ is given by this combined probability density. As such,
each of the unweighted PDFs is equally likely.  
This is reflected in Eq. (\ref{eqn:pdf;qcdsm;density}), 
as the probability density function
of the observable is the average of the response function over
the unweighted PDFs.
Finally, we have to generate the set $\{{\cal F}_i^u\}$. 
One method is to use the Gaussian approximation for 
Eq. (\ref{eqn:pdf;qcdsm;transform}) simplifying the 
generation of the set $\{{\cal F}_i^u\}$ \cite{Giele:1998gw}. 

Another more general approach is to apply a Metropolis Algorithm on the
combined probability density function 
$P_{prior}({\cal F})\times P_{exp}^{input}({\cal F})$.
This approach will handle any probability function. Furthermore, in a
Metropolis Algorithm approach convergence does not depend on the number
of parameters used in the PDF parametrization. Also, complicated 
non-linear parameter subspaces which have constant probability will
be modelled correctly. These properties make it possible to use
large number of parameters and explore the issue of parametrization
dependence of the PDFs.

Once the set $\{{\cal F}_i^u\}$ is generated we can predict the probability
density function for any observable by averaging over the experimental
response function using Eq. (\ref{eqn:pdf;qcdsm;density}).

\subsubsection{Lagrange multiplier method}\label{sec:pdf;qcdsm;lagrange}

The values of the fit parameters that minimize $\chi^{2}$ provide by 
definition the best fit to the global set of data.  The dependence of 
$\chi^{2}$ on those parameters in the neighborhood 
of the minimum can be characterized in quadratic approximation by the 
matrix of second derivatives, which is known as the Hessian matrix.
The inverse of the Hessian matrix is the error matrix, and it 
forms the basis for traditional estimates of the uncertainties.

The traditional assumption of quadratic behavior of $\chi^{2}$ in all 
directions in parameter space is not necessarily a very good one in 
the case of PDF fitting, because there are ``flat directions'' in which 
$\chi^{2}$ changes very slowly, so large changes in certain combinations 
of parameters are not ruled out.  This difficulty is always present in 
the case of PDF fitting, because as more data become available to pin 
down the PDFs better, more flexible parametrizations of them are 
introduced to allow ever-finer details to be determined.

To some extent, the flat directions can be allowed for by an 
{\it iterative method} \cite{Pumplin:2000vx,Pumplin:2001ct}, whereby the 
eigenvector directions of the error matrix and their eigenvalues 
are found using a convergent series of successive approximations.
This iterative method is implemented as an extension to the standard 
minimization program {\tt minuit}.
The result is a collection of PDFs (currently 40 in number) that 
probe the allowed range of possibilities along each of the 
eigenvector directions.  The PDF uncertainty on any physical 
quantity---or on some feature of the PDFs themselves---can easily 
be computed after that quantity is evaluated for each of the 
eigenvector sets.  This method has been applied, for example, 
to find the uncertainty range for predictions of $W$ and $Z$ cross 
sections and their correlations \cite{Pumplin:2001ct}.

To completely overcome the need to assume quadratic behavior for 
$\chi^{2}$ as a function of the fitting parameters, one can use a 
{\it Lagrange Multiplier method} \cite{Pumplin:2000vx,Stump:2001gu} to directly 
examine how $\chi^{2}$ varies as a function of any particular variable 
that is of interest.  
The method is a classical mathematical technique that is more fully 
called Lagrange's ``method of undetermined multipliers.''  To explain 
it by example, suppose we want to find the effect of PDF uncertainties 
on the prediction for the Higgs boson production cross section.  
Instead of finding the PDF parameters that minimize $\chi^{2}$ (which 
measures the quality of fit to the global data), one minimizes 
$\chi^{2} + \lambda \sigma_H$, where $\sigma_H$ is the predicted 
Higgs cross section---which is of course also a function of the PDF
parameters.  The minimization is carried out for a variety of values 
of the Lagrange Multiplier constant $\lambda$.
Each minimization provides one point on the curve of $\chi^{2}$ versus 
predicted $\sigma_H$.  Once that curve has been mapped out, the 
uncertainty range is defined by the region for which the increase in 
$\chi^{2}$ above its minimum value is acceptable.

The essential feature of the Lagrange Multiplier method is that it 
finds the largest possible range for the predictions of a given 
physical quantity, such as $\sigma_H$, that are consistent with 
any given assumed maximum increase $\Delta\chi^{2}$ above the 
best-fit value of $\chi^{2}$.  This method has been applied, 
for example, to study the possible range of variation of the 
rapidity distribution for $W$ production, by extremizing various 
moments of that distribution \cite{Pumplin:2000vx,Stump:2001gu}.

\subsubsection{Covariance matrix method}\label{sec:pdf;qcdsm;covariance}

The covariance matrix method  
is based on the Bayesian approach to the treatment of the systematic errors, 
when the latter are considered as a random fluctuation like the 
statistical errors. 
Having no place for the detailed discussion of the advantages of this 
approach we refer to the 
introduction into this scope given in Ref.~\cite{D'Agostini:1995fv};
the only point we would like to underline here is 
that application of the Bayesian approach is especially justified 
in the analysis of the data sets with the numerous
sources of independent systematic errors, which is the case for the 
extraction of PDFs from existing experimental data.

Let the data sample have one common
additive systematic error\footnote{We consider the case of one source 
of systematic errors, generalization
on the many sources case is straightforward.}.
In this case following the Bayesian approach the measured values 
$y_i$ are given by   
\begin{equation}
y_i=f_i+\mu_i \sigma_i+\lambda s_i,
\label{eqn:pdf;qcdsm;CADDSET}
\end{equation}
where $f_i(\theta^0)$ is the theoretical model describing the data,
$\theta^0$  -- the fitted parameter of this model,
$\sigma_i$ -- statistical errors, 
$s_i$  -- systematic errors for each point, $\mu_i$ and $\lambda$
-- the independent random variables, and the index $i$ runs 
through the data points from 1 to $N$.
The only assumption we make
is that the average and the dispersion of these variables are zero and unity 
respectively. 
It is natural to assume that the $\mu_i$ are Gaussian distributed when the
data points are obtained from large statistical samples. Such an
assumption is often not justified for the distribution of $\lambda$

Within the covariance matrix approach the estimate of the 
fitted parameter $\hat\theta$ is obtained 
from a minimization of the functional 
\begin{equation}
\chi^2(\theta)=\sum_{i,j=1}^{N} (f_i(\theta)-y_i) E_{ij} (f_j(\theta)-y_j),
\label{eqn:pdf;qcdsm;CORCHI}
\end{equation}
where 
\begin{displaymath}
E_{ij}=\frac{1}{\sigma_i \sigma_j}
\Bigl(\delta_{ij} -
\frac{\rho_i\rho_j}{1+\rho^2}\Bigr)\ ,
\end{displaymath}
is the inverse of the covariance matrix  
\begin{equation}
C_{ij} =  s_i s_j+\delta_{ij}\sigma_i\sigma_j,
\label{eqn:pdf;qcdsm;COVA}
\end{equation}
$\rho$ is modulus of the vector $\vec\rho$ with $N$ components equal to 
$\rho_i={s_i}/{\sigma_i}$ 
and $\delta_{ij}$ is the Kronecker symbol. 
We underline that with the data model given by Eq. (\ref{eqn:pdf;qcdsm;CADDSET})
one does not need to assume a specific form for the distribution of $y_i$
in order to obtain the central value of this estimate.
In a linear approximation for $f_i(\theta)$ one also does not need such 
assumptions to estimate the dispersion.
In this approximation the dispersion reads \cite{Alekhin:2000es}
\begin{equation}
D_{\rm C}(\hat\theta)=\frac{1}{\phi^2}
\left[1+\frac{\rho^2 z^2}{1+\rho^2(1-z^2)}\right],
\label{eqn:pdf;qcdsm;DISPCA}
\end{equation}
where $\phi$ is modulus of the vector $\vec\phi$ 
with $N$ components equal to 
$\phi^i={f_i'(\theta_0)}/{\sigma_i}$, the symbol $'$ denotes the derivative
with respect to $\theta$, and 
$z$ is the cosine of the angle between $\vec\rho$ and $\vec\phi$.

To obtain the distribution of $\hat\theta$ one needs to know 
the complete set of its moments, which, in turn, requires similar information 
for the moments of $y_i$. At the same time 
from considerations similar to the proof of 
the Central Limit theorem of statistics (see. Ref.~\cite{James:1971})
one can conclude that the distribution of $\hat\theta$ is Gaussian, 
if all independent systematic errors are comparable in value and 
their number is large enough. More educated guesses of the form of 
the distribution can be performed with the help of the general approach 
described in subsection \ref{sec:pdf;qcdsm;ransam}.

The dispersion of fitted parameters obtained by the 
covariance method is 
different from the one obtained by the offset method described 
in subsection \ref{sec:pdf;qcdsm;offset}. Indeed, the dispersion of the 
parameter estimate obtained by the offset method
applied to the analysis of the data set given by 
Eq. (\ref{eqn:pdf;qcdsm;CADDSET})
is equal \cite{Alekhin:2000es}
\begin{equation}
D_{\rm O}(\hat\theta)=\frac{1}{\phi^2}\left(1+\rho^2 z^2\right).
\label{eqn:pdf;qcdsm;DISPSA}
\end{equation}
One can see that $D_{\rm O}$ is generally larger than $D_{\rm C}$.
The difference is especially visible for the case $N\gg 1$, when 
$\rho \gg 1$, if the systematic errors are not negligible as compared 
to the statistical ones. In this case and if $z\ne 1$
\begin{equation}
D_{\rm C}(\hat\theta)\approx\frac{1}{\phi^2(1-z^2)}\ ,
\label{eqn:pdf;qcdsm;dispmr}
\end{equation}
while
\begin{equation}
D_{\rm O}(\hat\theta)\approx\frac{\rho^2 z^2}{\phi^2}.
\label{eqn:pdf;qcdsm;disp0r}
\end{equation}
One can see that the standard deviation of the offset method
estimator rises linearly with
the increase of the systematics, while the 
covariance matrix dispersion saturates and the difference between them 
may be very large.

Some peculiarities arise in the case when the systematic errors 
are multiplicative, i.e. when the $s_i$ in
Eq. (\ref{eqn:pdf;qcdsm;COVA})
are given by $\eta_i y_i$, where $\eta_i$ are constants. 
As it was noted in Ref.~\cite{D'Agostini:1994uj} in this case 
the covariance matrix estimator may be biased.
The manifestation of this bias is that 
the fitted curve lays lower the data points on average, 
which is reflected by a distortion of the fitted parameters.
In order to minimize such bias one has to calculate
the covariance matrix of Eq. (\ref{eqn:pdf;qcdsm;COVA})  using the 
relation $s_i=\eta_i f_i(\theta)$. In this approach
the covariance matrix depends on the fitted parameters and hence 
has to be iteratively re-calculated during the fit.
This certainly makes calculation more time-consuming and 
difficult, but in this case the bias of the estimator is non-negligible 
as compared to the value of its standard deviation if the systematic
error on the fitted parameter is an order of magnitude
larger than the statistical error \cite{Alekhin:2000es}.
\subsubsection{Offset Method}\label{sec:pdf;qcdsm;offset}
With the offset method~\cite{Botje:2001fx}, already mentioned above,
the systematic errors are incorporated in the model prediction
%
\begin{equation}\label{eqn:pdf;qcdsm;modelpred}
 t_i(\theta,\lambda) = f_i(\theta) + \sum_k \lambda_k\; s_{ik}\ ,
\end{equation}
%
where we allow for several sources of systematic error $(k)$.  The
$\chi^2$, to be minimized in a fit, is defined as
%
\begin{equation}\label{eqn:pdf;qcdsm;xpsdef}
  \chi^2(\theta,\lambda)
   = \sum_i \left( \frac{y_i-t_i(\theta,\lambda)}{\sigma_i} \right)^2 
  + \sum_k \lambda_k^2.   
\end{equation}
%
It can be shown~\cite{Stump:2001gu} that leaving both $\theta$ and
$\lambda$ free in the fit is mathematically equivalent to the
covariance method described in the previous section.  However, there
is also the choice to fix the systematic parameters to their central
values $\lambda_k = 0$ which results in minimizing
%
\begin{equation}\label{eqn:pdf;qcdsm;offsetchi}
 \chi^2(\theta) = \sum_i \left( \frac{y_i-f_i(\theta)}{\sigma_i} \right)^2\ ,
\end{equation}
%
where only statistical errors are taken into account to get the best
value $\hat{\theta}$ of the parameters.  Because systematic errors are
ignored in the $\chi^2$ such a fit forces the theory prediction to be
as close as possible to the data.

The systematic errors on $\theta$ are estimated from fits where each
systematic parameter $\lambda_k$ is offset by its assumed error
($\pm1$) after which the resulting deviations $\Delta \theta$ are
added in quadrature.  To first order this lengthy procedure can be
replaced by a calculation of two Hessian matrices $M$ and $C$, defined
by
%
\begin{equation}\label{eqn:pdf;qcdsm;mcijdef}
  M_{ij} = \frac{1}{2} \frac{\partial^2 \chi^2}
  {\partial \theta_i \partial \theta_j} \qquad
  C_{ij} = \frac{1}{2} \frac{\partial^2 \chi^2}
  {\partial \theta_i \partial \lambda_j}\ .
\end{equation}
%
The statistical covariance matrix of the fitted parameters is then
given by
%
\begin{equation}\label{eqn:pdf;qcdsm;cstatdef}
  V_{\rm stat}^{\theta} = M^{-1}\ ,
\end{equation}
%
while a systematic covariance matrix can be defined
by~\cite{Pascaud:1995qs}
%
\begin{equation}\label{eqn:pdf;qcdsm;csysdef}
  V_{\rm syst}^{\theta} = M^{-1} C C^T M^{-1}\ ,
\end{equation}
%
where $C^T$ is the transpose of $C$.

Having obtained the best values and the covariance matrix of the
parameters, the covariance of any two functions $F(\theta)$ and
$G(\theta)$ can be calculated with the standard formula for linear
error propagation
%
\begin{equation}\label{eqn:pdf;qcdsm;fgmat}  
  V_{FG} = \sum_{ij} \frac{\partial F(\theta)}{\partial \theta_i}\;
  V_{ij}^{\theta}\;
  \frac{\partial G(\theta)}{\partial \theta_j}\ ,
\end{equation}
%
where $V^{\theta}$ is the statistical, systematic or, if the total
error is to be calculated, the sum of both covariance matrices.

Comparing Eqs. (\ref{eqn:pdf;qcdsm;xpsdef})  
and (\ref{eqn:pdf;qcdsm;offsetchi}) it is
clear that the parameter values obtained by the covariance and offset
methods will, in general, be different. This difference is accounted
for by the difference in the error estimates, those of the offset
method being larger in most cases. In statistical language this means
that the parameter estimation of the offset method is not {\em
efficient}. The method has a further disadvantage that the goodness of
fit cannot be directly judged from the $\chi^2$ which is calculated
from statistical errors only.

For a global QCD analysis of deep inelastic scattering data which uses
the offset method to propagate the systematic errors, we refer
to~\cite{Botje:1999dj} (see {\tt http://www.nikhef.nl/user/h24/qcdnum}
for the corresponding PDF set with full error information).

\subsection{The LHAPDF interface
\protect\footnote{Contributing authors:
S.~Alekhin, W.~Giele, J.~Pumplin}}

\subsubsection{Introduction}

The Les Houches Accord PDF (LHAPDF) interface package is designed to work with
PDF sets. A PDF set can consist of many individual member PDFs. While
the interpretation of the member PDFs depends on the particular set,
the LHAPDF interface is designed to accommodate PDFs with uncertainties
as well as ``central fit'' PDFs.
For PDFs with uncertainties the PDF set represents one ``fit'' to the
data. For instance, a random sampling PDF set would use Eq.
(\ref{eqn:pdf;qcdsm;density}).
In other words for each PDF in the set the observable is calculated. 
The set of resulting predictions of the observable build up the
probability density.
The individual member PDFs of the set are needed to calculate
the PDF uncertainty on the observable. All PDF sets are defined through
external files. This means that a new set can be added by simply downloading
its file while the LHAPDF interface code does not change. The evolution
code is not part of LHAPDF. The current default choice included and 
interfaced is QCDNUM \cite{QCDNUM}. 
Each group that contributes PDF sets can provide 
their own evolution code; or they can employ QCDNUM, which is available in
the package.

\subsubsection{The philosophy}

The Les Houches Accord Parton Distribution Function interface
was conceived at the Les Houches 2001 workshop in the PDF working group to
enable the usage of Parton Distribution Functions with uncertainties
in a uniform manner. 
When PDFs with uncertainties are considered, a ``fit'' to
the data no longer is described by a single PDF. Instead in its
most flexible implementation, a fit is represented by a PDF set
consisting of many individual PDF members. Calculating the observable
for all the PDF members enables one to reconstruct the uncertainty
on the observable. The LHAPDF interface was made with this in mind 
and manipulates PDF
sets. 

The LHAPDF interface can be viewed as a successor to PDFLIB and improvements
were added. To list some of the features:
\begin{itemize}
\item
  The handling of PDF sets to enable PDF fits that include uncertainties.
\item
  The default evolution codes are benchmarked and compared for accuracy.
  Apart from accuracy another important feature of the evolution code
  is speed. Currently the default for evolution program is QCDNUM.
\item
  All PDF sets are defined through external files in parametrized form.
  This means the files are compact compared to storing the PDFs in a 
  grid format. Also
  new PDF sets can be defined by constructing the PDF defining files.
  The actual LHAPDF code does not have to be changed.
\item
  The LHAPDF code is modular and default choices like the QCDNUM
  evolution code can be easily altered.
\end{itemize}

Note that the current ``best fit'' PDFs can be viewed as PDF sets with
one member PDF and can be easily defined through the PDF set external
file definition. Alternatively one can group these ``fits'' is single
sets (e.g. MRST98 set) as they often represent best fits given different
model assumptions and as such reflect theoretical modelling uncertainties.

The first version of the code is available in Fortran from 
{\bf\tt http://pdf.fnal.gov}.

\subsubsection{Interfacing with LHAPDF}
The interface of LHAPDF with an external code is easy. We will
describe the basic steps sufficient for most applications. The 
web site contains more detailed information about additional
function calls. The function calls described here will be not
be altered in any way in future versions. 
Including the LHAPDF evolution code into a program involves 
three steps:
\begin{enumerate}
\item First one has to setup the LHAPDF interface code:
\begin{center}{\bf\tt call InitPDFset({\it name})}\end{center}
  It is called only once at the beginning of the code.
  The string variable {\it name} is the file name of the 
  external PDF file that defines the PDF set.
  For the standard evolution code QCDNUM it will 
  either calculate or read from file the LO/NLO splitting 
  function weights. The calculation of the weights might
  take some time depending on the chosen grid size. However,
  after the first run a grid file is created. Subsequent
  runs will use this file to read in the weights so that the lengthy 
  calculation of these weights is avoided. The file depends
  on the grid parameters and flavor thresholds. This means different
  PDF sets can have different grid files. The name of the grid file
  is specified in the PDF setup file. 
\item To use a individual PDF member it has to be initialized:
  \begin{center}{\bf\tt call InitPDF({\it mem)}}\end{center}
  The integer {\it mem} specifies the member PDF number. This
  routine only needs to be called when changing to a new PDF
  member. The time needed to setup depends on the evolution code used
  For QCDNUM the grid size is the determining factor. Note that 
  {\it mem}=0 selects the ``best fit'' PDF.
\item Once the PDF member is initialized one can call the evolution
  codes which will use the selected member. 

  The function call
  \begin{center}{\bf\tt function alphasPDF({\it Q})}\end{center}
  returns the value of $\alpha_S(Q)$ at double precision scale $Q$. 
  Note that its value can change between different PDF members.

  The subroutine call
  \begin{center}{\bf\tt call evolvePDF({\it x,Q,f})}\end{center}
  returns the PDF momentum densities $f$ (i.e. $x\times$ PDF number density)
  at double precision momentum fraction $x$ and double precision scale $Q$.
  The double precision array {\it f(-6:6)} will contain the momentum PDFs
  using the labelling convention of table \ref{tab:pdf;qcdsm;labeling}.
  \begin{table}[t]\hspace{2cm}\label{tab:pdf;qcdsm;labeling}\begin{center}
    \begin{tabular}{|l|r|r|r|r|r|r|r|r|r|r|r|r|r|}\hline
        & -6 & -5 & -4 & -3 & -2 & -1 & 0 &
      1 & 2 & 3 & 4 & 5 & 6 \\ \hline
      {\it P} & $\bar t$ & $\bar b$ & $\bar c$ & $\bar s$ 
      & $\bar u$ & $\bar d$ &
      $g$ & $d$ & $u$ & $s$ & $c$ & $b$ & $t$ \\
      {\it $\overline{\mbox{P}}$} & $t$ & $b$ & $c$ & $s$ & $u$ & $d$ & $g$ &
      $\bar d$ & $\bar u$ & $\bar s$ & $\bar c$ & $\bar b$ & $\bar t$ \\  
      \hline
    \end{tabular}
    \caption{The flavor enumeration convention used in the LHAPDF interface.
      (Note: CTEQ code use a different labelling scheme internally,
      with 1 $\leftrightarrow$ 2 and -1 $\leftrightarrow$ -2,
      but will adopt the above standard in the Les Houches interface.)}
  \end{center}\end{table}
  As long as the member PDF is not changed (by the 
  {\bf\tt call InitPDF} of step 2)
  the evolution calls will always use the same PDF member.
\end{enumerate}
A few additional calls can be useful:
\begin{itemize}
\item 
  To get the number of PDF members in the set: 
  \begin{center}{\bf\tt call numberPDF({\it Nmem})}.\end{center}
  The integer {\it Nmem} will contain the number of PDF members (excluding
  the special ``best fit'' member, i.e. the member numbers run
  from 0 to {\it Nmem}). 
\item
  Optionally the different PDF members can have weights which are obtained
  by:
  \begin{center}{\bf\tt call weightPDF({\it wgt})}.\end{center}
  The double precision variable {\it wgt} is for unweighted PDFs
  set to $1/Nmem$ such that the sum of all PDF member weights is unity.
  For weighted sets the use of the weights has to be defined by the
  method description.
\item 
  To get the evolution order of the PDFs:
  \begin{center}{\bf\tt call GetOrderPDF({\it order})}.\end{center}
  The integer variable {\it order} is 0 for Leading Order, 1 for
  Next-to-Leading Order, etc.
\item 
  To get the evolution order of $\alpha_S$:
  \begin{center}{\bf\tt call GetOrderAs({\it order})}.\end{center}
  The integer variable {\it order} is 0 for Leading Order, 1 for
  Next-to-Leading Order, etc.
\item 
  It is possible that during the PDF fitting the renormalization
  scale was chosen different from the factorization scale. The
  ratio of the renormalization scale over the factorization
  scale used in the ``fit'' can be obtained by
  \begin{center}{\bf\tt call GetRenFac(muf)}.\end{center}
  The double precision variable {\it muf} contains the ratio.
  Usually {\it muf} is equal to unity.
\item
  To get a description of the PDF set:
  \begin{center}{\bf\tt call GetDesc()}.\end{center}
  This call will print the PDF description to the standard output stream
\item
  The quark masses can be obtained by:
  \begin{center}{\bf\tt call GetQmass(nf,mass)}.\end{center}
  The mass {\it mass} is returned for quark flavor {\it nf}.
  The quark masses are used in the $\alpha_S$ evolution.
\item
  The flavor thresholds in the PDF evolution can be obtained by:
  \begin{center}{\bf\tt call GetThreshold(nf,Q)}.\end{center}
  The flavor threshold {\it Q} is returned for flavor {\it nf}.
  If {\it Q=-1d0} flavor is not in the evolution (e.g. the top quark
  is usually not included in the evolution).
  If {\it Q=0d0} flavor is parametrized at the parametrization scale.
  For positive non-zero values of {\it Q} the value is set to the
  flavor threshold at which the PDF starts to evolve.
\item
  The call returns the number of flavors used in the PDF:
  \begin{center}{\bf\tt call GetNf(nfmax)}.\end{center}
  Usually the returned value for {\it nfmax} 
  is equal to five as the top quark is usually not considered in the PDFs.
\end{itemize}

\subsubsection{An example}
A very simple example is given below. It accesses all member PDFs in the
set {\tt mypdf.LHpdf} and print out the $\alpha_S(M_Z)$ value and the
gluon PDF at several $(x,Q)$ points.

\begin{verbatim}
      program example
      implicit real*8(a-h,o-z)
      character*32 name
      real*8 f(-6:6)
*
      name='mypdf.LHpdf'
      call InitPDFset(name)
*
      QMZ=91.71d0
      write(*,*)
      call numberPDF(N)
      do i=1,N
         write(*,*) '---------------------------------------------'
         call InitPDF(i)
         write(*,*) 'PDF set ',i
         write(*,*)
         a=alphasPDF(QMZ)
         write(*,*) 'alphaS(MZ) = ',a
         write(*,*)
         write(*,*) 'x*Gluon'
         write(*,*) '   x     Q=10 GeV     Q=100 GeV    Q=1000 GeV'
         do x=0.01d0,0.095d0,0.01d0
            Q=10d0
            call evolvePDF(x,Q,f)
            g1=f(0)
            Q=100d0
            call evolvePDF(x,Q,f)
            g2=f(0)
            Q=1000d0
            call evolvePDF(x,Q,f)
            g3=f(0)
            write(*,*) x,g1,g2,g3
         enddo
      enddo
*
      end
\end{verbatim}

% ----------------------------------------------------------------------------

\subsection{Reference results for the evolution of parton distributions%
\protect\footnote{Contributing authors: G. Salam, A. Vogt}} 

In this section we provide a new set of benchmark tables for the evolution of 
unpolarized parton distributions of hadrons in perturbative QCD. Unlike the 
only comparable previous study \cite{Blumlein:1996rp}, we include results for 
unequal factorization and renormalization scales, $\mu_{\rm f} \!\neq\! 
\mu_{\rm r}$, and for the evolution with a variable number of partonic flavours 
$N_{\rm f}$. Besides the standard LO and NLO approximations, we also present 
the evolution including the (still approximate) NNLO splitting functions and 
the corresponding non-trivial second-order matching conditions at the 
heavy-quark thresholds. Our reference results are computed using two entirely 
independent and conceptually different evolution programs which, however, agree 
to better than 1 part in $10^5$ for momentum fractions $10^{-8} < x < 0.9$.

\subsubsection{Evolution equations and their solutions}
\label{gsav-s1}

At N$^{\rm m}$LO the scale dependence (`evolution') of the parton distributions 
$ f_p (x,\mu_{\rm f}^2) \equiv p(x,\mu_{\rm f}^2) $, where $ p = q_i, 
\,\bar{q}_i\, , g $ with $ i = 1, \,\ldots , N_{\rm f} $, is governed by the 
$ 2N_{\rm f}\! +\! 1$ coupled integro-differential equations
%
\begin{equation}
\label{gsav-eq1}
 \frac {d\, f_p(x,\mu_{\rm f}^2)} {d \ln \mu_{\rm f}^2} \: = \: 
 \sum_{l=0}^m \, a_{\rm s}^{l+1}(\mu_{\rm r}^2)\:
 {\cal P}_{pp'}^{(l)}\bigg(x,\frac{\mu_{\rm f}^2}{\mu_{\rm r}^2}\bigg) 
 \otimes f_{p'}(x,\mu_{\rm f}^2) \:\: .
\end{equation}
%
Here $\otimes$ denotes the Mellin convolution in the fractional-momentum 
variable $x$,
%
%\begin{equation}
%  a(x,\,\ldots ) \otimes b(x\,\ldots) \:\equiv\: \int_x^1 \! \frac{dy}{y} \: 
%  a(y,\,\ldots )\, b\bigg(\frac{x}{y},\,\ldots \bigg) \:\: ,
%\end{equation}
%
and summation over $p'$ is understood. % in Eq.~(\ref{gsav-eq1}). 
The scale dependence of the strong coupling $a_{\rm s} \equiv \alpha_{\rm s}/
(4\pi)$ is given by 
%
\begin{equation}
\label{gsav-eq2}
  \frac{d\, a_{\rm s}}{d \ln \mu_r^2} \: = \: \beta^{}_{\rm N^mLO}(a_{\rm s})
  \: = \: - \sum_{l=0}^m \, a_{\rm s}^{l+2} \,\beta_l \:\: .
\end{equation}
%
The general splitting functions ${\cal P}^{(l)}$ in Eq.~(\ref{gsav-eq1}) can be
reduced to the simpler expressions $P^{(l)}(x)$ at $\mu_{\rm r} = \mu_{\rm f}$. 
Up to NNLO ($\equiv$ N$^{\rm 2\,}$LO) the corresponding relations read
%
\begin{eqnarray}
\label{gsav-eq3}
  {\cal P}^{(0)} \bigg(x, \frac{\mu_{\rm f}^2}{\mu_{\rm r}^2} \bigg) 
  &\! =\! & P^{(0)}(x) \nonumber \\
  {\cal P}^{(1)} \bigg(x, \frac{\mu_{\rm f}^2}{\mu_{\rm r}^2} \bigg)
  &\! =\! & P^{(1)}(x) 
  - \beta_0 P^{(0)}(x)\, \ln \frac{\mu_{\rm f}^2}{\mu_{\rm r}^2} \\
  {\cal P}^{(2)} \bigg(x, \frac{\mu_{\rm f}^2}{\mu_{\rm r}^2} \bigg)
  &\! =\! & P^{(2)}(x) - \bigg\{ \beta_1 P^{(0)}(x) + 2\beta_0 P^{(1)}(x) 
    \bigg\} \ln \frac{\mu_{\rm f}^2}{\mu_{\rm r}^2} 
    + \beta_0^2 P^{(0)}(x)\, \ln^2 \frac{\mu_{\rm f}^2}{\mu_{\rm r}^2} \:\: .
  \nonumber
\end{eqnarray}
%
The generalization to higher orders is straightforward but not required for the
present calculations.

The LO and NLO coefficients $\beta_0$ and $\beta_1$ of the $\beta$-function in
Eq.~(\ref{gsav-eq2}) and the corresponding splitting functions $P^{(0)}(x)$ and 
$P^{(1)}(x)$ in Eq.~(\ref{gsav-eq3}) have been known for a long time, 
see Ref.~\cite{Furmanski:1982cw} and references therein. 
In the $\overline{\rm MS}$ scheme adopted here, the coefficient $\beta_2$ has 
been calculated in Refs.~\cite{Tarasov:1980au,Larin:1993tp}. The NNLO 
quantities $P^{(2)}(x)$ have not been computed so far, but approximate 
expressions have been constructed \cite{vanNeerven:1999ca,vanNeerven:2000uj,%
vanNeerven:2000wp} from the available partial results \cite{Larin:1994vu,%
Larin:1997wd,Retey:2000nq,Catani:1994sq,Blumlein:1996jp,Fadin:1998py,%
Ciafaloni:1998gs,Gracey:1994nn,Bennett:1998ch}.

An obvious and widespread approach to Eq.~(\ref{gsav-eq1}) is the direct 
numerical solution by discretization in both $x$ and $\mu_{\rm f}^2$. This 
method is also used by one of us (G.S.).
%  
The parton distributions are represented on a grid with $n$ points
uniformly spaced in $\ln 1/x$.  Their values at arbitrary $x$ are 
defined to be equal to a $p^\mathrm{th}$ order interpolation of 
neighbouring grid points.  The splitting functions can then be 
represented as sparse matrices acting on the vector of grid points. At
initialization the program takes a set of subroutines for the splitting
functions and calculates the corresponding matrices with a Gaussian 
adaptive integrator.  At each value of $\mu_{\rm f}^2$ the derivatives 
of the parton distributions are calculated through matrix 
multiplication and the evolution is carried out with a Runge-Kutta 
method. The algorithm has partial overlap with those of 
Refs.~\cite{QCDNUM,Santorelli:1998yt,Ratcliffe:2000kp,Pascaud:2001bi} 
and is described in more detail in Appendix~F of 
Ref.~\cite{Dasgupta:2001eq}.

For the reference tables presented below, the program has been run
with $7^{th}$ order interpolation in $x$ and multiple $x$-grids: one 
for $10^{-8}<x<1$ with $750$ points, another for $0.135<x<1$ with $240$
points and a third for $0.6<x<1$ with $180$ points. A grid uniform in
$\ln \mu_{\rm f}^2$ has been used with $220$ points in the range 
$ 2\,\mbox{ GeV}^2 < \mu_{\rm f}^2 < 10^6 \mbox{ GeV}^2 $. Halving the 
density of points in both $x$ and $\mu_{\rm f}^2$ leaves the results 
unchanged at the level of better than $1$ part in $10^{5}$ in the range 
$10^{-8}<x<0.9$ (except close to sign-changes of parton distributions).

An important alternative to this direct numerical treatment of 
Eq.~(\ref{gsav-eq1}) is the Mellin-$N$ moment solution in terms of a power 
expansion. This method is employed by the second author (A.V.). Here
Eq.~(\ref{gsav-eq1}) is transformed to $N$-space (reducing the convolution to a 
simple product) and $\mu_{\rm f}^2$ is replaced by $a_{\rm s}$ as the 
independent variable, assuming that $\mu_{\rm f}^{}/\mu_{\rm r}^{}$ is a fixed 
number. Expanding the resulting r.h.s.\ into a power series in $a_{\rm s}$, 
one arrives at
%
\begin{equation}
\label{gsav-eq4}
  \frac {d\, f_p(N,a_{\rm s})} {d\, a_{\rm s}} \: = \: - \sum_{l=0}^\infty 
  \, a_{\rm s}^{l-1}\, R_{pp'}^{(l)}(N) \, f_{p'}(N,a_{\rm s}) 
\end{equation}
%
with
%
\begin{equation}
\label{gsav-eq5}
  R_{pp'}^{(0)}(N)  \: = \: \frac{1}{\beta_0}\, P_{pp'}^{(0)}(N)
  \:\: , \quad\quad
  R_{pp'}^{\, (k\geq 1)} \: = \: \frac{1}{\beta_0}\, {\cal P}_{pp'}^{(k)}(N)
  - \sum_{l=1}^{k} \,\frac{\beta_l}{\beta_0}\, R_{pp'}^{(k-l)}(N) \:\: .
\end{equation}
%
At N$^{\rm m}$LO only the coefficients $\beta_{\, l\leq m}$ and 
${\cal P}^{(l\leq m)}$ are retained in Eq.~(\ref{gsav-eq5}). The solution of 
Eq.~(\ref{gsav-eq4}) can be expressed as an expansion around the LO result
%
\begin{equation}
\label{gsav-eq6}
  f(N,\mu_{\rm f}^2) \: = \: 
  \bigg[ 1 + \sum_{k=1}^\infty \, a^k_{\rm s}\, U_k(N) \bigg] 
  \left( \frac{a_{\rm s}}{a^{}_0} \right)^{-R_0(N)}
  \bigg[ 1 + \sum_{k=1}^\infty \, a^k_0 \, U_k(N) \bigg]^{-1} 
  f(N,\mu_{\rm f,0}^2) \:\: ,
\end{equation}
%
where $\mu_{\rm f,0}^2$ is the initial scale for the evolution, and $a^{}_0 
\equiv a_{\rm s}(\mu_{\rm r}^2(\mu_{\rm f,0}^2))$.
It is understood in Eq.~(\ref{gsav-eq6}) that the matrix structure is 
simplified by switching to the appropriate flavour singlet and non-singlet
combinations. For the explicit construction of the remaining $2\!\times\! 2$
matrices $U_k$ the reader is referred to Section 5 of 
Ref.~\cite{Blumlein:1998em}. Finally the solutions $f_p(N,\mu_{\rm f}^2)$ are 
transformed back to $x$-space by
%
\begin{equation}
\label{gsav-eq7}
  f_p (x,\mu_{\rm f}^2) \: = \: \frac{1}{\pi}\, \int_0^\infty \! dz \:
  {\rm Im} \, \Big[ \, e^{i\phi} x^{- c - z \exp(i\phi)} 
  f_p (N\! =\! c\! +\! z \exp(i\phi),\, \mu_{\rm f}^2) \, \Big] \:\: .
\end{equation}
%

The Mellin inversions (\ref{gsav-eq7}) can be performed with a sufficient 
accuracy using a fixed chain of Gauss quadratures. Hence the quantities 
$R_0(N)$ and  $U_k(N)$ in Eq.~(\ref{gsav-eq6}) have to the computed only once 
for the corresponding support points at the initialization of the program,
rendering the $N$-space evolution competitive in speed with fast $x$-space
codes.
Except where the parton distributions become very small, an accuracy of 1 part
in $10^5$ or better is achieved by including contributions up to $k = 15$ in 
Eq.~(\ref{gsav-eq6}) and using at most 18 eight-point Gauss quadratures with,
e.g., $z_{\rm max} = 70$, $c = 2$, $\phi = 3\pi/4$ in Eq.~(\ref{gsav-eq7}).

The two methods for solving Eq.~(\ref{gsav-eq1}) discussed above are completely
equivalent, i.e., they do not differ even by terms beyond N$^{\rm m}$LO, 
provided that the coupling $a_{\rm s}$ evolves exactly according to 
Eq.~(\ref{gsav-eq2}). This condition is fulfilled in the present calculations.
Thus the results of our two programs can be compared directly, yielding a
powerful check of the correctness and accuracy of the codes. Note that the 
only previously published high-precision comparison of NLO evolution programs 
\cite{Blumlein:1996rp} used the truncated expansion of $a_{\rm s}^2
(\mu_{\rm r})$ in terms of inverse powers of $\ln \mu_{\rm r}^2 / \Lambda^2$ 
which does not exactly conform to Eq.~(\ref{gsav-eq2}). Consequently such 
direct comparisons were not possible in Ref.~\cite{Blumlein:1996rp}.

Following Ref.~\cite{Furmanski:1982cw}, the $N$-space solution (\ref{gsav-eq6})
has usually been subjected to a further expansion in the coupling constants, 
retaining only the terms up to order $m$ in the product of the $U$-matrices. 
Only the terms thus kept in Eq.~(\ref{gsav-eq6}) are free from contributions by 
the higher-order coefficients $\beta_{k>m}$ and $P^{(k>m)}$, and only these 
terms combine to factorization-scheme independent quantities when the parton 
distributions are convoluted with the N$^{\rm m}$LO partonic cross sections.
Reference results for the truncated solution will be presented elsewhere 
\cite{GSAV}.


\subsubsection{Initial conditions and heavy-quark treatment}
\label{gsav-s2}

The following initial conditions for the reference results have been set up
at the Les Houches meeting: The evolution is started at 
%
\begin{equation}
\label{gsav-eq8}
  \mu_{\rm f,0}^2 \: = \: 2 \mbox{ GeV}^2 \:\: .
\end{equation}
%
Roughly along the lines of the CTEQ5M parametrization \cite{Lai:1999wy}, the 
input distributions are chosen as
%
\begin{eqnarray}
\label{gsav-eq9}
  xu_v(x,\mu_{\rm f,0}^2)       &\! =\! & 5.107200\: x^{0.8}\: (1-x)^3  
    \nonumber \\
  xd_v(x,\mu_{\rm f,0}^2)       &\! =\! & 3.064320\: x^{0.8}\: (1-x)^4  
    \nonumber \\
  xg\,(x,\mu_{\rm f,0}^2)       &\! =\! & 1.700000\, x^{-0.1} (1-x)^5 
    \\
  x\bar{d}\,(x,\mu_{\rm f,0}^2) &\! =\! & .1939875\, x^{-0.1} (1-x)^6
    \nonumber\\
  x\bar{u}\,(x,\mu_{\rm f,0}^2) &\! =\! & (1-x)\: x\bar{d}\,(x,\mu_{\rm f,0}^2)
    \nonumber\\
  xs\,(x,\mu_{\rm f,0}^2)       &\! =\! & x\bar{s}\,(x,\mu_{\rm f,0}^2) 
    \: = \: 0.2\, x(\bar{u}+\bar{d}\,)(x,\mu_{\rm f,0}^2) 
    \nonumber  
\end{eqnarray}
%
where, as usual, $q_{i,v} \equiv q_i - \bar{q}_i$. The running couplings are 
specified via
% 
\begin{equation}
\label{gsav-eq10}
  \alpha_{\rm s}(\mu_{\rm r}^2\! =\! 2\mbox{ GeV}^2) \: = \: 0.35 \:\: .
\end{equation}
%
For simplicity these initial conditions are employed regardless of the order of 
the evolution and the ratio of the renormalization and factorization scales. 
At LO this ratio is fixed to unity, beyond LO we use
%
\begin{equation}
\label{gsav-eq11}
  \mu_{\rm r}^2 \: = \: k_{\rm r} \mu_{\rm f}^2 \:\: , \quad\quad 
  k_{\rm r} \: = \:  0.5\, , \:\: 1\, , \:\: 2 \:\: .
\end{equation}
 
For the evolution with a fixed number $N_{\rm f} > 3$ of quark flavours the 
quark distributions not specified in Eq.~(\ref{gsav-eq9}) are assumed to vanish 
at $\mu_{\rm f,0}^2$, and Eq.~(\ref{gsav-eq10}) is understood to refer to the 
chosen value of~$N_{\rm f}$. For the evolution with a variable $ N_{\rm f} = 3 
\ldots 6 $, Eqs.~(\ref{gsav-eq9}) and (\ref{gsav-eq10}) always refer to three 
flavours. $N_{\rm f}$ is then increased by one unit at the heavy-quark pole 
masses taken as
%
\begin{equation}
\label{gsav-eq12}
  m_{\rm c} \: = \: \mu_{\rm f,0}      \: , \quad 
  m_{\rm b} \: = \: 4.5 \mbox{ GeV}^2  \: , \quad
  m_{\rm t} \: = \: 175 \mbox{ GeV}^2  \:\: ,
\end{equation}
%
i.e., the evolution is performed as discussed in Section \ref{gsav-s1} between 
these thresholds, and the respective matching conditions are invoked at 
$\mu_{\rm f}^2 = m_{\rm h}^2$, $h = c,\, b,\, t$. For the parton distributions 
these conditions have been derived in Ref.~\cite{Buza:1998wv}. Up to 
N$^{\rm m=2\,}$LO they read 
%
\begin{equation}
\label{gsav-eq13}
  l_i^{\,(N_{\rm f}+1)}(x,m_h^2) \: = \:  l_i^{\,(N_{\rm f})}(x,m_h^2) + 
  \delta_{m2} \: a_{\rm s}^2\: A^{\rm NS,(2)}_{qq,h}(x) \otimes 
  l_i^{\, (N_{\rm f})}(x,m_h^2) 
\end{equation}
%
where $l = q,\, \bar{q}$ and $i = 1,\ldots N_{\rm f}$, and
\begin{eqnarray}
\label{gsav-eq14}
  g^{\, (N_{\rm f}+1)}(x,m_h^2) \:\:\: &\! = \!\! & 
    g^{\, (N_{\rm f})}(x,m_h^2) + \nonumber \\[0.5mm] & & 
    \delta_{m2} \, a_{\rm s}^2\, \Big[  
    A_{\rm gq,h}^{S,(2)}(x) \otimes \Sigma^{\, (N_{\rm f})}(x,m_h^2) +
    A_{\rm gg,h}^{S,(2)}(x) \otimes g^{\, (N_{\rm f})}(x,m_h^2) \Big ] 
  \\[0.5mm]
  (h+\bar{h})^{\, (N_{\rm f}+1)}(x,m_h^2) \! &\! =\!\! & 
    \delta_{m2} \, a_{\rm s}^2\: \Big [ 
    \tilde{A}_{\rm hq}^{S,(2)}(x) \otimes \Sigma^{\, (N_{\rm f})}(x,m_h^2) +
    \tilde{A}_{\rm hg}^{S,(2)}(x) \otimes g^{\, (N_{\rm f})}(x,m_h^2) \Big  ] 
  \quad \nonumber
\end{eqnarray}
%
with $\Sigma^{\,(N_{\rm f})} \equiv \sum_{i=1}^{N_{\rm f}} (q_i + \bar{q}_i)$
and $h = \bar{h}$.  The coefficients $A^{(2)}$ can be found in Appendix B of 
ref.~\cite{Buza:1998wv} -- due to our choice of $\mu_{\rm f}^2 = m_{\rm h}^2$
for the thresholds only the scale-independent parts of the expressions are 
needed here -- from where the notation for these coefficients has been taken 
over. 
%
The corresponding N$^{\rm m\,}$LO relation for the coupling constant 
\cite{Larin:1995va,Chetyrkin:1997sg} is given by
%
\begin{equation}
\label{gsav-eq15}
  a_{\rm s}^{\, (N_{\rm f}+1)}(k_{\rm r} m_h^2) \: = \: 
  a_{\rm s}^{\, (N_{\rm f})}(k_{\rm r} m_h^2) + \sum_{n=1}^{m} \, 
  \Big( a_{\rm s}^{\, (N_{\rm f})}(k_{\rm r} m_h^2) \Big)^{n+1} 
  \sum_{l=0}^{n} \, c^{}_{n,l} \, \ln k_{\rm r} \:\: .
\end{equation}
%
The pole-mass coefficients $c_{n,l}$ in Eq.~(\ref{gsav-eq15}) can be inferred 
from Eq.~(9) of Ref.~\cite{Chetyrkin:1997sg}, where $4\, a_{\rm s}^{\, 
(N_{\rm f}-1)}$ is expressed in terms of $4\,a_{\rm s}^{\, (N_{\rm f})}$.
Note that we use $a_{\rm s}^{\, (N_{\rm f}+1)}(k_{\rm r} m_h^2)$ on the r.h.s.\
of Eq.~(\ref{gsav-eq14}).

\subsubsection{The benchmark results}
\label{gsav-s3}

We have compared the results of our two evolution programs, under the
conditions specified in Section \ref{gsav-s2}, at 500 $x$-$\mu_{\rm f}^2$
points covering the range $10^{-8} \leq x \leq 0.9$ and $2 \mbox{ GeV}^2 \leq
\mu_{\rm f}^2 \leq 10^6 \mbox{ GeV}^2 $. A representative subset of our results 
at $\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^4$, a scale relevant to high-$E_T$ jets
at {\sc Tevatron} and close to $m_{\rm W}^2$, $m_{\rm Z}^2$ and, possibly,
$m_{\rm Higgs}^2$, is presented in Tables \ref{gsav-t1}--\ref{gsav-t5}.
These results are given in terms of the valence distributions, defined below 
Eq.~(\ref{gsav-eq9}), $L_{\pm} \equiv \bar{d} \pm \bar{u}$, and the quark-%
antiquark sums $q_+ \! \equiv q \! -\!\bar{q}$ for $q = s,\, c$ and, for the 
variable-$N_{\rm f}$ case, $b$.

For compactness an abbreviated notation is employed throughout the tables, 
i.e., all numbers $a\cdot 10^b$ are written as $a^b$. In the vast majority of 
the $x$-$\mu_{\rm f}^2$ points our results are found to agree to all five 
figures displayed, except for the tiny NLO and NNLO sea-quark distributions at 
$x = 0.9$, in the tables. In fact, the differences for $x < 0.9$ are not larger 
than $\pm 1$ in the sixth digit, i.e, the offsets are smaller than 1 part in 
$10^5$. Entries where these residual offsets lead to a different fifth digit 
after rounding are indicated by the subscript '$*$'. The number with the 
smaller modulus is then given in the tables, e.g., $1.1111_*^1$ should be read 
as $1.11115 \cdot 10^1$ with an uncertainty of $\pm 1$ in the last figure. 

As mentioned in Section \ref{gsav-s1}, the three-loop (NNLO) splitting 
functions $P^{(2)}(x)$ in Eq.~(\ref{gsav-eq3}) are not yet exactly known. For
the NNLO reference results in Tables \ref{gsav-t4} and \ref{gsav-t5} we have
resorted to the average of the two extreme approximations constructed in 
Ref.~\cite{vanNeerven:2000wp}. The remaining uncertainties of these results 
can be estimated by instead employing the extreme approximations itselves.
Their relative effects are illustrated, also at $\mu_{\rm f}^2 = 10^4 
\mbox{ GeV}^4$, in Fig.~\ref{gsav-f1}. The uncertainties of the NNLO evolution
to this scale turn out to be virtually negligible at $x \geq 0.05$ and to 
amount to less than $\pm 1\%$ down to $x \simeq 10^{-4}$.

\begin{figure}[hbt]
\vspace{2mm}
\begin{center}
\mbox{\includegraphics[height=12.0cm]{pdffig.eps}}
\end{center}
\vspace{-7mm}
\caption{The relative effects of the present uncertainties of the NNLO 
 splitting functions on the evolution of the input (\ref{gsav-eq8}) -- 
 (\ref{gsav-eq10}) for $\mu_{\rm r}=\mu_{\rm f}$, estimated by using the 
 extreme approximations `A' and `B' of Ref.~\cite{vanNeerven:2000wp} instead of 
 their average.}
\label{gsav-f1}
\end{figure}  

\begin{table}[htp]
\caption{Reference results for the $N_{\rm f}=4$ (FFN) and the 
 variable-$N_{\rm f}$ (VFN) leading-order evolution for the initial conditions
 (\ref{gsav-eq8}) -- (\ref{gsav-eq10}), shown together with the input parton
 distributions (\ref{gsav-eq9}). The respective values for $\alpha_{\rm s}
 (\mu_{\rm r}^2\! =\! \mu_{\rm f}^2\! =\! 10^4 \mbox{ GeV}^2)$ read $0.117374$ 
 (FFN) and $0.122306$ (VFN). For the notation see the first two paragraphs of 
 Section \ref{gsav-s3}.} 
\label{gsav-t1}
%
\begin{center}
\vspace{5mm}
\begin{tabular}{||c||r|r|r|r|r|r|r|r||}
\hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{1}{||c||}{$x$} &
\multicolumn{1}{c|} {$xu_v$} &
\multicolumn{1}{c|} {$xd_v$} &
\multicolumn{1}{c|} {$xL_-$} &
\multicolumn{1}{c|} {$xL_+$} &
\multicolumn{1}{c|} {$xs_+$} &
\multicolumn{1}{c|} {$xc_+$} &
\multicolumn{1}{c|} {$xb_+$} &
\multicolumn{1}{c||}{$xg$} \\[0.5mm] \hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{Input, $\,\mu_{\rm f}^2 = 2 \mbox{ GeV}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.2829^{-5}\!$ &$\! 7.6972^{-6}\!$ &$\! 9.7224^{-8}\!$ &$\! 3.8890^{+0}\!$ &
$\! 7.7779^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 8.5202^{+0}\!$\\
$\! 10^{-6}\!$ &
$\! 8.0943^{-5}\!$ &$\! 4.8566^{-5}\!$ &$\! 7.7227^{-7}\!$ &$\! 3.0891^{+0}\!$ &
$\! 6.1782^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 6.7678^{+0}\!$\\
$\! 10^{-5}\!$ &
$\! 5.1070^{-4}\!$ &$\! 3.0642^{-4}\!$ &$\! 6.1341^{-6}\!$ &$\! 2.4536^{+0}\!$ &
$\! 4.9072^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 5.3756^{+0}\!$\\
$\! 10^{-4}\!$ &
$\! 3.2215^{-3}\!$ &$\! 1.9327^{-3}\!$ &$\! 4.8698^{-5}\!$ &$\! 1.9478^{+0}\!$ &
$\! 3.8957^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 4.2681^{+0}\!$\\
$\! 10^{-3}\!$ &
$\! 2.0271^{-2}\!$ &$\! 1.2151^{-2}\!$ &$\! 3.8474^{-4}\!$ &$\! 1.5382^{+0}\!$ &
$\! 3.0764^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 3.3750^{+0}\!$\\
$\! 10^{-2}\!$ &
$\! 1.2448^{-1}\!$ &$\! 7.3939^{-2}\!$ &$\! 2.8946^{-3}\!$ &$\! 1.1520^{+0}\!$ &
$\! 2.3041^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 2.5623^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.9008^{-1}\!$ &$\! 3.1864^{-1}\!$ &$\! 1.2979^{-2}\!$ &$\! 4.9319^{-1}\!$ &
$\! 9.8638^{-2}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 1.2638^{+0}\!$\\
$\! 0.3    \!$ &
$\! 6.6861^{-1}\!$ &$\! 2.8082^{-1}\!$ &$\! 7.7227^{-3}\!$ &$\! 8.7524^{-2}\!$ &
$\! 1.7505^{-2}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 3.2228^{-1}\!$\\
$\! 0.5    \!$ &
$\! 3.6666^{-1}\!$ &$\! 1.1000^{-1}\!$ &$\! 1.6243^{-3}\!$ &$\! 9.7458^{-3}\!$ &
$\! 1.9492^{-3}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 5.6938^{-2}\!$\\
$\! 0.7    \!$ &
$\! 1.0366^{-1}\!$ &$\! 1.8659^{-2}\!$ &$\! 1.0259^{-4}\!$ &$\! 3.8103^{-4}\!$ &
$\! 7.6207^{-5}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 4.2810^{-3}\!$\\
$\! 0.9    \!$ &
$\! 4.6944^{-3}\!$ &$\! 2.8166^{-4}\!$ &$\! 1.7644^{-7}\!$ &$\! 4.3129^{-7}\!$ &
$\! 8.6259^{-8}\!$ &$\! 0.0   ^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 1.7180^{-5}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{LO, $\, N_{\rm f} = 4\,$, 
                       $\,\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 5.7722^{-5}\!$ &$\! 3.4343^{-5}\!$ &$\! 7.6527^{-7}\!$ &$\! 9.9465^{+1}\!$ &
$\! 4.8642^{+1}\!$ &$\! 4.7914^{+1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 1.3162^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 3.3373^{-4}\!$ &$\! 1.9800^{-4}\!$ &$\! 5.0137^{-6}\!$ &$\! 5.0259^{+1}\!$ &
$\! 2.4263^{+1}\!$ &$\! 2.3685^{+1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 6.0008^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 1.8724^{-3}\!$ &$\! 1.1065^{-3}\!$ &$\! 3.1696^{-5}\!$ &$\! 2.4378^{+1}\!$ &
$\! 1.1501^{+1}\!$ &$\! 1.1042^{+1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 2.5419^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.0057^{-2}\!$ &$\! 5.9076^{-3}\!$ &$\! 1.9071^{-4}\!$ &$\! 1.1323^{+1}\!$ &
$\! 5.1164^{+0}\!$ &$\! 4.7530^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 9.7371^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.0392^{-2}\!$ &$\! 2.9296^{-2}\!$ &$\! 1.0618^{-3}\!$ &$\! 5.0324^{+0}\!$ &
$\! 2.0918^{+0}\!$ &$\! 1.8089^{+0}\!$ &$\! 0.0   ^{+0}\!$ &$\! 3.2078^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.1955^{-1}\!$ &$\! 1.2433^{-1}\!$ &$\! 4.9731^{-3}\!$ &$\! 2.0433^{+0}\!$ &
$\! 7.2814^{-1}\!$ &$\! 5.3247^{-1}\!$ &$\! 0.0   ^{+0}\!$ &$\! 8.0546^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.7267^{-1}\!$ &$\! 2.8413^{-1}\!$ &$\! 1.0470^{-2}\!$ &$\! 4.0832^{-1}\!$ &
$\! 1.1698^{-1}\!$ &$\! 5.8864^{-2}\!$ &$\! 0.0   ^{+0}\!$ &$\! 8.8766^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.7925^{-1}\!$ &$\! 1.4186^{-1}\!$ &$\! 3.3029^{-3}\!$ &$\! 4.0165^{-2}\!$ &
$\! 1.0516^{-2}\!$ &$\! 4.1379_*^{-3}\!$&$\!0.0   ^{+0}\!$ &$\! 8.2676^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.3476^{-1}\!$ &$\! 3.5364^{-2}\!$ &$\! 4.2815^{-4}\!$ &$\! 2.8624^{-3}\!$ &
$\! 7.3138^{-4}\!$ &$\! 2.6481^{-4}\!$ &$\! 0.0   ^{+0}\!$ &$\! 7.9240^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.3123^{-2}\!$ &$\! 3.5943^{-3}\!$ &$\! 1.5868^{-5}\!$ &$\! 6.8961^{-5}\!$ &
$\! 1.7725^{-5}\!$ &$\! 6.5549^{-6}\!$ &$\! 0.0   ^{+0}\!$ &$\! 3.7311^{-4}\!$\\
$\! 0.9    \!$ &
$\! 4.3443^{-4}\!$ &$\! 2.2287^{-5}\!$ &$\! 1.1042^{-8}\!$ &$\! 3.6293^{-8}\!$ &
$\! 1.0192^{-8}\!$ &$\! 4.8893_*^{-9}\!$&$\!0.0   ^{+0}\!$ &$\! 1.0918^{-6}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{LO, $\, N_{\rm f} = 3\ldots 5\,$, 
                        $\,\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 5.8771^{-5}\!$ &$\! 3.4963^{-5}\!$ &$\! 7.8233^{-7}\!$ &$\! 1.0181^{+2}\!$ &
$\! 4.9815^{+1}\!$ &$\! 4.9088^{+1}\!$ &$\! 4.6070^{+1}\!$ &$\! 1.3272^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 3.3933^{-4}\!$ &$\! 2.0129^{-4}\!$ &$\! 5.1142^{-6}\!$ &$\! 5.1182^{+1}\!$ &
$\! 2.4725^{+1}\!$ &$\! 2.4148^{+1}\!$ &$\! 2.2239^{+1}\!$ &$\! 6.0117^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 1.9006^{-3}\!$ &$\! 1.1229^{-3}\!$ &$\! 3.2249^{-5}\!$ &$\! 2.4693^{+1}\!$ &
$\! 1.1659^{+1}\!$ &$\! 1.1201^{+1}\!$ &$\! 1.0037^{+1}\!$ &$\! 2.5282^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.0186^{-2}\!$ &$\! 5.9819^{-3}\!$ &$\! 1.9345^{-4}\!$ &$\! 1.1406^{+1}\!$ &
$\! 5.1583^{+0}\!$ &$\! 4.7953^{+0}\!$ &$\! 4.1222^{+0}\!$ &$\! 9.6048^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.0893^{-2}\!$ &$\! 2.9576^{-2}\!$ &$\! 1.0730^{-3}\!$ &$\! 5.0424^{+0}\!$ &
$\! 2.0973^{+0}\!$ &$\! 1.8147^{+0}\!$ &$\! 1.4582^{+0}\!$ &$\! 3.1333^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.2080^{-1}\!$ &$\! 1.2497^{-1}\!$ &$\! 4.9985^{-3}\!$ &$\! 2.0381^{+0}\!$ &
$\! 7.2625^{-1}\!$ &$\! 5.3107^{-1}\!$ &$\! 3.8106^{-1}\!$ &$\! 7.7728^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.7166^{-1}\!$ &$\! 2.8334^{-1}\!$ &$\! 1.0428^{-2}\!$ &$\! 4.0496^{-1}\!$ &
$\! 1.1596^{-1}\!$ &$\! 5.8288^{-2}\!$ &$\! 3.5056^{-2}\!$ &$\! 8.4358^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.7597^{-1}\!$ &$\! 1.4044^{-1}\!$ &$\! 3.2629^{-3}\!$ &$\! 3.9592^{-2}\!$ &
$\! 1.0363^{-2}\!$ &$\! 4.0740^{-3}\!$ &$\! 2.2039^{-3}\!$ &$\! 7.8026^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.3284^{-1}\!$ &$\! 3.4802^{-2}\!$ &$\! 4.2031^{-4}\!$ &$\! 2.8066^{-3}\!$ &
$\! 7.1707^{-4}\!$ &$\! 2.5958^{-4}\!$ &$\! 1.3522^{-4}\!$ &$\! 7.4719^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.2643^{-2}\!$ &$\! 3.5134^{-3}\!$ &$\! 1.5468^{-5}\!$ &$\! 6.7201^{-5}\!$ &
$\! 1.7278^{-5}\!$ &$\! 6.3958^{-6}\!$ &$\! 3.3996^{-6}\!$ &$\! 3.5241^{-4}\!$\\
$\! 0.9    \!$ &
$\! 4.2047^{-4}\!$ &$\! 2.1529^{-5}\!$ &$\! 1.0635^{-8}\!$ &$\! 3.4998^{-8}\!$ &
$\! 9.8394^{-9}\!$ &$\! 4.7330^{-9}\!$ &$\! 2.8903^{-9}\!$ &$\! 1.0307^{-6}\!$\\
%
\hline \hline
%
\end{tabular}
\end{center}
\end{table}

\begin{table}[htp]
\caption{Reference results for the $N_{\rm f}=4$ next-to-leading-order 
 evolution for the initial conditions (\ref{gsav-eq8}) -- (\ref{gsav-eq10}). 
 The corresponding value of the strong coupling is $\alpha_{\rm s}
 (\mu_{\rm r}^2 \! =\! 10^4 \mbox{ GeV}^2) = 0.110902$. As in the leading-order 
 case, the valence distributions $s_v$ and $c_v$ vanish for the input 
 (\ref{gsav-eq9}). The notation is explained in the first two paragraphs of 
 Section \ref{gsav-s3}.}
\label{gsav-t2}
%
\vspace{4mm}
\begin{center}
\begin{tabular}{||c||r|r|r|r|r|r|r||}
\hline \hline
\multicolumn{8}{||c||}{} \\[-3mm]
\multicolumn{8}{||c||}{NLO, $\, N_{\rm f} = 4$, 
                       $\, \mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
\multicolumn{8}{||c||}{} \\[-0.35cm] \hline \hline
 & & & & & & & \\[-0.3cm]
\multicolumn{1}{||c||}{$x$} &
\multicolumn{1}{c|} {$xu_v$} &
\multicolumn{1}{c|} {$xd_v$} &
\multicolumn{1}{c|} {$xL_-$} &
\multicolumn{1}{c|} {$xL_+$} &
\multicolumn{1}{c|} {$xs_+$} &
\multicolumn{1}{c|} {$xc_+$} &
\multicolumn{1}{c||}{$xg$} \\ \hline \hline
\multicolumn{8}{||c||}{} \\[-3mm]
\multicolumn{8}{||c||}{$\mu_{\rm r}^2 = \mu_{\rm f}^2$} \\
\multicolumn{8}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
 $10^{-7}$ & $1.0616^{-4}$ & $6.2328^{-5}$ & $4.2440^{-6}$ & $1.3598^{+2}$
           & $6.6913_*^{+1}$&$6.6195^{+1}$ & $1.1483_*^{+3}$ \\
 $10^{-6}$ & $5.4177^{-4}$ & $3.1719^{-4}$ & $1.9241^{-5}$ & $6.8396^{+1}$
           & $3.3342^{+1}$ & $3.2771^{+1}$ & $5.3911^{+2}$ \\
 $10^{-5}$ & $2.6870^{-3}$ & $1.5677^{-3}$ & $8.3575^{-5}$ & $3.2728^{+1}$
           & $1.5685^{+1}$ & $1.5231^{+1}$ & $2.3528^{+2}$ \\
 $10^{-4}$ & $1.2841^{-2}$ & $7.4558^{-3}$ & $3.4911^{-4}$ & $1.4746^{+1}$ 
           & $6.8355^{+0}$ & $6.4769^{+0}$ & $9.2872_*^{+1}$ \\
 $10^{-3}$ & $5.7926^{-2}$ & $3.3337^{-2}$ & $1.4162^{-3}$ &$6.1648_*^{+0}$ 
           & $2.6659^{+0}$ & $2.3878^{+0}$ & $3.1502^{+1}$ \\
 $10^{-2}$ & $2.3026^{-1}$ & $1.2928^{-1}$ & $5.3251^{-3}$ & $2.2527^{+0}$  
           & $8.4220_*^{-1}$&$6.5246^{-1}$ & $8.1066^{+0}$ \\
 $0.1$     & $5.5452^{-1}$ & $2.7336^{-1}$ & $1.0011^{-2}$ & $3.9336^{-1}$
           & $1.1489^{-1}$ & $6.0351^{-2}$ & $8.9867^{-1}$ \\
 $0.3$     & $3.5393^{-1}$ & $1.3158^{-1}$ & $3.0362^{-3}$ & $3.5848^{-2}$
           & $9.2030^{-3}$ & $3.3890^{-3}$ & $8.3451^{-2}$ \\
 $0.5$     & $1.2271^{-1}$ & $3.1967^{-2}$ & $3.8265^{-4}$ & $2.4126^{-3}$
           & $5.8424^{-4}$ & $1.6955^{-4}$ & $8.0473^{-3}$ \\
 $0.7$     & $2.0429^{-2}$ & $3.1473_*^{-3}$&$1.3701^{-5}$ & $5.3622^{-5}$
           & $1.2393^{-5}$ & $2.7807^{-6}$ & $3.8721^{-4}$ \\
 $0.9$     & $3.6096^{-4}$ & $1.8317^{-5}$ & $8.923 ^{-9}$ & $2.092 ^{-8}$
           & $4.039 ^{-9}$ &$-2.405^{-10}$ & $1.2127^{-6}$ \\
\hline \hline
%
\multicolumn{8}{||c||}{} \\[-3mm]
\multicolumn{8}{||c||}{$\mu_{\rm r}^2 = 2\, \mu_{\rm f}^2$} \\
\multicolumn{8}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
 $10^{-7}$ & $9.2960^{-5}$ & $5.4699^{-5}$ & $3.3861^{-6}$ & $1.2214^{+2}$
           & $5.9987^{+1}$ & $5.9265^{+1}$ & $1.0911^{+3}$ \\
 $10^{-6}$ & $4.8463^{-4}$ & $2.8440^{-4}$ & $1.5820^{-5}$ & $6.1831^{+1}$
           & $3.0056^{+1}$ & $2.9483^{+1}$ & $5.1456^{+2}$ \\
 $10^{-5}$ & $2.4578^{-3}$ & $1.4374^{-3}$ & $7.1265^{-5}$ & $2.9845^{+1}$
           & $1.4240^{+1}$ & $1.3785^{+1}$ & $2.2580^{+2}$ \\
 $10^{-4}$ & $1.2018^{-2}$ & $6.9946^{-3}$ & $3.1111^{-4}$ & $1.3618^{+1}$
           & $6.2690^{+0}$ & $5.9088^{+0}$ & $8.9753^{+1}$ \\
 $10^{-3}$ & $5.5483^{-2}$ & $3.2009^{-2}$ & $1.3254^{-3}$ & $5.8076^{+0}$
           & $2.4848^{+0}$ & $2.2050^{+0}$ & $3.0729^{+1}$ \\
 $10^{-2}$ & $2.2595^{-1}$ & $1.2720^{-1}$ & $5.2141^{-3}$ & $2.1896^{+0}$
           & $8.0746^{-1}$ & $6.1564^{-1}$ & $8.0188^{+0}$ \\
 $0.1$     & $5.6007^{-1}$ & $2.7697^{-1}$ & $1.0180^{-2}$ & $3.9945^{-1}$
           & $1.1570^{-1}$ & $5.9661^{-2}$ & $9.1201^{-1}$ \\
 $0.3$     & $3.6474^{-1}$ & $1.3612^{-1}$ & $3.1588^{-3}$ & $3.7501^{-2}$
           & $9.6302^{-3}$ & $3.5499^{-3}$ & $8.6368^{-2}$ \\
 $0.5$     & $1.2843^{-1}$ & $3.3610^{-2}$ & $4.0510^{-4}$ & $2.5822^{-3}$
           & $6.3044^{-4}$ & $1.8999^{-4}$ & $8.4178^{-3}$ \\
 $0.7$     & $2.1779^{-2}$ & $3.3725^{-3}$ & $1.4798^{-5}$ & $5.9125^{-5}$
           & $1.3961^{-5}$ & $3.5593^{-6}$ & $4.0836^{-4}$ \\
 $0.9$     & $3.9817^{-4}$ & $2.0321^{-5}$ & $9.987 ^{-9}$ & $2.555 ^{-8}$
           & $5.586 ^{-9}$ & $7.930^{-10}$ & 1.3008$^{-6}$ \\
\hline \hline
%
\multicolumn{8}{||c||}{} \\[-3mm]
\multicolumn{8}{||c||}{$\mu_{\rm r}^2 = 1/2\, \mu_{\rm f}^2$} \\
\multicolumn{8}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
 $10^{-7}$ & $1.2438^{-4}$ &$ 7.2817^{-5}$ & $5.5568^{-6}$ & $1.4556^{+2}$
           & $7.1706^{+1}$ &$ 7.0990^{+1}$ & $1.1468^{+3}$ \\
 $10^{-6}$ & $6.1759^{-4}$ &$ 3.6051^{-4}$ & $2.4322^{-5}$ & $7.3406^{+1}$
           & $3.5851^{+1}$ &$ 3.5282_*^{+1}$&$5.4041^{+2}$ \\
 $10^{-5}$ & $2.9770^{-3}$ &$ 1.7316^{-3}$ & $1.0121_*^{-5}$&$3.5158_*^{+1}$
           & $1.6903^{+1}$ &$ 1.6452^{+1}$ & $2.3663^{+2}$ \\
 $10^{-4}$ & $1.3820^{-2}$ &$ 7.9998^{-3}$ & $4.0093^{-4}$ & $1.5795^{+1}$
           & $7.3626^{+0}$ &$ 7.0057^{+0}$ & $9.3640^{+1}$ \\
 $10^{-3}$ & $6.0585^{-2}$ &$ 3.4766^{-2}$ & $1.5300^{-3}$ & $6.5284^{+0}$
           & $2.8504^{+0}$ &$ 2.5740^{+0}$ & $3.1795^{+1}$ \\
 $10^{-2}$ & $2.3422_*^{-1}$&$1.3114^{-1}$ & $5.4411^{-3}$ & $2.3221^{+0}$
           & $8.8022^{-1}$ &$ 6.9260^{-1}$ & $8.1613^{+0}$ \\
 0.1       & $5.4824^{-1}$ &$ 2.6954^{-1}$ & $9.8435^{-2}$ & $3.8787^{-1}$
           & $1.1419^{-1}$ &$ 6.0997^{-2}$ & $8.9361^{-1}$ \\
 0.3       & $3.4425^{-1}$ &$ 1.2760^{-1}$ & $2.9317^{-3}$ & $3.4294^{-2}$
           & $8.7599^{-3}$ &$ 3.1681^{-3}$ & $8.2031^{-2}$ \\
 0.5       & $1.1794^{-1}$ &$ 3.0618^{-2}$ & $3.6454^{-4}$ & $2.2530^{-3}$
           & $5.3541^{-4}$ &$ 1.4134^{-4}$ & $7.8595^{-3}$ \\
 0.7       & $1.9356^{-2}$ &$ 2.9698^{-3}$ & $1.2847^{-5}$ & $4.8328^{-5}$
           & $1.0666^{-5}$ &$ 1.6668^{-6}$ & $3.7624_*^{-4}$ \\
 0.9       & $3.3264^{-4}$ &$ 1.6800^{-5}$ & $8.124 ^{-9}$ & $1.573 ^{-8}$
           & $2.024 ^{-9}$ &$-1.870 ^{-9}$ & $1.1647^{-6}$ \\
\hline \hline
%
\end{tabular}
\end{center}
\end{table}

\begin{table}[htp]
\caption{As Table \ref{gsav-t2}, but for the variable-$N_{\rm f}$ evolution 
 using Eqs.~(\ref{gsav-eq12}) -- (\ref{gsav-eq15}).
 The corresponding values for the strong coupling $\alpha_{\rm s}(\mu_{\rm r}^2 
 \! =\! 10^4 \mbox{ GeV}^2)$ are given by $0.116461$, $0.116032$ and $0.115663$
 for $\mu_{\rm r}^2 / \mu_{\rm f}^2 = 0.5$, $1$ and $2$, respectively.} 
\label{gsav-t3}
%
\vspace{4mm}
\begin{center}
\begin{tabular}{||c||r|r|r|r|r|r|r|r||}
\hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{NLO, $\, N_{\rm f} = 3\,\ldots 5\,$,
                       $\,\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
\multicolumn{1}{||c||}{$x$} &
\multicolumn{1}{c|} {$xu_v$} &
\multicolumn{1}{c|} {$xd_v$} &
\multicolumn{1}{c|} {$xL_-$} &
\multicolumn{1}{c|} {$xL_+$} &
\multicolumn{1}{c|} {$xs_+$} &
\multicolumn{1}{c|} {$xc_+$} &
\multicolumn{1}{c|} {$xb_+$} &
\multicolumn{1}{c||}{$xg$} \\[0.5mm] \hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.0927^{-4}\!$ &$\! 6.4125^{-5}\!$ &$\! 4.3925^{-6}\!$ &$\! 1.3787^{+2}\!$ &
$\! 6.7857^{+1}\!$ &$\! 6.7139^{+1}\!$ &$\! 6.0071^{+1}\!$ &$\! 1.1167^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 5.5533^{-4}\!$ &$\! 3.2498^{-4}\!$ &$\! 1.9829^{-5}\!$ &$\! 6.9157^{+1}\!$ &
$\! 3.3723^{+1}\!$ &$\! 3.3153^{+1}\!$ &$\! 2.8860^{+1}\!$&$\!5.2289_*^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 2.7419^{-3}\!$ &$\! 1.5989^{-3}\!$ &$\! 8.5701^{-5}\!$ &$\! 3.2996^{+1}\!$ &
$\! 1.5819^{+1}\!$ &$\! 1.5367^{+1}\!$ &$\! 1.2892^{+1}\!$ &$\! 2.2753^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.3039^{-2}\!$ &$\! 7.5664^{-3}\!$ &$\! 3.5582^{-4}\!$ &$\! 1.4822^{+1}\!$ &
$\! 6.8739^{+0}\!$ &$\! 6.5156^{+0}\!$ &$\! 5.1969^{+0}\!$ &$\! 8.9513^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.8507^{-2}\!$ &$\! 3.3652^{-2}\!$ &$\! 1.4329^{-3}\!$ &$\! 6.1772^{+0}\!$ &
$\! 2.6726^{+0}\!$ &$\! 2.3949^{+0}\!$ &$\! 1.7801^{+0}\!$ &$\! 3.0245^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3128^{-1}\!$ &$\! 1.2978^{-1}\!$ &$\! 5.3472^{-3}\!$ &$\! 2.2500^{+0}\!$ &
$\! 8.4161^{-1}\!$ &$\! 6.5235^{-1}\!$ &$\! 4.3894^{-1}\!$ &$\! 7.7491^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.5324^{-1}\!$ &$\! 2.7252^{-1}\!$ &$\! 9.9709^{-3}\!$ &$\! 3.9099^{-1}\!$ &
$\! 1.1425^{-1}\!$ &$\! 6.0071^{-2}\!$ &$\! 3.5441^{-2}\!$ &$\! 8.5586^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.5129^{-1}\!$ &$\! 1.3046^{-1}\!$ &$\! 3.0061^{-3}\!$ &$\! 3.5463^{-2}\!$ &
$\! 9.1084^{-3}\!$ &$\! 3.3595^{-3}\!$ &$\! 1.9039^{-3}\!$ &$\! 7.9625^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.2130^{-1}\!$ &$\! 3.1564^{-2}\!$ &$\! 3.7719^{-4}\!$ &$\! 2.3775^{-3}\!$ &
$\! 5.7606^{-4}\!$ &$\! 1.6761^{-4}\!$ &$\! 1.0021^{-4}\!$ &$\! 7.7265^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.0101^{-2}\!$&$\! 3.0932^{-3}\!$ &$\! 1.3440^{-5}\!$ &$\! 5.2605^{-5}\!$ &
$\! 1.2166^{-5}\!$ &$\! 2.7408^{-6}\!$ &$\! 2.0095^{-6}\!$ &$\! 3.7574^{-4}\!$\\
$\! 0.9    \!$ &
$\!3.5232_*^{-4}\!$&$\! 1.7855^{-5}\!$ &$\! 8.680 ^{-9}\!$ &$\! 2.028 ^{-8}\!$ &
$\! 3.896 ^{-9}\!$ &$\!-2.666^{-10}\!$ &$\! 5.819^{-10}\!$ &$\! 1.1954^{-6}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = 2\, \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 9.5154^{-5}\!$ &$\! 5.5970^{-5}\!$ &$\! 3.4869^{-6}\!$ &$\! 1.2301^{+2}\!$ &
$\! 6.0424^{+1}\!$ &$\! 5.9703^{+1}\!$ &$\! 5.3916^{+1}\!$ &$\! 1.0568^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 4.9433^{-4}\!$ &$\! 2.8998^{-4}\!$ &$\! 1.6229^{-5}\!$ &$\! 6.2149^{+1}\!$ &
$\! 3.0215^{+1}\!$ &$\! 2.9643^{+1}\!$ &$\! 2.6100^{+1}\!$ &$\! 4.9744^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 2.4974^{-3}\!$ &$\! 1.4600^{-3}\!$ &$\! 7.2776^{-5}\!$ &$\! 2.9936^{+1}\!$ &
$\! 1.4286^{+1}\!$ &$\! 1.3831^{+1}\!$ &$\! 1.1768^{+1}\!$ &$\! 2.1783^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.2161^{-2}\!$ &$\! 7.0751^{-3}\!$ &$\! 3.1597^{-4}\!$ &$\! 1.3631^{+1}\!$ &
$\! 6.2759^{+0}\!$ &$\! 5.9160^{+0}\!$ &$\! 4.7997^{+0}\!$ &$\! 8.6372^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.5907^{-2}\!$ &$\! 3.2239^{-2}\!$ &$\!1.3377_*^{-3}\!$&$\! 5.8020^{+0}\!$ &
$\! 2.4824^{+0}\!$ &$\! 2.2029^{+0}\!$ &$\! 1.6701^{+0}\!$ &$\! 2.9488^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.2670^{-1}\!$ &$\! 1.2756^{-1}\!$ &$\! 5.2304^{-3}\!$ &$\! 2.1840^{+0}\!$ &
$\! 8.0522^{-1}\!$ &$\! 6.1376^{-1}\!$ &$\! 4.2146^{-1}\!$ &$\! 7.6713^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.5915^{-1}\!$ &$\! 2.7636^{-1}\!$ &$\! 1.0151^{-2}\!$ &$\! 3.9744^{-1}\!$ &
$\! 1.1509^{-1}\!$ &$\! 5.9318^{-2}\!$ &$\! 3.5492^{-2}\!$ &$\! 8.7075^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.6281^{-1}\!$ &$\! 1.3530^{-1}\!$ &$\! 3.1367^{-3}\!$ &$\! 3.7201^{-2}\!$ &
$\! 9.5529^{-3}\!$ &$\! 3.5209^{-3}\!$ &$\! 1.9869^{-3}\!$ &$\! 8.2697^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.2739^{-1}\!$ &$\! 3.3311^{-2}\!$ &$\! 4.0102^{-4}\!$ &$\! 2.5552^{-3}\!$ &
$\! 6.2392^{-4}\!$ &$\! 1.8815^{-4}\!$ &$\! 1.0842^{-4}\!$ &$\! 8.1131^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.1534^{-2}\!$ &$\! 3.3317^{-3}\!$ &$\! 1.4600^{-5}\!$ &$\! 5.8333^{-5}\!$ &
$\! 1.3778^{-5}\!$ &$\! 3.5193^{-6}\!$ &$\! 2.2960^{-6}\!$ &$\! 3.9788^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.9150^{-4}\!$ &$\! 1.9963^{-5}\!$ &$\! 9.797 ^{-9}\!$ &$\! 2.502 ^{-8}\!$ &
$\! 5.459 ^{-9}\!$ &$\! 7.574^{-10}\!$ &$\! 1.002 ^{-9}\!$ &$\! 1.2895^{-6}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = 1/2\, \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.2937^{-4}\!$ &$\! 7.5695^{-5}\!$ &$\! 5.8161^{-6}\!$ &$\! 1.4923^{+2}\!$ &
$\! 7.3543^{+1}\!$ &$\! 7.2829^{+1}\!$ &$\! 6.4225^{+1}\!$ &$\! 1.1239^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 6.3890^{-4}\!$ &$\! 3.7272^{-4}\!$ &$\! 2.5317^{-5}\!$ &$\! 7.4973^{+1}\!$ &
$\! 3.6635^{+1}\!$ &$\! 3.6068^{+1}\!$ &$\! 3.0846^{+1}\!$ &$\! 5.2763^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 3.0615^{-3}\!$ &$\! 1.7796^{-3}\!$ &$\! 1.0470^{-4}\!$ &$\! 3.5762^{+1}\!$ &
$\! 1.7205^{+1}\!$ &$\! 1.6755^{+1}\!$ &$\! 1.3756^{+1}\!$ &$\! 2.3006^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.4120^{-2}\!$ &$\! 8.1674^{-3}\!$ &$\! 4.1155^{-4}\!$ &$\! 1.5994^{+1}\!$ &
$\! 7.4628^{+0}\!$ &$\! 7.1064^{+0}\!$ &$\!5.5232_*^{+0}\!$&$\! 9.0599^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 6.1458^{-2}\!$ &$\! 3.5239^{-2}\!$ &$\! 1.5557^{-3}\!$ &$\! 6.5774^{+0}\!$ &
$\! 2.8758^{+0}\!$ &$\! 2.5999^{+0}\!$ &$\! 1.8773^{+0}\!$ &$\! 3.0589^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3574^{-1}\!$ &$\! 1.3187^{-1}\!$ &$\! 5.4739^{-3}\!$ &$\! 2.3254^{+0}\!$ &
$\! 8.8301^{-1}\!$ &$\! 6.9613^{-1}\!$ &$\! 4.5587^{-1}\!$ &$\! 7.7996^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.4630^{-1}\!$ &$\! 2.6827^{-1}\!$ &$\! 9.7828^{-3}\!$ &$\! 3.8485^{-1}\!$ &
$\! 1.1348^{-1}\!$ &$\! 6.0825^{-2}\!$ &$\! 3.5524^{-2}\!$ &$\! 8.4801^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.4035^{-1}\!$ &$\! 1.2596^{-1}\!$ &$\! 2.8877^{-3}\!$ &$\! 3.3752^{-2}\!$ &
$\! 8.6313^{-3}\!$ &$\! 3.1345^{-3}\!$ &$\! 1.8308^{-3}\!$ &$\! 7.7817^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.1589^{-1}\!$ &$\! 3.0034^{-2}\!$ &$\! 3.5665^{-4}\!$ &$\! 2.2032^{-3}\!$ &
$\! 5.2398^{-4}\!$ &$\! 1.3891^{-4}\!$ &$\! 9.2307^{-5}\!$ &$\! 7.4943^{-3}\!$\\
$\! 0.7    \!$ &
$\! 1.8886^{-2}\!$ &$\! 2.8923^{-3}\!$ &$\! 1.2475^{-5}\!$ &$\! 4.6898^{-5}\!$ &
$\! 1.0350^{-5}\!$ &$\! 1.6179^{-6}\!$ &$\! 1.7069^{-6}\!$ &$\! 3.6221^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.2053^{-4}\!$ &$\! 1.6155^{-5}\!$ &$\! 7.787 ^{-9}\!$ &$\! 1.489 ^{-8}\!$ &
$\! 1.848 ^{-9}\!$ &$\!-1.884 ^{-9}\!$ &$\! 6.129^{-11}\!$ &$\! 1.1353^{-6}\!$\\
%
\hline \hline
%
\end{tabular}
\end{center}
\end{table}

\begin{table}[htp]
\caption{Reference results for the $N_{\rm f}\!=\!4$ next-next-to-leading-order
 evolution for the initial conditions (\ref{gsav-eq8}) -- (\ref{gsav-eq10}).
 The corresponding value of the strong coupling is $\alpha_{\rm s}
 (\mu_{\rm r}^2 \! =\! 10^4 \mbox{ GeV}^2) = 0.110141$. The valence 
 distributions $s_v$ and $c_v$ are equal for the input (\ref{gsav-eq9}). The 
 notation is explained in the first two paragraphs of Section \ref{gsav-s3}.}
\label{gsav-t4}
%
\vspace{4mm}
\begin{center}
\begin{tabular}{||c||r|r|r|r|r|r|r|r||}
\hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{NNLO, $\, N_{\rm f} = 4\,$, 
                       $\,\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
\multicolumn{1}{||c||}{$x$} &
\multicolumn{1}{c|} {$xu_v$} &
\multicolumn{1}{c|} {$xd_v$} &
\multicolumn{1}{c|} {$xL_-$} &
\multicolumn{1}{c|} {$xL_+$} &
\multicolumn{1}{c|} {$xs_v$} &
\multicolumn{1}{c|} {$xs_+$} &
\multicolumn{1}{c|} {$xc_+$} &
\multicolumn{1}{c||}{$xg$} \\[0.5mm] \hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.4425^{-4}\!$ &$\! 8.9516_*^{-5}\!$&$\!6.0836^{-6}\!$ &$\! 1.3492^{+2}\!$ &
$\! 1.3116^{-5}\!$ &$\! 6.6385^{+1}\!$ &$\! 6.5669^{+1}\!$ &$\! 1.0105^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 6.7325^{-4}\!$ &$\! 4.1192^{-4}\!$ &$\! 2.6060^{-5}\!$ &$\! 6.9636^{+1}\!$ &
$\! 4.8617^{-5}\!$ &$\! 3.3965^{+1}\!$ &$\! 3.3396^{+1}\!$ &$\! 4.9543^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 3.0782^{-3}\!$ &$\! 1.8510^{-3}\!$ &$\! 1.0652^{-4}\!$ &$\! 3.3826^{+1}\!$ &
$\! 1.5280^{-4}\!$ &$\! 1.6236^{+1}\!$ &$\! 1.5784^{+1}\!$ &$\! 2.2406^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.3746^{-2}\!$ &$\! 8.0979^{-3}\!$ &$\! 4.1608^{-4}\!$ &$\! 1.5278^{+1}\!$ &
$\! 3.3973^{-4}\!$ &$\! 7.1041^{+0}\!$ &$\! 6.7468^{+0}\!$&$\!9.0854_*^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.9222^{-2}\!$ &$\! 3.4137^{-2}\!$ &$\! 1.5718^{-3}\!$ &$\! 6.3249^{+0}\!$ &
$\! 2.6800^{-4}\!$ &$\! 2.7483^{+0}\!$ &$\! 2.4713^{+0}\!$ &$\! 3.1353^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3074^{-1}\!$ &$\! 1.2920^{-1}\!$ &$\! 5.5461^{-3}\!$ &$\! 2.2736^{+0}\!$ &
$\!-5.2053^{-4}\!$ &$\! 8.5420^{-1}\!$ &$\! 6.6542^{-1}\!$ &$\! 8.1341^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.5178^{-1}\!$ &$\! 2.7165^{-1}\!$ &$\! 1.0024^{-2}\!$ &$\! 3.9020^{-1}\!$ &
$\!-3.0680^{-4}\!$ &$\! 1.1386^{-1}\!$ &$\! 5.9780^{-2}\!$ &$\! 9.0567^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.5071^{-1}\!$ &$\! 1.3025^{-1}\!$ &$\! 3.0098^{-3}\!$ &$\! 3.5358^{-2}\!$ &
$\!-3.1905^{-5}\!$ &$\! 9.0479^{-3}\!$ &$\! 3.3060^{-3}\!$ &$\! 8.4184^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.2117^{-1}\!$ &$\! 3.1528^{-2}\!$ &$\! 3.7743^{-4}\!$ &$\! 2.3866^{-3}\!$ &
$\!-2.7207^{-6}\!$ &$\! 5.7966^{-4}\!$ &$\! 1.7171^{-4}\!$&$\!8.1127_*^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.0078^{-2}\!$ &$\! 3.0886^{-3}\!$ &$\! 1.3442^{-5}\!$ &$\! 5.4226^{-5}\!$ &
$\!-1.0121^{-7}\!$ &$\! 1.2936^{-5}\!$ &$\! 3.5305^{-6}\!$ &$\! 3.8946^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.5111^{-4}\!$ &$\! 1.7783^{-5}\!$ &$\! 8.870 ^{-9}\!$ &$\! 2.630 ^{-8}\!$ &
$\!-1.450^{-10}\!$ &$\! 7.115 ^{-9}\!$ &$\! 2.973 ^{-9}\!$ &$\! 1.2147^{-6}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = 2\, \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.2819^{-4}\!$ &$\! 7.8560^{-5}\!$ &$\! 5.2370^{-6}\!$ &$\! 1.3229^{+2}\!$ &
$\! 8.8691^{-6}\!$ &$\! 6.5070^{+1}\!$ &$\! 6.4352^{+1}\!$ &$\! 1.0296^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 6.1508^{-4}\!$ &$\! 3.7250^{-4}\!$ &$\! 2.2897^{-5}\!$ &$\! 6.7704^{+1}\!$ &
$\! 3.3358^{-5}\!$ &$\! 3.2997^{+1}\!$ &$\! 3.2427^{+1}\!$ &$\! 4.9907^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 2.8891^{-3}\!$ &$\! 1.7248^{-3}\!$ &$\! 9.5768^{-5}\!$ &$\! 3.2745^{+1}\!$ &
$\! 1.0595^{-4}\!$ &$\! 1.5694^{+1}\!$ &$\! 1.5241^{+1}\!$ &$\! 2.2382^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.3223^{-2}\!$ &$\! 7.7642^{-3}\!$ &$\! 3.8437^{-4}\!$ &$\! 1.4803^{+1}\!$ &
$\! 2.3744^{-4}\!$ &$\! 6.8647^{+0}\!$ &$\! 6.5065^{+0}\!$ &$\! 9.0293^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.8091^{-2}\!$ &$\! 3.3494^{-2}\!$ &$\! 1.4978^{-3}\!$ &$\! 6.1716^{+0}\!$ &
$\! 1.8990^{-4}\!$ &$\! 2.6702^{+0}\!$ &$\! 2.3924^{+0}\!$ &$\! 3.1117^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.2927^{-1}\!$ &$\! 1.2858^{-1}\!$ &$\! 5.4465_*^{-3}\!$&$\!2.2481^{+0}\!$ &
$\!-3.6267^{-4}\!$ &$\! 8.4001^{-1}\!$ &$\! 6.5031^{-1}\!$ &$\! 8.0964^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.5429^{-1}\!$ &$\! 2.7326^{-1}\!$ &$\! 1.0073^{-2}\!$ &$\! 3.9298^{-1}\!$ &
$\!-2.1631^{-4}\!$ &$\! 1.1440^{-1}\!$ &$\! 5.9718^{-2}\!$ &$\! 9.0854^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.5501^{-1}\!$ &$\! 1.3205^{-1}\!$ &$\! 3.0556^{-3}\!$ &$\! 3.6008^{-2}\!$ &
$\!-2.2642^{-5}\!$ &$\! 9.2226^{-3}\!$ &$\! 3.3770^{-3}\!$&$\!8.5020_*^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.2340^{-1}\!$ &$\! 3.2166^{-2}\!$ &$\! 3.8590^{-4}\!$ &$\! 2.4458^{-3}\!$ &
$\!-1.9414^{-6}\!$ &$\! 5.9487^{-4}\!$ &$\! 1.7699^{-4}\!$&$\!8.2294_*^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.0597^{-2}\!$ &$\! 3.1751^{-3}\!$ &$\! 1.3854^{-5}\!$ &$\! 5.5708^{-5}\!$ &
$\!-7.2724^{-8}\!$ &$\! 1.3244^{-5}\!$ &$\! 3.5361^{-6}\!$ &$\! 3.9686^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.6527^{-4}\!$ &$\! 1.8544^{-5}\!$ &$\! 9.204 ^{-9}\!$ &$\! 2.616 ^{-8}\!$ &
$\!-1.057^{-10}\!$ &$\! 6.700 ^{-9}\!$ &$\! 2.364 ^{-9}\!$ &$\! 1.2497^{-6}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = 1/2\, \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.6566^{-4}\!$ &$\! 1.0501^{-4}\!$ &$\! 7.1003^{-6}\!$ &$\! 1.3181^{+2}\!$ &
$\! 2.0996^{-5}\!$ &$\! 6.4832^{+1}\!$ &$\! 6.4117^{+1}\!$ &$\! 9.4582^{+2}\!$\\
$\! 10^{-6}\!$ &
$\! 7.4521^{-4}\!$ &$\! 4.6428^{-4}\!$ &$\! 2.9744^{-5}\!$ &$\! 6.9457^{+1}\!$ &
$\! 7.6503^{-5}\!$ &$\! 3.3877^{+1}\!$&$\!3.3309_*^{+1}\!$&$\!4.7543_*^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 3.2879^{-3}\!$ &$\! 2.0041^{-3}\!$ &$\! 1.1872^{-4}\!$ &$\! 3.4193^{+1}\!$ &
$\! 2.3809^{-4}\!$ &$\! 1.6421^{+1}\!$ &$\! 1.5970^{+1}\!$ &$\! 2.1922^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.4223^{-2}\!$ &$\! 8.4396^{-3}\!$ &$\! 4.5146^{-4}\!$ &$\! 1.5524^{+1}\!$ &
$\! 5.2716^{-4}\!$ &$\! 7.2281^{+0}\!$ &$\! 6.8715^{+0}\!$ &$\! 9.0112^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.9889^{-2}\!$ &$\! 3.4553^{-2}\!$ &$\! 1.6545^{-3}\!$ &$\! 6.4073^{+0}\!$ &
$\! 4.1428^{-4}\!$ &$\! 2.7908^{+0}\!$ &$\! 2.5144^{+0}\!$ &$\! 3.1340^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3100^{-1}\!$ &$\! 1.2915_*^{-1}\!$&$\!5.6663^{-3}\!$ &$\! 2.2851^{+0}\!$ &
$\!-8.0979^{-4}\!$ &$\! 8.6079^{-1}\!$ &$\! 6.7251^{-1}\!$ &$\! 8.1526^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.5039^{-1}\!$ &$\! 2.7076^{-1}\!$ &$\! 1.0032^{-2}\!$&$\!3.8851_*^{-1}\!$ &
$\!-4.7547^{-4}\!$ &$\! 1.1333^{-1}\!$ &$\! 5.9500^{-2}\!$ &$\! 9.0801^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.4890^{-1}\!$ &$\! 1.2949^{-1}\!$ &$\! 2.9942^{-3}\!$ &$\! 3.5090^{-2}\!$ &
$\!-4.9326^{-5}\!$ &$\! 8.9666^{-3}\!$ &$\! 3.2669^{-3}\!$ &$\! 8.4307^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.2026^{-1}\!$ &$\! 3.1269^{-2}\!$&$\!3.7428_*^{-4}\!$ &$\! 2.3728^{-3}\!$ &
$\!-4.1969^{-6}\!$ &$\! 5.7784^{-4}\!$ &$\! 1.7391^{-4}\!$ &$\! 8.1102^{-3}\!$\\
$\! 0.7    \!$ &
$\! 1.9867^{-2}\!$ &$\! 3.0534^{-3}\!$ &$\! 1.3285^{-5}\!$ &$\! 5.4607^{-5}\!$ &
$\!-1.5564^{-7}\!$&$\!1.3275_*^{-5}\!$ &$\! 3.9931^{-6}\!$&$\!3.8822_*^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.4524^{-4}\!$ &$\! 1.7466^{-5}\!$ &$\! 8.828 ^{-9}\!$ &$\! 2.930 ^{-8}\!$ &
$\!-2.216^{-10}\!$ &$\! 8.837 ^{-9}\!$ &$\! 4.777 ^{-9}\!$ &$\! 1.2041^{-6}\!$\\
%
\hline \hline
%
\end{tabular}
\end{center}
\end{table}

\begin{table}[htp]
\caption{As Table \ref{gsav-t4}, but for the variable-$N_{\rm f}$ evolution
 using Eqs.~(\ref{gsav-eq12}) -- (\ref{gsav-eq15}).
 The corresponding values for the strong coupling $\alpha_{\rm s}(\mu_{\rm r}^2
 \! =\! 10^4 \mbox{ GeV}^2)$ are given by $0.115818$, $0.115605$ and $0.115410$
 for $\mu_{\rm r}^2 / \mu_{\rm f}^2 = 0.5$, $1$ and $2$, respectively. For
 brevity the small, but non-vanishing valence distributions $s_v$, $c_v$ and 
 $b_v$ are not displayed.}
\label{gsav-t5}
%
\vspace{3mm}
\begin{center}
\begin{tabular}{||c||r|r|r|r|r|r|r|r||}
\hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{NNLO, $\, N_{\rm f} = 3\,\ldots 5\,$,
                       $\,\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
\multicolumn{1}{||c||}{$x$} &
\multicolumn{1}{c|} {$xu_v$} &
\multicolumn{1}{c|} {$xd_v$} &
\multicolumn{1}{c|} {$xL_-$} &
\multicolumn{1}{c|} {$xL_+$} &
\multicolumn{1}{c|} {$xs_+$} &
\multicolumn{1}{c|} {$xc_+$} &
\multicolumn{1}{c|} {$xb_+$} &
\multicolumn{1}{c||}{$xg$} \\[0.5mm] \hline \hline
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.5053^{-4}\!$ &$\! 9.3396^{-5}\!$ &$\! 6.3350^{-6}\!$ &$\! 1.4236^{+2}\!$ &
$\! 7.0107^{+1}\!$ &$\! 6.8526^{+1}\!$ &$\! 5.8645^{+1}\!$ &$\! 1.0101^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 6.9800^{-4}\!$ &$\! 4.2692^{-4}\!$ &$\! 2.7015^{-5}\!$ &$\! 7.2905^{+1}\!$ &
$\! 3.5601^{+1}\!$ &$\! 3.4447^{+1}\!$ &$\! 2.8749^{+1}\!$ &$\! 4.9242^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 3.1691^{-3}\!$ &$\! 1.9048^{-3}\!$ &$\! 1.0987^{-4}\!$ &$\! 3.5102^{+1}\!$ &
$\! 1.6875^{+1}\!$ &$\! 1.6059^{+1}\!$ &$\! 1.3011^{+1}\!$ &$\! 2.2122^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.4048^{-2}\!$ &$\! 8.2711^{-3}\!$ &$\! 4.2655^{-4}\!$ &$\! 1.5694^{+1}\!$ &
$\! 7.3126^{+0}\!$ &$\! 6.7623^{+0}\!$ &$\! 5.2796^{+0}\!$ &$\! 8.8987^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 6.0085^{-2}\!$ &$\! 3.4610^{-2}\!$ &$\! 1.5988^{-3}\!$ &$\! 6.4190^{+0}\!$ &
$\! 2.7963^{+0}\!$ &$\! 2.4502^{+0}\!$ &$\! 1.8138^{+0}\!$ &$\! 3.0407^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3240^{-1}\!$ &$\! 1.3001^{-1}\!$ &$\! 5.5869^{-3}\!$ &$\! 2.2762^{+0}\!$ &
$\! 8.5667^{-1}\!$ &$\! 6.6665^{-1}\!$ &$\! 4.5032^{-1}\!$ &$\! 7.7873^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.4993^{-1}\!$ &$\! 2.7035^{-1}\!$ &$\! 9.9604^{-3}\!$ &$\! 3.8527^{-1}\!$ &
$\! 1.1231^{-1}\!$ &$\! 6.4471^{-2}\!$ &$\! 3.7283^{-2}\!$ &$\! 8.5270^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.4622^{-1}\!$ &$\! 1.2833^{-1}\!$ &$\! 2.9571^{-3}\!$ &$\! 3.4600^{-2}\!$ &
$\! 8.8409^{-3}\!$ &$\! 4.0131^{-3}\!$ &$\! 2.1046^{-3}\!$ &$\! 7.8896^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.1868^{-1}\!$ &$\! 3.0811^{-2}\!$ &$\! 3.6760^{-4}\!$ &$\! 2.3198^{-3}\!$ &
$\! 5.6309^{-4}\!$ &$\! 2.3748^{-4}\!$ &$\! 1.2003^{-4}\!$ &$\! 7.6399^{-3}\!$\\
$\! 0.7    \!$ &
$\! 1.9486^{-2}\!$ &$\! 2.9902^{-3}\!$ &$\! 1.2964^{-5}\!$ &$\! 5.2335^{-5}\!$ &
$\! 1.2504^{-5}\!$ &$\! 5.6019^{-6}\!$ &$\! 2.8883^{-6}\!$ &$\! 3.7079^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.3522^{-4}\!$ &$\! 1.6933^{-5}\!$ &$\! 8.4201^{-9}\!$ &$\! 2.5098^{-8}\!$ &
$\! 6.8378^{-9}\!$ &$\! 4.3190^{-9}\!$ &$\! 2.6761^{-9}\!$ &$\! 1.1731^{-6}\!$\\
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = 2\, \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.3305^{-4}\!$ &$\! 8.1537^{-5}\!$ &$\! 5.4285^{-6}\!$ &$\! 1.3770^{+2}\!$ &
$\! 6.7776^{+1}\!$ &$\! 6.6452^{+1}\!$ &$\! 5.7708^{+1}\!$ &$\! 1.0172^{+3}\!$\\
$\! 10^{-6}\!$ &
$\! 6.3465^{-4}\!$ &$\! 3.8427^{-4}\!$ &$\! 2.3638^{-5}\!$ &$\! 7.0043^{+1}\!$ &
$\! 3.4167^{+1}\!$ &$\! 3.3187^{+1}\!$ &$\! 2.8117^{+1}\!$ &$\! 4.9102^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 2.9622^{-3}\!$ &$\! 1.7679^{-3}\!$ &$\! 9.8415^{-5}\!$ &$\! 3.3642^{+1}\!$ &
$\! 1.6144^{+1}\!$ &$\! 1.5435^{+1}\!$ &$\! 1.2685^{+1}\!$ &$\! 2.1913^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.3469^{-2}\!$ &$\! 7.9047^{-3}\!$ &$\! 3.9275^{-4}\!$ &$\! 1.5088^{+1}\!$ &
$\! 7.0083^{+0}\!$ &$\! 6.5145^{+0}\!$ &$\! 5.1475^{+0}\!$ &$\! 8.7878^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 5.8793^{-2}\!$ &$\! 3.3879^{-2}\!$ &$\! 1.5195^{-3}\!$ &$\! 6.2333^{+0}\!$ &
$\! 2.7018^{+0}\!$ &$\! 2.3753^{+0}\!$ &$\! 1.7742^{+0}\!$ &$\! 3.0063^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3060^{-1}\!$ &$\! 1.2923^{-1}\!$ &$\! 5.4792^{-3}\!$ &$\! 2.2478^{+0}\!$ &
$\! 8.4084^{-1}\!$ &$\! 6.5026^{-1}\!$ &$\! 4.4324^{-1}\!$ &$\! 7.7467^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.5280^{-1}\!$ &$\! 2.7222^{-1}\!$ &$\!1.0021_*^{-2}\!$&$\! 3.8898^{-1}\!$ &
$\! 1.1313^{-1}\!$ &$\! 6.2921^{-2}\!$ &$\! 3.7050^{-2}\!$ &$\! 8.5900^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.5141^{-1}\!$ &$\! 1.3051^{-1}\!$ &$\! 3.0134^{-3}\!$ &$\! 3.5398^{-2}\!$ &
$\! 9.0559^{-3}\!$ &$\! 3.8724^{-3}\!$ &$\! 2.0992^{-3}\!$ &$\! 8.0225^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.2140^{-1}\!$ &$\! 3.1590^{-2}\!$ &$\! 3.7799^{-4}\!$ &$\! 2.3919^{-3}\!$ &
$\! 5.8148^{-4}\!$ &$\! 2.2373^{-4}\!$ &$\! 1.1917^{-4}\!$ &$\! 7.8099^{-3}\!$\\
$\! 0.7    \!$ &
$\! 2.0120^{-2}\!$ &$\! 3.0955^{-3}\!$ &$\! 1.3467^{-5}\!$ &$\! 5.4182^{-5}\!$ &
$\! 1.2896^{-5}\!$ &$\! 5.0316^{-6}\!$ &$\! 2.8148^{-6}\!$ &$\! 3.8098^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.5230^{-4}\!$ &$\! 1.7849^{-5}\!$ &$\! 8.835 ^{-9}\!$ &$\! 2.523 ^{-8}\!$ &
$\! 6.499 ^{-9}\!$ &$\! 3.377 ^{-9}\!$ &$\! 2.404 ^{-9}\!$ &$\! 1.2196^{-6}\!$\\
%
\hline \hline
%
\multicolumn{9}{||c||}{} \\[-3mm]
\multicolumn{9}{||c||}{$\mu_{\rm r}^2 = 1/2\, \mu_{\rm f}^2$} \\
\multicolumn{9}{||c||}{} \\[-0.3cm] \hline \hline
 & & & & & & & \\[-0.3cm]
$\! 10^{-7}\!$ &
$\! 1.7462^{-4}\!$ &$\! 1.1066^{-4}\!$ &$\! 7.4685^{-6}\!$ &$\! 1.4238^{+2}\!$ &
$\! 7.0121^{+1}\!$ &$\! 6.8041^{+1}\!$ &$\! 5.7115^{+1}\!$ &$\! 9.6557^{+2}\!$\\
$\! 10^{-6}\!$ &
$\! 7.7945^{-4}\!$ &$\! 4.8542^{-4}\!$ &$\!3.1111_*^{-5}\!$&$\! 7.4239^{+1}\!$ &
$\! 3.6270^{+1}\!$ &$\! 3.4778^{+1}\!$ &$\! 2.8492^{+1}\!$ &$\! 4.8135^{+2}\!$\\
$\! 10^{-5}\!$ &
$\! 3.4099^{-3}\!$ &$\! 2.0775^{-3}\!$ &$\! 1.2343^{-4}\!$ &$\! 3.6106^{+1}\!$ &
$\! 1.7379^{+1}\!$ &$\! 1.6353^{+1}\!$ &$\! 1.3052^{+1}\!$ &$\! 2.1979^{+2}\!$\\
$\! 10^{-4}\!$ &
$\! 1.4619^{-2}\!$ &$\! 8.6691^{-3}\!$ &$\! 4.6591^{-4}\!$ &$\! 1.6161^{+1}\!$ &
$\! 7.5479^{+0}\!$ &$\! 6.8868^{+0}\!$ &$\! 5.3331^{+0}\!$ &$\! 8.9300^{+1}\!$\\
$\! 10^{-3}\!$ &
$\! 6.1014^{-2}\!$ &$\! 3.5170^{-2}\!$ &$\! 1.6915^{-3}\!$ &$\! 6.5570^{+0}\!$ &
$\! 2.8669^{+0}\!$ &$\! 2.4816^{+0}\!$ &$\! 1.8360^{+0}\!$&$\!3.0619_*^{+1}\!$\\
$\! 10^{-2}\!$ &
$\! 2.3320^{-1}\!$ &$\! 1.3024^{-1}\!$ &$\! 5.7229^{-3}\!$ &$\! 2.2924^{+0}\!$ &
$\! 8.6595^{-1}\!$ &$\! 6.7561^{-1}\!$ &$\! 4.5538^{-1}\!$ &$\! 7.8182^{+0}\!$\\
$\! 0.1    \!$ &
$\! 5.4799^{-1}\!$ &$\! 2.6905^{-1}\!$ &$\!9.9481_*^{-3}\!$&$\! 3.8193^{-1}\!$ &
$\! 1.1125^{-1}\!$ &$\! 6.7100^{-2}\!$ &$\! 3.7702^{-2}\!$ &$\! 8.4914^{-1}\!$\\
$\! 0.3    \!$ &
$\! 3.4291^{-1}\!$ &$\! 1.2693^{-1}\!$ &$\! 2.9238^{-3}\!$ &$\! 3.4069^{-2}\!$ &
$\! 8.6866^{-3}\!$ &$\! 4.3918^{-3}\!$ &$\! 2.1434^{-3}\!$&$\!7.8107_*^{-2}\!$\\
$\! 0.5    \!$ &
$\! 1.1694^{-1}\!$ &$\! 3.0310^{-2}\!$ &$\! 3.6112^{-4}\!$ &$\! 2.2827^{-3}\!$ &
$\! 5.5538^{-4}\!$ &$\! 2.7737^{-4}\!$ &$\! 1.2414^{-4}\!$ &$\! 7.5373^{-3}\!$\\
$\! 0.7    \!$ &
$\! 1.9076^{-2}\!$ &$\! 2.9217^{-3}\!$ &$\! 1.2646^{-5}\!$ &$\! 5.2035^{-5}\!$ &
$\! 1.2677^{-5}\!$ &$\! 7.2055^{-6}\!$ &$\! 3.0900^{-6}\!$&$\!3.6439_*^{-4}\!$\\
$\! 0.9    \!$ &
$\! 3.2404^{-4}\!$ &$\! 1.6333^{-5}\!$ &$\! 8.231 ^{-9}\!$ &$\! 2.752 ^{-8}\!$ &
$\! 8.380 ^{-9}\!$ &$\! 6.769 ^{-9}\!$ &$\! 3.204 ^{-9}\!$ & $\!1.1424^{-6}\!$\\
%
\hline \hline
%
\end{tabular}
\end{center}
\end{table}

%\bibliography{pdfref}
%\end{document}








