rmn_function(N,n,p)
{
        l_length(p)
        m_rbinom(N,n,p[1])
        if (l == 2) 
                { cbind(m,-m+n) } 
        else 
                { cbind(m,rmn(N,-m+n,p[2:l]/sum(p[2:l]))) }
}