  We present the bivariate neutral atomic hydrogen (HI)---stellar mass function
(HISMF) (phi(M_HI, M_*)) for massive (log M_*/M_sun > 10) galaxies derived from
a sample of 480 local (0.025 < z < 0.050) galaxies observed in HI at Arecibo as
part of the GALEX Arecibo SDSS Survey (GASS). We fit six different models to
the HISMF and find that a Schechter function that extends down to a 1% HI gas
fraction, with an additional fractional contribution below that limit, is the
best parametrization of the HISMF. We calculate Omega_{HI, M_* >10^10} and find
that massive galaxies contribute 41% of the HI density in the local universe.
In addition to the binned HISMF we derive a continuous bivariate fit, which
reveals that the Schechter parameters only vary weakly with stellar mass:
M_HI^*, the characteristic HI mass, scales as M_*^0.39, alpha, the slope of the
HISMF at moderate HI masses, scales as M_*^0.07, and f, the fraction of
galaxies with HI gas fraction greater than 1%, scales as M_*^-0.24. The
variation of f with stellar mass should be a strong constraint for numerical
simulations. To understand the physical mechanisms that produce the shape of
the HISMF we redefine the parameters of the Schechter function as explicit
functions of stellar mass and star formation rate to produce a trivariate fit.
This analysis reveals strong trends with SFR. While M_HI^* varies weakly with
stellar mass and SFR, alpha is a stronger function of both stellar mass and
especially star formation rate. The HISMF is a crucial tool that can be used to
constrain cosmological galaxy simulations, test observational predictions of
the HI content of populations of galaxies, and identify galaxies whose
properties deviate from average trends.
