∑MathAdvanced

Geodesics & Exponential Map

Key Points

•
Geodesics are the “straightest possible” paths on curved spaces (manifolds) and locally minimize distance.
•
The exponential map $ $exp_{p}$ $ shoots out from a point $p$ along a geodesic in the direction of a tangent vector, landing at time 1.
•
On a sphere, geodesics are great circles; $ $exp_{p}$ (v) = $cos$ \ $∣ v ∥$ \,p + $sin$ \ $∣ v ∥$ \, $\frac{v}{∥ v ∥}$ $ provides a closed-form formula.
•
Numerically, geodesics satisfy a second-order ODE using Christoffel symbols: $ $\overset{x}{¨}^{k}$ + $Γ_{ij}^{k}$ $\overset{x}{˙}^{i}$ $\overset{x}{˙}^{j}$ = 0$.
•
The exponential and logarithm maps provide local coordinate systems (normal coordinates) and are key to optimization on manifolds.
•
Computationally, closed forms are O(1), while generic geodesic integration with RK methods is O(N) in the number of steps.
•
Beware coordinate singularities and step-size instability when integrating the geodesic ODE; robust normalization and small-angle handling are essential.
•
Applications range from robotics and computer graphics to statistics on manifolds and geometric machine learning.

Prerequisites

→Multivariable calculus — Gradients, Jacobians, and Taylor expansions are needed to understand metrics, derivatives, and local coordinates.
→Linear algebra — Inner products, matrices, and positive-definiteness underpin the metric tensor and Christoffel symbols.
→Ordinary differential equations — Geodesics satisfy a nonlinear second-order ODE; numerical integration methods are essential.
→Vector calculus in 3D — For embedded surfaces, dot products and partial derivatives define the induced metric.
→Basic differential geometry — Concepts like manifolds, tangent spaces, and connections provide the formal framework for geodesics and exp/log.

Detailed Explanation

Tap terms for definitions

01Overview

Imagine you are walking on the surface of the Earth. If you keep pointing in the same compass direction without turning, you trace a great circle—this is a geodesic. In curved spaces (manifolds), geodesics generalize straight lines from flat Euclidean geometry: they are locally shortest paths and the straightest possible trajectories constrained to the surface. The exponential map formalizes the idea of “shooting” from a point along a given direction (a tangent vector) to reach another point on the manifold. Given a point p and a tangent vector v at p, the exponential map exp_p(v) returns the point reached by following the unique geodesic starting at p with initial velocity v for unit time. Together, geodesics and the exponential map let us measure distances, transport information, and build coordinates without flattening the surface. Practically, we use geodesics to compute shortest routes on curved surfaces, to interpolate rotations and shapes, and to optimize functions where variables live on manifolds (like unit spheres or spaces of rotations). On some manifolds (e.g., the sphere or rotation group SO(3)), exp and log have elegant closed forms that are numerically stable with care. On general manifolds, we integrate the geodesic differential equations using the metric and Christoffel symbols. These tools bridge intuition (“walk straight ahead”) and computation (ODE solvers), enabling algorithms in graphics, robotics, and machine learning that respect intrinsic geometry.

02Intuition & Analogies

Think of driving on a perfectly smooth, hilly landscape without ever touching the steering wheel. Your car will follow the road that feels straight based only on the surface beneath you—no turning left or right. That path, which keeps your heading constant relative to the surface, is a geodesic. In flat parking lots, these are just straight lines. On a globe, they’re great circles like the equator or any meridian. Now imagine a launcher at a point p on the surface. You set a direction (a tangent vector v) and a speed (its length). If the launcher fires for exactly one second, where does the projectile land if it must stay glued to the surface and always keep its wheel alignment straight? That landing point is the exponential map exp_p(v). If you halve the firing time, you get exp_p(0.5 v), which lies halfway along the geodesic. If you reverse the direction, you run the geodesic backward: exp_p(-v). Close to p, distances on the surface look almost Euclidean: walking a small vector v from p along the geodesic travels roughly length ||v||. This is why exp_p provides excellent local coordinates (normal coordinates). But as you go farther, curvature matters: paths bend, and going “straight” may eventually circle back (as on a sphere). The exponential map captures both the local straightness and the global bending imposed by curvature.

03Formal Definition

Let (M, g) be a smooth Riemannian manifold with Levi–Civita connection \(

\nabla

\). A smooth curve \(

γ

: I

\to

M\) is a geodesic if its covariant acceleration vanishes: \(

\nabla_{\overset{γ}{˙}} \overset{γ}{˙}

= 0\). In local coordinates \((

x^{1}

\dots

x^{n}

)\), this is the geodesic equation \(

\overset{x}{¨}^{k}

Γ_{ij}^{k}

(x)\,

\overset{x}{˙}^{i}

\overset{x}{˙}^{j}

= 0\), where \(

Γ_{ij}^{k}

\) are Christoffel symbols of g. Geodesics are critical points of the energy functional \(E(

γ

) =

\frac{1}{2} \int

g_{ij}

(

γ

)\,

\overset{x}{˙}^{i}

\overset{x}{˙}^{j}

\, dt\). For sufficiently short intervals, they also minimize the length functional. For any point \(p

\in

M\) and tangent vector \(v

\in

T_{p} M

\), there exists a unique geodesic \(

γ_{v}

\) with \(

γ_{v}

(0)=p\) and \(

\overset{γ}{˙}_{v}

(0)=v\). The exponential map at p is defined by \(

exp_{p}

(v) =

γ_{v}

(1)\). On a neighborhood of the origin in \(

T_{p} M

\), \(

exp_{p}

\) is a diffeomorphism onto a normal neighborhood of p, and \(d(

exp_{p}

)_0 =

Id_{T_{p} M}

\). The injectivity radius at p is the largest radius for which \(

exp_{p}

\) is a diffeomorphism on the ball of that radius in \(

T_{p} M

\). The cut locus is where geodesics cease to minimize distance or lose uniqueness.

