The Gamma distribution belongs to the exponential families. Therefore, the Fisher information metric is $I(\theta)=\nabla^2 F(\theta)$. However, integrating the square root of the information matrix is difficult (no closed form solution). The author proceeds by characterizing the Riemannian geodesic using the differential equation relying on Christoffel symbols. Geodesics on the Gamma manifold are unique since the manifold is simply connected, complete with all sectional curvatures nonpositive. The authors come up with a Newton-like numerical optimization algorithm that depends on a good initialization. First, they show that the metric is bounded by Poincare metrics for which closed form equations of the geodesics are known. This yields a good starting tangent vector.
It is quite impressive to look at the formula of the closed-form equation of the Poincare geodesics. Those formula are surprisingly quite complicated.
The authors implemented their algorithm in FORTRAN and show that the algorithm always convergence on the domain examples, with high numerical precisions.