  A collection of algorithms is described for numerically computing with smooth
functions defined on the unit sphere. Functions are approximated to essentially
machine precision by using a structure-preserving iterative variant of Gaussian
elimination together with the double Fourier sphere method. We show that this
procedure allows for stable differentiation, reduces the oversampling of
functions near the poles, and converges for certain analytic functions.
Operations such as function evaluation, differentiation, and integration are
particularly efficient and can be computed by essentially one-dimensional
algorithms. A highlight is an optimal complexity direct solver for Poisson's
equation on the sphere using a spectral method. Without parallelization, we
solve Poisson's equation with $100$ million degrees of freedom in one minute on
a standard laptop. Numerical results are presented throughout. In a companion
paper (part II) we extend the ideas presented here to computing with functions
on the disk.
