In this investigation, a new singularity-free formulation of a three-dimensional Euler–Bernoulli beam with large deformations and large rotations is developed. The position of the centroid line of the beam is integrated from its slope, which can be easily expressed by Euler parameters. The hyperspherical interpolation function is used to guarantee that the normalization constraint equation of Euler parameters is always satisfied. Each node of a beam element has only four nodal coordinates, which are significantly fewer than those in an absolute node coordinate formulation (ANCF) and the finite element method (FEM). Governing equations of the beam and constraint equations are derived using Lagrange's equations for systems with constraints, which are solved by a differential-algebraic equation (DAE) solver. The current formulation can be used to calculate the static equilibrium and linear and nonlinear dynamics of an Euler–Bernoulli beam under arbitrary, concentrated, and distributed forces. While the mass matrix is more complex than that in the ANCF, the stiffness matrix and generalized forces are simpler, which is amenable for calculating the equilibrium of the beam. Several numerical examples are presented to demonstrate the performance of the current formulation. It is shown that the current formulation can achieve the same accuracy as the ANCF and FEM with a fewer number of coordinates.