04When to Use

Shortest paths on surfaces: computing routes on terrains, globe navigation, or inside curved spaces where straight lines are not available.
Optimization on manifolds: gradient descent on spheres, Stiefel/Grassmann manifolds, or rotation groups uses exp maps for retractions and log maps for measuring errors.
Interpolation and averaging: spherical linear interpolation (slerp) and Frechét/Karcher means on manifolds need exp/log to move between points and tangent spaces.
Computer graphics and geometry processing: parameterizing meshes, computing geodesic distances, and tracing intrinsic curves on surfaces.
Robotics and control: state spaces like SO(3)×R^3 use geodesics to plan motions respecting configuration space geometry.
Statistics on manifolds: diffusion, kernel methods, and distributions (e.g., von Mises–Fisher on the sphere) rely on geodesic distances and exponential coordinates. Use closed-form exp/log when available (sphere, SO(3), SPD matrices) for stability and speed. Use numerical ODE integration (e.g., RK4) with Christoffel symbols for general embedded surfaces or metrics. Prefer normal coordinates via exp for local linear approximations, but beware of injectivity radius when moving far from the base point.

⚠️Common Mistakes

Confusing “shortest” with “straightest” globally: geodesics are locally minimizing, not always globally shortest (e.g., antipodal points on a sphere have many geodesics, and great circles stop minimizing past the cut locus).
Ignoring coordinate singularities: latitude–longitude charts have singularities at the poles; integrating geodesics there can explode unless you change charts or use embeddings/orthonormal frames.
Neglecting normalization in closed forms: on the sphere, failing to normalize points to unit length or to handle small angles leads to large numerical errors (e.g., dividing by (\sin \theta) when (\theta \approx 0)).
Large time steps in ODE solvers: geodesic equations are nonlinear; too large a step in RK methods can drift off the surface or accumulate energy error; project back to the surface and/or adapt step sizes.
Misusing exp/log beyond injectivity radius: the log map becomes multivalued past the cut locus (e.g., antipodes on S^2); optimization steps should be kept within regions where exp/log are well-defined and unique.
Mixing ambient and intrinsic quantities: computing with Euclidean distances between embedded points instead of geodesic distances on the manifold leads to incorrect results. Always use the metric-induced operations.

Key Formulas

Curve length

L (γ) = \int_{a}^{b} g_{ij} (γ (t)) \overset{x}{˙}^{i} (t) \overset{x}{˙}^{j} (t) d t

Explanation: The geodesic length is computed using the metric along the curve. Geodesics locally minimize this length among nearby curves with the same endpoints.

Energy functional

E (γ) = \frac{1}{2} \int_{a}^{b} g_{ij} (γ (t)) \overset{x}{˙}^{i} (t) \overset{x}{˙}^{j} (t) d t

Explanation: Geodesics are stationary points of the energy. For fixed endpoints, minimizing energy is equivalent to minimizing length among constant-speed curves.

Geodesic equation

\overset{x}{¨}^{k} + Γ_{ij}^{k} (x) \overset{x}{˙}^{i} \overset{x}{˙}^{j} = 0

Explanation: This second-order ODE expresses zero covariant acceleration. Given initial position and velocity, it determines a unique geodesic.

Christoffel symbols

Γ_{ij}^{k} = \frac{1}{2} g^{k ℓ} (\partial_{i} g_{j ℓ} + \partial_{j} g_{i ℓ} - \partial_{ℓ} g_{ij})

Explanation: Christoffel symbols are computed from the metric and its first derivatives. They encode how coordinates curve and are used in geodesic ODEs.

Exponential map definition

exp_{p} (v) = γ_{v} (1), where \nabla_{\overset{γ}{˙}_{v}} \overset{γ}{˙}_{v} = 0, γ_{v} (0) = p, \overset{γ}{˙}_{v} (0) = v

Explanation: Starting at p with initial velocity v, follow the geodesic for unit time. The endpoint is the exponential map value.

Differential at zero

d (exp_{p})_{0} = Id_{T_{p} M}

Explanation: Infinitesimally, the exponential map behaves like the identity from the tangent space to the manifold. This underlies normal coordinates.

Local distance via exp

dist (p, exp_{p} (v)) = ∥ v ∥ + O (∥ v ∥^{3})

Explanation: For small v, geodesic distance from p to ex $p_{p}$ (v) equals the tangent-space norm up to cubic error. This justifies linear approximations locally.

Sphere exponential

exp_{p} (v) = cos ∥ v ∥ p + sin ∥ v ∥ \frac{v}{∥ v ∥}, v ⊥ p, p \in S^{2}

Explanation: On the unit sphere, shooting along the great circle has a closed form. It’s fast and stable if small-angle cases are handled carefully.

