Dedicated to Doron Zeilberger on the occasion of his 60th birthday.Abstract. With the help of computer algebra we study the diagonal matrix elements Or p , where O = {1, β, iαnβ} are the standard Dirac matrix operators and the angular brackets denote the quantum-mechanical average for the relativistic Coulomb problem. Using Zeilberger's extension of Gosper's algorithm and a variant to it, three-term recurrence relations for each of these expectation values are derived together with some transformation formulas for the corresponding generalized hypergeometric series. In addition, the virial recurrence relations for these integrals are also found and proved algorithmically.