  We adapt the Faddeev-LeVerrier algorithm for the computation of
characteristic polynomials to the computation of the Pfaffian of a
skew-symmetric matrix. This yields a very simple, easy to implement and
parallelize algorithm of computational cost $O(n^{\beta+1})$ where $n$ is the
size of the matrix and $O(n^{\beta})$ is the cost of multiplying $n\times
n$-matrices, $\beta\in[2,2.37286)$. We compare its performance to that of other
algorithms and show how it can be used to compute the Euler form of a
Riemannian manifold using computer algebra.
