1. I looked up forward LU error bounds in Higham's Accuracy and Stability of Numerical Algorithms, Theorem 9.15 (citing Barrlund and Sun), for the LU decomposition
$$A=LU,\quad A+\Delta A=(L+\Delta L)(U+\Delta U)$$
it gives the norm-wise error bound
$$\begin{gather}
\frac{\|\Delta U\|_F}{\|U\|_2} \leq \frac{\|L^{-1}\|_2\|U^{-1}\|_2\|A\|_2}{1-\|L^{-1}\|_2\|U^{-1}\|_2\|\Delta A\|_2}\frac{\|\Delta A\|_F}{\|A\|_F} =: M, \\
\|\Delta A\|_\infty \leq n^2\gamma_{3n}\rho_n \|A\|_\infty,
\end{gather}$$
where $\gamma_k = k\epsilon/(1-k\epsilon)$ ($\epsilon$ is the unit roundoff), and $\rho_n$ is the growth factor, defined as
$$\max_{i,j,k} |a^{(k)}_{i,j}|/\max_{i,j} |a_{i,j}|$$
with $a^{(k)}_{i,j}$ being the matrix elements at $k$-th stage of LU factorization by Gaussian elimination.
So, in principle, so long as you can estimate all of the above numbers, the procedure should produce the correct determinant sign so long as the smallest value on the diagonal of $U$ is not so close to zero that its sign could have been influenced by numerical errors:
$$ \min_k |u_{kk}| \geq M\|U\|_2 \geq \|\Delta U\|_F. $$
Of course, this is not very elegant. Higham also gives a component-wise error bound, which should be stricter.
2. Not likely, but if you know something about how the numbers $a,b$ appear in off-diagonal elements, it is possible there might be a lower bound $m = \min_A |\det A|$. If it so happens that it is easy to calculate, then the determinant's sign can be accurately determined so long as the numerical errors are at most $m$.
3. Eigendecomposition, although more expensive, might also be helpful. For example, GSL (https://www.gnu.org/software/gsl/manual/html_node/Real-Symmetric-Matrices.html#Real-Symmetric-Matrices) promises that
The computed eigenvalues are accurate to an absolute accuracy of \epsilon ||A||_2, where \epsilon is the machine precision.
This would directly address the issue of whether the determinant's sign is correct - that would require all eigenvalues to satisfy $|\lambda| > \epsilon \|A\|_2$.
4. In general, computation of a sign of some quantity is ill-conditioned when that quantity is small - the determinant's sign would be very sensitive to small changes in the input matrix. So the common solution to this problem is just to look for a way to avoid computing the determinant's sign at all, but I don't know how feasible that is.
5. Arbitrary-precision floating-point arithmetic might also help. Although not with GSL, there are libraries (e.g., Eigen) that implement linear algebra in a way that can work with, for example, mpfr.
6. (Edit.) The $LDL^\top$ decomposition is even easier I think: from Barrlund (Eq. 2.1b), ($D$ is the diagonal from $U$, so I don't think it matters if you actually compute $LU$ instead because this is a perturbation analysis)
$$ \|\Delta D\|_F \leq \frac{\kappa_2(A)}{1-\kappa_2(A)\frac{\|\Delta A\|_2}{\|A\|_2}} \|\Delta A\|_F = \beta, $$
which now only needs estimates of $A$'s 2-norm condition number $\kappa_2(A)$ and the normwise error backward error estimates $\|\Delta A\|_2$, $\|\Delta A\|_F$ from above. The correct sign is now guaranteed by something like $|d_{kk}|\geq \beta$.