Sphere logarithm

lo g_{p} (q) = θ \frac{q - ( p \cdot q ) p}{∥ q - ( p \cdot q ) p ∥}, θ = arccos (p \cdot q), p, q \in S^{2}, q \neq = - p

Explanation: For points not antipodal, the log map returns the tangent vector that reaches q along the shortest great-circle arc from p.

Parallel transport

\nabla_{\overset{γ}{˙}} V = 0

Explanation: A vector field V along γ is parallel if its covariant derivative vanishes. Parallel transport preserves lengths and angles along geodesics.

Jacobi equation

\frac{D ^{2} J}{d t ^{2}} + R (J, \overset{γ}{˙}) \overset{γ}{˙} = 0

Explanation: Jacobi fields describe variations of geodesics and detect conjugate points. Curvature tensor R couples nearby geodesics.

Complexity Analysis

Closed-form geodesic computations on special manifolds (like

S^{2}

or SO(3)) run in constant time and space: O(1) arithmetic and a handful of trigonometric evaluations. These are numerically stable when inputs are normalized and small-angle cases use series expansions to avoid division by near-zero quantities. The memory footprint is O(1), holding a few vectors and scalars. For general manifolds specified by a parametric embedding X(u) or a metric g(u), computing geodesics requires numerically integrating the second-order geodesic ODE. A standard approach reduces it to a first-order system in positions and velocities and integrates with an explicit method like RK4. If N steps are used, and each step requires evaluating the metric, its derivatives, and inverting a small matrix (e.g., 2×2 for surfaces), the time complexity is O(N). The constant factor includes several evaluations of X and its partial derivatives (often via finite differences), building g and ∂g, assembling Christoffel symbols (O(

d^{3}

) with tiny d), and performing updates. Space complexity remains O(1) beyond storing the current state, although saving a trajectory of T samples costs O(T). Accuracy and stability depend on step size h: local truncation error is O(

h^{5}

) for RK4, yielding global error O(

h^{4}

). Curvature and coordinate singularities can amplify errors, possibly requiring adaptive step sizes or reparameterization. Projection back to the manifold (e.g., renormalizing to the sphere) can control drift when integrating in an embedding, at the cost of bias. When nearest-neighbor or mesh-based approximations (e.g., graph shortest paths) are used instead of intrinsic ODEs, complexity depends on graph size (Dijkstra O(E log V)), but these approximate true geodesics only coarsely unless the mesh is dense.

Code Examples

Closed-form exponential and logarithm maps on the 2-sphere (S^2) with geodesic evaluation

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Vec3 {
5     double x, y, z;
6     Vec3(double x_=0, double y_=0, double z_=0): x(x_), y(y_), z(z_) {}
7     Vec3 operator+(const Vec3& o) const { return Vec3(x+o.x, y+o.y, z+o.z); }
8     Vec3 operator-(const Vec3& o) const { return Vec3(x-o.x, y-o.y, z-o.z); }
9     Vec3 operator*(double s) const { return Vec3(x*s, y*s, z*s); }
10 };
11 
12 static inline double dot(const Vec3& a, const Vec3& b){ return a.x*b.x + a.y*b.y + a.z*b.z; }
13 static inline double norm(const Vec3& a){ return sqrt(max(0.0, dot(a,a))); }
14 static inline Vec3 normalize(const Vec3& a){ double n = norm(a); if(n==0) return a; return a*(1.0/n); }
15 
16 // Project a vector to the tangent space at p on S^2
17 static inline Vec3 tangent_project(const Vec3& p, const Vec3& v){
18     return v - p * dot(p,v); // remove radial component
19 }
20 
21 // Exponential map on unit sphere S^2 at point p for tangent vector v (v must be tangent at p)
22 Vec3 sphere_exp(const Vec3& p_unit, const Vec3& v_tan){
23     Vec3 p = normalize(p_unit);
24     Vec3 v = tangent_project(p, v_tan); // ensure tangency
25     double theta = norm(v);
26     if(theta < 1e-12) {
27         // exp_p(0) = p, and for very small theta use first-order approximation
28         return p + v; // practically the same; caller may renormalize
29     }
30     Vec3 vhat = v * (1.0/theta);
31     Vec3 q = p * cos(theta) + vhat * sin(theta);
32     return normalize(q); // stay on the sphere
33 }
34 
35 // Logarithm map on unit sphere S^2 at point p for point q (q != -p)
36 Vec3 sphere_log(const Vec3& p_unit, const Vec3& q_unit){
37     Vec3 p = normalize(p_unit);
38     Vec3 q = normalize(q_unit);
39     double c = std::clamp(dot(p,q), -1.0, 1.0);
40     double theta = acos(c);
41     if(theta < 1e-12){
42         return Vec3(0,0,0); // points coincide; zero tangent vector
43     }
44     if(fabs(M_PI - theta) < 1e-9){
45         // Antipodal: log is not unique; return one of many choices (undefined direction)
46         // Here we pick an arbitrary tangent vector orthogonal to p
47         Vec3 arbitrary = fabs(p.x) < 0.9 ? Vec3(1,0,0) : Vec3(0,1,0);
48         Vec3 u = normalize(tangent_project(p, arbitrary));
49         return u * theta; // magnitude pi, arbitrary direction
50     }
51     Vec3 u = q - p * c; // component of q orthogonal to p
52     double s = norm(u);
53     Vec3 uhat = u * (1.0 / s);
54     return uhat * theta;
55 }
56 
57 // Evaluate the geodesic starting at p with tangent v at parameter t: gamma(t) = exp_p(t v)
58 Vec3 sphere_geodesic_point(const Vec3& p, const Vec3& v_tan, double t){
59     return sphere_exp(p, v_tan * t);
60 }
61 
62 int main(){
63     // Example: start at p on the equator, shoot east with unit-speed (unit tangent length)
64     Vec3 p = normalize(Vec3(1,0,0));
65     Vec3 east = tangent_project(p, Vec3(0,1,0)); // already orthogonal to p
66     east = east * (1.0 / max(1e-12, norm(east))); // unit tangent
67 
68     // Move 90 degrees (pi/2) along the great circle
69     Vec3 q = sphere_geodesic_point(p, east, M_PI/2);
70 
71     cout.setf(ios::fixed); cout<<setprecision(6);
72     cout << "q = (" << q.x << ", " << q.y << ", " << q.z << ")\n";
73 
74     // Recover the tangent using log
75     Vec3 v = sphere_log(p, q);
76     cout << "||v|| (should be ~pi/2) = " << norm(v) << "\n";
77 
78     // Verify round-trip exp(log)
79     Vec3 q2 = sphere_exp(p, v);
80     cout << "round-trip dot(q,q2) = " << dot(q, q2) << " (1 means perfect)\n";
81     return 0;
82 }
83

