This paper presents a high precision integration method for the dynamic response analysis of structures with holonomic constraints. A detail recursive scheme suitable for algebraic and differential equations (ADEs) which incorporates generalized forces is established. The matrix exponential involved in the scheme is calculated precisely using $2^N$ algorithm. The Taylor expansions of the nonlinear term concerned with state variables of the structure and the generalized constraint forces of the ADEs are derived and consequently, their particular integrals are obtained. The accuracy and effectiveness of the present method is demonstrated by two numerical examples, a plane truss with circular slot at its tip point and a slewing flexible cantilever beam which is currently interesting in optimal control of robot manipulators.