  Aims. We use Hermite splines to interpolate pressure and its derivatives
simultaneously, thereby preserving mathematical relations between the
derivatives. The method therefore guarantees that thermodynamic identities are
obeyed even between mesh points. In addition, our method enables an estimation
of the precision of the interpolation by comparing the Hermite-spline results
with those of frequent cubic (B-) spline interpolation.
  Methods. We have interpolated pressure as a function of temperature and
density with quintic Hermite 2D-splines. The Hermite interpolation requires
knowledge of pressure and its first and second derivatives at every mesh point.
To obtain the partial derivatives at the mesh points, we used tabulated values
if given or else thermodynamic equalities, or, if not available, values
obtained by differentiating B-splines.
  Results. The results were obtained with the grid of the SAHA-S
equation-of-state (EOS) tables. The maximum $lg P$ difference lies in the range
from $10^{-9}$ to $10^{-4}$, and $\Gamma_1$ difference varies from $10^{-9}$ to
$10^{-3}$. Specifically, for the points of a solar model, the maximum
differences are one order of magnitude smaller than the aforementioned values.
The poorest precision is found in the dissociation and ionization regions,
occurring at $T \sim 1.5\cdot 10^3 - 10^5$ K. The best precision is achieved at
higher temperatures, $T>10^5$ K. To discuss the significance of the
interpolation errors we compare them with the corresponding difference between
two different equation-of-state formalisms, SAHA-S and OPAL 2005. We find that
the interpolation errors of the pressure are a few orders of magnitude less
than the differences from between the physical formalisms, which is
particularly true for the solar-model points.
