# include <vector>
# include <cstddef>

namespace CppAD { // BEGIN CppAD namespace
template <class Scalar>
Scalar det_of_minor(
    const std::vector<Scalar>& a  ,
    size_t                     m  ,
    size_t                     n  ,
    std::vector<size_t>&       r  ,
    std::vector<size_t>&       c  )
{
    const size_t R0 = r[m]; // R(0)
    size_t       Cj = c[m]; // C(j)    (case j = 0)
    size_t       Cj1 = m;   // C(j-1)  (case j = 0)

    // check for 1 by 1 case
    if( n == 1 ) return a[ R0 * m + Cj ];

    // initialize determinant of the minor M
    Scalar detM = Scalar(0);

    // initialize sign of factor for next sub-minor
    int s = 1;

    // remove row with index 0 in M from all the sub-minors of M
    r[m] = r[R0];

    // for each column of M
    for(size_t j = 0; j < n; j++)
    {   // element with index (0,j) in the minor M
        Scalar M0j = a[ R0 * m + Cj ];

        // remove column with index j in M to form next sub-minor S of M
        c[Cj1] = c[Cj];

        // compute determinant of the current sub-minor S
        Scalar detS = det_of_minor(a, m, n - 1, r, c);

        // restore column Cj to represenation of M as a minor of A
        c[Cj1] = Cj;

        // include this sub-minor term in the summation
        if( s > 0 )
            detM = detM + M0j * detS;
        else
            detM = detM - M0j * detS;

        // advance to next column of M
        Cj1 = Cj;
        Cj  = c[Cj];
        s   = - s;
    }

    // restore row zero to the minor representation for M
    r[m] = R0;

    // return the determinant of the minor M
    return detM;
}
} // END CppAD namespace