We present the bivariate neutral atomic hydrogen (H I)-stellar mass function (HISMF) φ(MHI, M*) for massive (log M */M⊙> 10) galaxies derived from a sample of 480 local (0.025 <z <0.050) galaxies observed in H I at Arecibo as part of the GALEX Arecibo SDSS Survey. We fit six different models to the HISMF and find that a Schechter function that extends down to a 1% H I gas fraction, with an additional fractional contribution below that limit, is the best parameterization of the HISMF. We calculate ωHIM*> 1010 and find that massive galaxies contribute 41% of the H I 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 M0.39 *; α, the slope of the HISMF at moderate HI masses, scales as M0.07; and f, the fraction of galaxies with H I gas fraction greater than 1%, scales as . 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 (SFR) to produce a trivariate fit. This analysis reveals strong trends with SFR. While varies weakly with stellar mass and SFR (M* HI∞M0.22 *, M* HI∞ SFR -0.03), α is a stronger function of both stellar mass and especially SFR (α ∞M* 0.47, α ∞ SFR0.95). The HISMF is a crucial tool that can be used to constrain cosmological galaxy simulations, test observational predictions of the H I content of populations of galaxies, and identify galaxies whose properties deviate from average trends. © 2013. The American Astronomical Society. All rights reserved.