The example you cited appears to be generating random Householder vectors and multiplying them out using backwards accumulation.
Another simple thing to do would be to generate a random matrix $\mathbf A$, then compute its $\mathbf A=\mathbf Q \mathbf R$ decomposition and discard the $\mathbf R$ factor. The two LAPACK functions that you need are [geqrf] (to factor $\mathbf A = \mathbf Q \mathbf R$ using Householder reflections) and [orgqr] (to reconstruct $\mathbf Q$ from the implicit reflector form). Although this is more work (about 2x) than the numpy algorithm, it might be faster in practice for large N because it will use BLAS3 kernels (whereas one-by-one backwards accumulation is only BLAS2).
There are a lot of C++ algebraic libraries, just look for one that happens to wrap these two functions. I happen to author/maintain such a library that has them (myramath), see below for a test program to generate a random orthogonal matrix:
myra::Matrix<double> rvs(int N)
auto A = myra::Matrix<double>::random(N,N);
auto tau = myra::geqrf_inplace(A);
// Form Q.
int N = 10;
auto Q = rvs(N);
std::cout << "Q = " << Q << std::endl;
// Check Q is orthogonal.
auto I = myra::Matrix<double>::identity(N);
std::cout << "|Q'Q-I| = " << myra::frobenius(myra::gemm(Q,'T',Q)-I) << std::endl;
std::cout << "|QQ'-I| = " << myra::frobenius(myra::gemm(Q,Q,'T')-I) << std::endl;