This program implements closed-form exponential and logarithm maps on the unit sphere S^2 and uses them to evaluate points along great-circle geodesics. The exponential map rotates the base point p by angle ||v|| in the tangent direction v, while the logarithm map recovers the tangent vector that reaches a nearby point q along the shortest arc. The code carefully handles tangency, normalization, and small-angle/antipodal edge cases. Because S^2 admits analytic formulas, the computation is fast and stable.

Time: O(1) per exp/log evaluation (a few dot products and trig calls)Space: O(1)

Numerical geodesic integration on a parametric surface via Christoffel symbols (paraboloid)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Parametric surface: paraboloid X(u,v) = (u, v, 0.5*(u*u + v*v))
5 struct Surface {
6     static array<double,3> X(double u, double v){
7         return {u, v, 0.5*(u*u + v*v)};
8     }
9 };
10 
11 struct MetricData {
12     // metric g_ij and its partial derivatives dg_ij/du, dg_ij/dv at (u,v)
13     double g[2][2];
14     double du[2][2];
15     double dv[2][2];
16 };
17 
18 static inline array<double,3> operator+(const array<double,3>& a, const array<double,3>& b){ return {a[0]+b[0], a[1]+b[1], a[2]+b[2]}; }
19 static inline array<double,3> operator-(const array<double,3>& a, const array<double,3>& b){ return {a[0]-b[0], a[1]-b[1], a[2]-b[2]}; }
20 static inline array<double,3> operator*(const array<double,3>& a, double s){ return {a[0]*s, a[1]*s, a[2]*s}; }
21 static inline double dot3(const array<double,3>& a, const array<double,3>& b){ return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
22 
23 // Finite differences for partial derivatives of X
24 static inline array<double,3> Xu(double u, double v, double h){
25     auto Xm = Surface::X(u-h, v);
26     auto Xp = Surface::X(u+h, v);
27     return (Xp - Xm) * (1.0/(2*h));
28 }
29 static inline array<double,3> Xv(double u, double v, double h){
30     auto Xm = Surface::X(u, v-h);
31     auto Xp = Surface::X(u, v+h);
32     return (Xp - Xm) * (1.0/(2*h));
33 }
34 
35 // Second derivatives
36 static inline array<double,3> Xuu(double u, double v, double h){
37     auto Xm = Surface::X(u-h, v);
38     auto X0 = Surface::X(u, v);
39     auto Xp = Surface::X(u+h, v);
40     return (Xp - X0*2.0 + Xm) * (1.0/(h*h));
41 }
42 static inline array<double,3> Xvv(double u, double v, double h){
43     auto Xm = Surface::X(u, v-h);
44     auto X0 = Surface::X(u, v);
45     auto Xp = Surface::X(u, v+h);
46     return (Xp - X0*2.0 + Xm) * (1.0/(h*h));
47 }
48 static inline array<double,3> Xuv(double u, double v, double h){
49     auto Xmm = Surface::X(u-h, v-h);
50     auto Xmp = Surface::X(u-h, v+h);
51     auto Xpm = Surface::X(u+h, v-h);
52     auto Xpp = Surface::X(u+h, v+h);
53     return (Xpp - Xpm - Xmp + Xmm) * (1.0/(4*h*h));
54 }
55 
56 // Build metric tensor and its partial derivatives via Xu, Xv and their derivatives
57 static inline MetricData buildMetric(double u, double v){
58     const double h = 1e-4; // finite-difference step
59     auto Xu_ = Xu(u,v,h), Xv_ = Xv(u,v,h);
60     // First fundamental form g = [E F; F G]
61     double E = dot3(Xu_, Xu_);
62     double F = dot3(Xu_, Xv_);
63     double G = dot3(Xv_, Xv_);
64 
65     // Derivatives of Xu, Xv
66     auto Xuu_ = Xuu(u,v,h), Xuv_ = Xuv(u,v,h), Xvv_ = Xvv(u,v,h);
67 
68     // Partial derivatives of metric entries
69     // dE/du = 2 <Xu, Xuu>, dE/dv = 2 <Xu, Xuv>
70     double dE_du = 2.0 * dot3(Xu_, Xuu_);
71     double dE_dv = 2.0 * dot3(Xu_, Xuv_);
72     // dF/du = <Xuu, Xv> + <Xu, Xuv>, dF/dv = <Xuv, Xv> + <Xu, Xvv>
73     double dF_du = dot3(Xuu_, Xv_) + dot3(Xu_, Xuv_);
74     double dF_dv = dot3(Xuv_, Xv_) + dot3(Xu_, Xvv_);
75     // dG/du = 2 <Xv, Xuv>, dG/dv = 2 <Xv, Xvv>
76     double dG_du = 2.0 * dot3(Xv_, Xuv_);
77     double dG_dv = 2.0 * dot3(Xv_, Xvv_);
78 
79     MetricData md{};
80     md.g[0][0] = E; md.g[0][1] = F; md.g[1][0] = F; md.g[1][1] = G;
81     md.du[0][0] = dE_du; md.du[0][1] = dF_du; md.du[1][0] = dF_du; md.du[1][1] = dG_du;
82     md.dv[0][0] = dE_dv; md.dv[0][1] = dF_dv; md.dv[1][0] = dF_dv; md.dv[1][1] = dG_dv;
83     return md;
84 }
85 
86 // Compute inverse g^{ij} and Christoffel symbols Gamma^k_{ij} for k,i,j in {0,1}
87 static inline void christoffel(double u, double v, double Gamma[2][2][2]){
88     MetricData md = buildMetric(u,v);
89     double E = md.g[0][0], F = md.g[0][1], G = md.g[1][1];
90     double det = E*G - F*F;
91     // Inverse metric
92     double inv[2][2];
93     inv[0][0] =  G/det; inv[0][1] = -F/det;
94     inv[1][0] = -F/det; inv[1][1] =  E/det;
95 
96     // Helpers: partials of g_ij; index 0->u, 1->v
97     double dg[2][2][2];
98     dg[0][0][0] = md.du[0][0]; dg[1][0][0] = md.dv[0][0];
99     dg[0][0][1] = md.du[0][1]; dg[1][0][1] = md.dv[0][1];
100     dg[0][1][0] = md.du[1][0]; dg[1][1][0] = md.dv[1][0];
101     dg[0][1][1] = md.du[1][1]; dg[1][1][1] = md.dv[1][1];
102 
103     // Christoffel: Gamma^k_{ij} = 1/2 g^{kl}( d_i g_{jl} + d_j g_{il} - d_l g_{ij} )
104     for(int k=0;k<2;k++){
105         for(int i=0;i<2;i++){
106             for(int j=0;j<2;j++){
107                 double sum = 0.0;
108                 for(int l=0;l<2;l++){
109                     double term = dg[i][j][l] + dg[j][i][l] - dg[l][i][j];
110                     sum += inv[k][l] * 0.5 * term;
111                 }
112                 Gamma[k][i][j] = sum;
113             }
114         }
115     }
116 }
117 
118 // State y = [u, v, du, dv]; RHS dy/dt
119 static inline array<double,4> geodesic_rhs(double u, double v, double du, double dv){
120     double Gamma[2][2][2];
121     christoffel(u,v,Gamma);
122     // Second derivatives from geodesic ODE: d2u = -Gamma^u_{ij} dxi dxj, d2v = -Gamma^v_{ij} dxi dxj
123     double d2u = - (Gamma[0][0][0]*du*du + 2.0*Gamma[0][0][1]*du*dv + Gamma[0][1][1]*dv*dv);
124     double d2v = - (Gamma[1][0][0]*du*du + 2.0*Gamma[1][0][1]*du*dv + Gamma[1][1][1]*dv*dv);
125     return {du, dv, d2u, d2v};
126 }
127 
128 // One RK4 step for the 4D system
129 static inline void rk4_step(array<double,4>& y, double h){
130     auto f = [](const array<double,4>& s){ return geodesic_rhs(s[0], s[1], s[2], s[3]); };
131     array<double,4> k1 = f(y);
132     array<double,4> y2{ y[0]+0.5*h*k1[0], y[1]+0.5*h*k1[1], y[2]+0.5*h*k1[2], y[3]+0.5*h*k1[3] };
133     array<double,4> k2 = f(y2);
134     array<double,4> y3{ y[0]+0.5*h*k2[0], y[1]+0.5*h*k2[1], y[2]+0.5*h*k2[2], y[3]+0.5*h*k2[3] };
135     array<double,4> k3 = f(y3);
136     array<double,4> y4{ y[0]+h*k3[0], y[1]+h*k3[1], y[2]+h*k3[2], y[3]+h*k3[3] };
137     array<double,4> k4 = f(y4);
138     for(int i=0;i<4;i++){
139         y[i] += (h/6.0) * (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i]);
140     }
141 }
142 
143 int main(){
144     // Initial param position and tangent (choose small region where chart is regular)
145     double u0 = 0.0, v0 = 0.0; // apex of paraboloid
146     // Initial direction in parameters: choose unit-speed w.r.t. metric approximately
147     // We roughly normalize using the metric at start
148     MetricData md0 = buildMetric(u0, v0);
149     // Speed^2 = [du dv] g [du dv]^T. Choose raw (1, 0.2) then scale.
150     double du0 = 1.0, dv0 = 0.2;
151     double speed2 = md0.g[0][0]*du0*du0 + 2*md0.g[0][1]*du0*dv0 + md0.g[1][1]*dv0*dv0;
152     double scale = 1.0 / sqrt(max(1e-12, speed2));
153     du0 *= scale; dv0 *= scale; // roughly unit-speed
154 
155     array<double,4> y{u0, v0, du0, dv0};
156 
157     double h = 0.01; // step size
158     int steps = 1000;
159 
160     cout.setf(ios::fixed); cout<<setprecision(6);
161     for(int i=0;i<steps;i++){
162         // Integrate one step
163         rk4_step(y, h);
164         // Get 3D point on surface
165         auto X3 = Surface::X(y[0], y[1]);
166         if(i % 50 == 0){
167             cout << "t=" << (i+1)*h << ": X=(" << X3[0] << ", " << X3[1] << ", " << X3[2] << ")\n";
168         }
169     }
170 
171     return 0;
172 }
173

