Computing the Pfaffian of a skew-symmetric matrix is a problem that arises in various fields of physics. Both computing the Pfaffian and a related problem, computing the canonical form of a skew-symmetric matrix under unitary congruence, can be solved easily once the skew-symmetric matrix has been reduced to skew-symmetric tridiagonal form. We develop efficient numerical methods for computing this tridiagonal form based on Gauss transformations, using a skew-symmetric, blocked form of the Parlett-Reid algorithm, or based on unitary transformations, using block Householder transformations and Givens rotations, that are applicable to dense and banded matrices, respectively. We also give a complete and fully optimized implementation of these algorithms in Python.