The use of Gaussian processes (GPs) is supported by efficient sampling algorithms, a rich methodological literature, and strong theoretical grounding. However, due to their prohibitive computation and storage demands, the use of exact GPs in Bayesian models is limited to problems containing at most several thousand observations. Sampling requires matrix operations that scale at $\mathcal{O}(n^3),$ where $n$ is the number of unique inputs. Storage of individual matrices scales at $\mathcal{O}(n^2),$ and can quickly overwhelm the resources of most modern computers. To overcome these bottlenecks, we develop a sampling algorithm using $\mathcal{H}$ matrix approximation of the matrices comprising the GP posterior covariance. These matrices can approximate the true conditional covariance matrix within machine precision and allow for sampling algorithms that scale at $\mathcal{O}(n \ \mbox{log}^2 n)$ time and storage demands scaling at $\mathcal{O}(n \ \mbox{log} \ n).$ We also describe how these algorithms can be used as building blocks to model higher dimensional surfaces at $\mathcal{O}(d \ n \ \mbox{log}^2 n)$, where $d$ is the dimension of the surface under consideration, using tensor products of one-dimensional GPs. Though various scalable processes have been proposed for approximating Bayesian GP inference when $n$ is large, to our knowledge, none of these methods show that the approximation's Kullback-Leibler divergence to the true posterior can be made arbitrarily small and may be no worse than the approximation provided by finite computer arithmetic. We describe $\mathcal{H}-$matrices, give an efficient Gibbs sampler using these matrices for one-dimensional GPs, offer a proposed extension to higher dimensional surfaces, and investigate the performance of this fast increased fidelity approximate GP, FIFA-GP, using both simulated and real data sets.