As Fernando Muro points out in the comments, Sylvester's law of inertia is probably the easiest way to determine the signature. You diagonalize the symmetric matrix by the Gram-Schmidt process. This is essentially as easy as Gaussian elimination.
A symmetric matrix is equivalent to a bilinear form, with basis $\{v_1, \ldots, v_n\}$, and matrix entries inner products $\langle v_i,v_j\rangle, i,j\leq n$. If $\langle v_1,v_1\rangle \neq 0$, then replace $v_i$ with $v_i-v_1 \langle v_1,v_i\rangle/\langle v_1,v_1\rangle$, and proceed by induction (this has the effect of performing the same row and column operations, so as to keep the matrix symmetric). The only variation on the standard Gram-Schmidt is what to do when $\langle v_1,v_1\rangle=0$, which clearly means the form is indefinite. If $\langle v_1, v_i\rangle=0$ for all $i\leq n$, then you can proceed by induction. Otherwise, by permuting indices, you may assume that $\langle v_1,v_2\rangle \neq 0$. Replace $v_1$ with $v_2+tv_1$, such that $\langle v_2,v_2\rangle+2t\langle v_1,v_2\rangle\neq 0$, and proceed as with the usual Gram-Schmidt. Once you have diagonalized in this way, Sylvester's law implies that the signature of the diagonal matrix is the same as that of the original one.
Note added: Since you claim you're interested in computing the signatures of 4-manifolds, there's some special points one might be able to take advantage of. Certainly then the form will be integral unimodular, so presumably your symmetric matrix has integer entries and determinant $\pm 1$. If the 4-manifolds are smooth, then by Donaldson's theorem, if the intersection form is (positive) definite, then it is integrally equivalent to the diagonal form $1\oplus 1\oplus \cdots\oplus 1=1^n$. Otherwise, if the form is indefinite, then it is either odd and integrally equivalent to a diagonal form $1^k\oplus (-1)^{n-k}$, or it is even and there is a unique equivalence class for a given $n$ and signature. One can then perform the diagonalization process over $\mathbb{Z}$ using some version of the Euclidean algorithm, except in the even case (in which case there is a standard tridiagonal integral form). Actually, in the even case one could take a sum $\oplus \pm 1$ to make the form odd, and then diagonalize over the integers (a similar trick taking a sum with a hyperbolic form makes things integrally diagonalizable in the general case if you are also considering non-smooth manifolds which might have non-diagonalizable forms). See Milnor-Husemoller or Conway-Sloane for more details.
To diagonalize an odd indefinite symmetric bilinear form $B$ over $\mathbb{Z}$, first solve the equation $B(v_1,v_1)=\pm 1$ (this is where you can use some kind of Euclidean algorithm), then do the orthogonalization process to get $v_2,\ldots, v_n$ such that $B(v_1,v_i)=0, i\geq 2$, and induct.