  This paper presents an algorithmic method for generating random orthogonal
matrices \(A\) that satisfy the property \(A^t S A = S\), where \(S\) is a
fixed real invertible symmetric or skew-symmetric matrix. This method is
significant as it generalizes the procedures for generating orthogonal matrices
that fix a general fixed symmetric or skew-symmetric bilinear form. These
include orthogonal matrices that fall to groups such as the symplectic group,
Lorentz group, Poincar\'e group, and more generally the indefinite orthogonal
group, to name a few. These classes of matrices play crucial roles in diverse
fields such as theoretical physics, where they are used to describe symmetries
and conservation laws, as well as in computational geometry, numerical
analysis, and number theory, where they are integral to the study of quadratic
forms and modular forms. The implementation of our algorithms can be
accomplished using standard linear algebra libraries.
