Hard-magnetic soft (HMS) beams made of soft polymer matrix embedded with hard-magnetic particles can generate large and fast deformation under magnetic stimulation. Dynamic modeling and simulation of HMS beams interacting with complex environment are challenging in terms of computational accuracy and efficiency. This paper presents a method for high-order modeling and efficient computation of HMS beams. The major contribution of the method is a new three-node HMS beam element of absolute nodal coordinate formulation (ANCF), which applies to two material models of nonlinear and linear elasticities (i.e. neoHookean and St. Venant-Kirchhoff) coupled with magnetic energy. To improve the efficiency of the method, the paper presents how to derive the generalized internal forces and their Jacobians via invariant tensors, and how to determine the generalized external forces to model dynamic loads and interactions including gravity, hydrodynamics in fluids, and frictional contact in pipelines. Afterwards, the paper gives both static and dynamic equations with Rayleigh damping and discusses the numerical algorithms. Finally, the paper makes a comparison of static analysis and the experimental observation to validate the accuracy of the proposed modeling method. The paper also discusses the dynamic simulations, including forced vibration, swimming motion, crawling locomotion, and navigating motion to demonstrate the predictive capability and efficacy of the proposed method for dynamic problems.