This paper presents a novel recursive formulation for the simulation of multibody system dynamics based on Hamilton's canonical equations. Although Hamilton's canonical equations exhibit many advantageous features compared to their acceleration-based counterparts, it appears that there is a lack of dedicated parallel algorithms for multirigid body system dynamics based on this formulation. Serial kinematic chains are considered in this paper. Initially, the standard set of Hamilton's canonical equations is joined together with constraint equations at the velocity level. The reformulation determines the system's joint velocities and constraint force impulses in a divide and conquer manner resulting in logarithmic numerical cost for parallel implementation. Subsequently, the equations of motion are rearranged in order to obtain the time derivatives of the total joint momenta. In the case of the sequential implementation, the entire algorithm exhibits linear computational cost. Presented numerical method is exact and noniterative. Numerical test cases reveal negligible energy drift without the use of any additional constraint stabilization techniques. The results are compared against more standard acceleration based formulation and the outcome from real-life physical experiment.