This paper presents a new highly efficient procedure for the determination of the dynamical equations of motion for complex multibody systems and their subsequent temporal integration using parallel computing. The method is applicable to general systems of rigid bodies which may contain arbitrary joint types, multiple branches, and/or close loops. The method is based on the explicit determination of constraint forces at key joint locations and the subsequent highly efficient determination of system state time derivatives. The algorithm uses a novel hybrid direct and iterative solution scheme which allows a substantially higher degree of parallelization than is generally obtainable using the more conventional recursive O(N) procedures. It is shown that at the coarsest level the parallelization obtainable easily exceeds that indicated by the topology of the system. The procedure can produce a theoretical and time optimal O(log2N) performance on computational throughput with a processor optimal O(N) processors on a MDMD distributed architecture processing system.