We present a new framework for statistical inference on Riemannian manifolds that achieves high-order accuracy, addressing the challenges posed by non-Euclidean parameter spaces frequently encountered in modern data science. Our approach leverages a novel and computationally efficient procedure to reach higher-order asymptotic precision. In particular, we develop a bootstrap algorithm on Riemannian manifolds that is both computationally efficient and accurate for hypothesis testing and confidence region construction. Although locational hypothesis testing can be reformulated as a standard Euclidean problem, constructing high-order accurate confidence regions necessitates careful treatment of manifold geometry. To this end, we establish high-order asymptotics under an appropriate coordinate representation induced by a second-order retraction, thereby enabling precise expansions that incorporate curvature effects. We demonstrate the versatility of this framework across various manifold settings, including spheres, the Stiefel manifold, fixed-rank matrix manifolds, and rank-one tensor manifolds; for Euclidean submanifolds, we also introduce a class of projection-like coordinate charts with strong consistency properties. Finally, numerical studies confirm the practical merits of the proposed procedure.