This program integrates geodesics on the paraboloid z = 0.5(x^2 + y^2), parameterized by (u,v) -> (u, v, 0.5(u^2+v^2)). It builds the induced metric from the embedding using finite-difference derivatives, computes Christoffel symbols, and advances the geodesic with an RK4 integrator. The state consists of parameters (u,v) and their velocities, which evolve according to the geodesic ODE. The code prints sampled 3D points along the intrinsic straightest path. Small step sizes and finite-difference epsilons improve accuracy; more sophisticated implementations would add adaptive stepping and error control.

Time: O(N) for N integration steps; each step performs a constant number of surface and derivative evaluationsSpace: O(1) for the integrator state; O(T) if storing T sampled points

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Vec3 {
5	double x, y, z;
6	Vec3(double x_=0, double y_=0, double z_=0): x(x_), y(y_), z(z_) {}
7	Vec3 operator+(const Vec3& o) const { return Vec3(x+o.x, y+o.y, z+o.z); }
8	Vec3 operator-(const Vec3& o) const { return Vec3(x-o.x, y-o.y, z-o.z); }
9	Vec3 operator(double s) const { return Vec3(xs, ys, zs); }
10	};
11
12	static inline double dot(const Vec3& a, const Vec3& b){ return a.xb.x + a.yb.y + a.z*b.z; }
13	static inline double norm(const Vec3& a){ return sqrt(max(0.0, dot(a,a))); }
14	static inline Vec3 normalize(const Vec3& a){ double n = norm(a); if(n==0) return a; return a*(1.0/n); }
15
16	// Project a vector to the tangent space at p on S^2
17	static inline Vec3 tangent_project(const Vec3& p, const Vec3& v){
18	return v - p * dot(p,v); // remove radial component
19	}
20
21	// Exponential map on unit sphere S^2 at point p for tangent vector v (v must be tangent at p)
22	Vec3 sphere_exp(const Vec3& p_unit, const Vec3& v_tan){
23	Vec3 p = normalize(p_unit);
24	Vec3 v = tangent_project(p, v_tan); // ensure tangency
25	double theta = norm(v);
26	if(theta < 1e-12) {
27	// exp_p(0) = p, and for very small theta use first-order approximation
28	return p + v; // practically the same; caller may renormalize
29	}
30	Vec3 vhat = v * (1.0/theta);
31	Vec3 q = p * cos(theta) + vhat * sin(theta);
32	return normalize(q); // stay on the sphere
33	}
34
35	// Logarithm map on unit sphere S^2 at point p for point q (q != -p)
36	Vec3 sphere_log(const Vec3& p_unit, const Vec3& q_unit){
37	Vec3 p = normalize(p_unit);
38	Vec3 q = normalize(q_unit);
39	double c = std::clamp(dot(p,q), -1.0, 1.0);
40	double theta = acos(c);
41	if(theta < 1e-12){
42	return Vec3(0,0,0); // points coincide; zero tangent vector
43	}
44	if(fabs(M_PI - theta) < 1e-9){
45	// Antipodal: log is not unique; return one of many choices (undefined direction)
46	// Here we pick an arbitrary tangent vector orthogonal to p
47	Vec3 arbitrary = fabs(p.x) < 0.9 ? Vec3(1,0,0) : Vec3(0,1,0);
48	Vec3 u = normalize(tangent_project(p, arbitrary));
49	return u * theta; // magnitude pi, arbitrary direction
50	}
51	Vec3 u = q - p * c; // component of q orthogonal to p
52	double s = norm(u);
53	Vec3 uhat = u * (1.0 / s);
54	return uhat * theta;
55	}
56
57	// Evaluate the geodesic starting at p with tangent v at parameter t: gamma(t) = exp_p(t v)
58	Vec3 sphere_geodesic_point(const Vec3& p, const Vec3& v_tan, double t){
59	return sphere_exp(p, v_tan * t);
60	}
61
62	int main(){
63	// Example: start at p on the equator, shoot east with unit-speed (unit tangent length)
64	Vec3 p = normalize(Vec3(1,0,0));
65	Vec3 east = tangent_project(p, Vec3(0,1,0)); // already orthogonal to p
66	east = east * (1.0 / max(1e-12, norm(east))); // unit tangent
67
68	// Move 90 degrees (pi/2) along the great circle
69	Vec3 q = sphere_geodesic_point(p, east, M_PI/2);
70
71	cout.setf(ios::fixed); cout<<setprecision(6);
72	cout << "q = (" << q.x << ", " << q.y << ", " << q.z << ")\n";
73
74	// Recover the tangent using log
75	Vec3 v = sphere_log(p, q);
76	cout << "\|\|v\|\| (should be ~pi/2) = " << norm(v) << "\n";
77
78	// Verify round-trip exp(log)
79	Vec3 q2 = sphere_exp(p, v);
80	cout << "round-trip dot(q,q2) = " << dot(q, q2) << " (1 means perfect)\n";
81	return 0;
82	}
83

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Parametric surface: paraboloid X(u,v) = (u, v, 0.5(uu + v*v))
5	struct Surface {
6	static array<double,3> X(double u, double v){
7	return {u, v, 0.5(uu + v*v)};
8	}
9	};
10
11	struct MetricData {
12	// metric g_ij and its partial derivatives dg_ij/du, dg_ij/dv at (u,v)
13	double g[2][2];
14	double du[2][2];
15	double dv[2][2];
16	};
17
18	static inline array<double,3> operator+(const array<double,3>& a, const array<double,3>& b){ return {a[0]+b[0], a[1]+b[1], a[2]+b[2]}; }
19	static inline array<double,3> operator-(const array<double,3>& a, const array<double,3>& b){ return {a[0]-b[0], a[1]-b[1], a[2]-b[2]}; }
20	static inline array<double,3> operator(const array<double,3>& a, double s){ return {a[0]s, a[1]s, a[2]s}; }
21	static inline double dot3(const array<double,3>& a, const array<double,3>& b){ return a[0]b[0]+a[1]b[1]+a[2]*b[2]; }
22
23	// Finite differences for partial derivatives of X
24	static inline array<double,3> Xu(double u, double v, double h){
25	auto Xm = Surface::X(u-h, v);
26	auto Xp = Surface::X(u+h, v);
27	return (Xp - Xm) * (1.0/(2*h));
28	}
29	static inline array<double,3> Xv(double u, double v, double h){
30	auto Xm = Surface::X(u, v-h);
31	auto Xp = Surface::X(u, v+h);
32	return (Xp - Xm) * (1.0/(2*h));
33	}
34
35	// Second derivatives
36	static inline array<double,3> Xuu(double u, double v, double h){
37	auto Xm = Surface::X(u-h, v);
38	auto X0 = Surface::X(u, v);
39	auto Xp = Surface::X(u+h, v);
40	return (Xp - X02.0 + Xm) (1.0/(h*h));
41	}
42	static inline array<double,3> Xvv(double u, double v, double h){
43	auto Xm = Surface::X(u, v-h);
44	auto X0 = Surface::X(u, v);
45	auto Xp = Surface::X(u, v+h);
46	return (Xp - X02.0 + Xm) (1.0/(h*h));
47	}
48	static inline array<double,3> Xuv(double u, double v, double h){
49	auto Xmm = Surface::X(u-h, v-h);
50	auto Xmp = Surface::X(u-h, v+h);
51	auto Xpm = Surface::X(u+h, v-h);
52	auto Xpp = Surface::X(u+h, v+h);
53	return (Xpp - Xpm - Xmp + Xmm) * (1.0/(4hh));
54	}
55
56	// Build metric tensor and its partial derivatives via Xu, Xv and their derivatives
57	static inline MetricData buildMetric(double u, double v){
58	const double h = 1e-4; // finite-difference step
59	auto Xu_ = Xu(u,v,h), Xv_ = Xv(u,v,h);
60	// First fundamental form g = [E F; F G]
61	double E = dot3(Xu_, Xu_);
62	double F = dot3(Xu_, Xv_);
63	double G = dot3(Xv_, Xv_);
64
65	// Derivatives of Xu, Xv
66	auto Xuu_ = Xuu(u,v,h), Xuv_ = Xuv(u,v,h), Xvv_ = Xvv(u,v,h);
67
68	// Partial derivatives of metric entries
69	// dE/du = 2 <Xu, Xuu>, dE/dv = 2 <Xu, Xuv>
70	double dE_du = 2.0 * dot3(Xu_, Xuu_);
71	double dE_dv = 2.0 * dot3(Xu_, Xuv_);
72	// dF/du = <Xuu, Xv> + <Xu, Xuv>, dF/dv = <Xuv, Xv> + <Xu, Xvv>
73	double dF_du = dot3(Xuu_, Xv_) + dot3(Xu_, Xuv_);
74	double dF_dv = dot3(Xuv_, Xv_) + dot3(Xu_, Xvv_);
75	// dG/du = 2 <Xv, Xuv>, dG/dv = 2 <Xv, Xvv>
76	double dG_du = 2.0 * dot3(Xv_, Xuv_);
77	double dG_dv = 2.0 * dot3(Xv_, Xvv_);
78
79	MetricData md{};
80	md.g[0][0] = E; md.g[0][1] = F; md.g[1][0] = F; md.g[1][1] = G;
81	md.du[0][0] = dE_du; md.du[0][1] = dF_du; md.du[1][0] = dF_du; md.du[1][1] = dG_du;
82	md.dv[0][0] = dE_dv; md.dv[0][1] = dF_dv; md.dv[1][0] = dF_dv; md.dv[1][1] = dG_dv;
83	return md;
84	}
85
86	// Compute inverse g^{ij} and Christoffel symbols Gamma^k_{ij} for k,i,j in {0,1}
87	static inline void christoffel(double u, double v, double Gamma[2][2][2]){
88	MetricData md = buildMetric(u,v);
89	double E = md.g[0][0], F = md.g[0][1], G = md.g[1][1];
90	double det = EG - FF;
91	// Inverse metric
92	double inv[2][2];
93	inv[0][0] = G/det; inv[0][1] = -F/det;
94	inv[1][0] = -F/det; inv[1][1] = E/det;
95
96	// Helpers: partials of g_ij; index 0->u, 1->v
97	double dg[2][2][2];
98	dg[0][0][0] = md.du[0][0]; dg[1][0][0] = md.dv[0][0];
99	dg[0][0][1] = md.du[0][1]; dg[1][0][1] = md.dv[0][1];
100	dg[0][1][0] = md.du[1][0]; dg[1][1][0] = md.dv[1][0];
101	dg[0][1][1] = md.du[1][1]; dg[1][1][1] = md.dv[1][1];
102
103	// Christoffel: Gamma^k_{ij} = 1/2 g^{kl}( d_i g_{jl} + d_j g_{il} - d_l g_{ij} )
104	for(int k=0;k<2;k++){
105	for(int i=0;i<2;i++){
106	for(int j=0;j<2;j++){
107	double sum = 0.0;
108	for(int l=0;l<2;l++){
109	double term = dg[i][j][l] + dg[j][i][l] - dg[l][i][j];
110	sum += inv[k][l] * 0.5 * term;
111	}
112	Gamma[k][i][j] = sum;
113	}
114	}
115	}
116	}
117
118	// State y = [u, v, du, dv]; RHS dy/dt
119	static inline array<double,4> geodesic_rhs(double u, double v, double du, double dv){
120	double Gamma[2][2][2];
121	christoffel(u,v,Gamma);
122	// Second derivatives from geodesic ODE: d2u = -Gamma^u_{ij} dxi dxj, d2v = -Gamma^v_{ij} dxi dxj
123	double d2u = - (Gamma[0][0][0]dudu + 2.0Gamma[0][0][1]dudv + Gamma[0][1][1]dv*dv);
124	double d2v = - (Gamma[1][0][0]dudu + 2.0Gamma[1][0][1]dudv + Gamma[1][1][1]dv*dv);
125	return {du, dv, d2u, d2v};
126	}
127
128	// One RK4 step for the 4D system
129	static inline void rk4_step(array<double,4>& y, double h){
130	auto f = [](const array<double,4>& s){ return geodesic_rhs(s[0], s[1], s[2], s[3]); };
131	array<double,4> k1 = f(y);
132	array<double,4> y2{ y[0]+0.5hk1[0], y[1]+0.5hk1[1], y[2]+0.5hk1[2], y[3]+0.5hk1[3] };
133	array<double,4> k2 = f(y2);
134	array<double,4> y3{ y[0]+0.5hk2[0], y[1]+0.5hk2[1], y[2]+0.5hk2[2], y[3]+0.5hk2[3] };
135	array<double,4> k3 = f(y3);
136	array<double,4> y4{ y[0]+hk3[0], y[1]+hk3[1], y[2]+hk3[2], y[3]+hk3[3] };
137	array<double,4> k4 = f(y4);
138	for(int i=0;i<4;i++){
139	y[i] += (h/6.0) * (k1[i] + 2.0k2[i] + 2.0k3[i] + k4[i]);
140	}
141	}
142
143	int main(){
144	// Initial param position and tangent (choose small region where chart is regular)
145	double u0 = 0.0, v0 = 0.0; // apex of paraboloid
146	// Initial direction in parameters: choose unit-speed w.r.t. metric approximately
147	// We roughly normalize using the metric at start
148	MetricData md0 = buildMetric(u0, v0);
149	// Speed^2 = [du dv] g [du dv]^T. Choose raw (1, 0.2) then scale.
150	double du0 = 1.0, dv0 = 0.2;
151	double speed2 = md0.g[0][0]du0du0 + 2md0.g[0][1]du0dv0 + md0.g[1][1]dv0*dv0;
152	double scale = 1.0 / sqrt(max(1e-12, speed2));
153	du0 = scale; dv0 = scale; // roughly unit-speed
154
155	array<double,4> y{u0, v0, du0, dv0};
156
157	double h = 0.01; // step size
158	int steps = 1000;
159
160	cout.setf(ios::fixed); cout<<setprecision(6);
161	for(int i=0;i<steps;i++){
162	// Integrate one step
163	rk4_step(y, h);
164	// Get 3D point on surface
165	auto X3 = Surface::X(y[0], y[1]);
166	if(i % 50 == 0){
167	cout << "t=" << (i+1)*h << ": X=(" << X3[0] << ", " << X3[1] << ", " << X3[2] << ")\n";
168	}
169	}
170
171	return 0;
172	}
173