This paper presents the mathematical framework for the efficient construction of the equations of motion of complex multibody problems when uncertainty exists in the systems' parameters and/or inputs. Herein, uncertainty is prorogated through the system dynamics by using the method of polynomial chaos expansion (PCE). In this scheme, states of the system are projected onto the space of appropriate orthogonal base functions. Furthermore, the method of Divide-and-Conquer Algorithm (DCA) is extended to construct the equations of motion of the resulting non-deterministic system. In this scheme, the mathematical formulation to generate the projected handle equations of motion of the bodies and the projected constraint equations of the connecting joints are constructed in terms of the PCE. Finally, these projected equations are used to perform the assembly and disassembly passes to form and solve the equations of motion. The proposed method is highly parallelizable and scales down the computational complexity as a linear and logarithmic function of the state variables in serial and parallel implementations, respectively.