feat(geometry): pn_geometry.hpp — Pn projective-metric substrate (jReality port)

Ports the small but pervasive de.jreality.math.Pn surface that the Java
algorithmic core uses across every not-yet-ported geometric phase:
HyperbolicLayout, SphericalLayout, FundamentalPolygon (9c), KoebePolyhedron
(10c'), quasi-isothermic (10e), hyperelliptic theta, CircleDomain (11c).
Porting it once unblocks all downstream consumers instead of re-deriving
the metric ad hoc per phase.

pn_geometry.hpp — header-only, Eigen, ~120 LOC:
  PnMetric { EUCLIDEAN=0, ELLIPTIC=+1, HYPERBOLIC=-1 } matching jReality.
  pn_inner_product  — bilinear form per signature.
  pn_norm / pn_dehomogenize / pn_set_to_length / pn_normalize.
  pn_distance_between — Euclidean (spatial), elliptic (acos), hyperbolic (acosh).
  pn_linear_interpolation — affine (E) and slerp (S/H constant-speed geodesic).

HYPERBOLIC uses the timelike-positive ("upper-sheet") convention, identical to
the already-verified projective_math.hpp::hyperbolicDistance; pn_distance_between
is regression-anchored against it in the tests.

test_pn_geometry.cpp — 6 tests, all GREEN:
  InnerProductSignatures, EuclideanDistance, EllipticDistanceIsAngle,
  HyperbolicDistanceClosedFormAndAnchor (+ projective_math.hpp anchor),
  NormAndScaling, LinearInterpolationGeodesic.

doc/roadmap/java-parity.md — new "Infrastructure / support layers" section:
  de.jreality.math.Pn   → pn_geometry.hpp (partial, 6 fns covering core usage)
  de.jreality.math.Rn   → Eigen (no separate port needed)
  MatrixBuilder         → not yet (only needed for hyperelliptic theta)
  conformallab XML types → not planned (GUI persistence, not algorithmic)

246/246 cgal tests pass.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

code/include/pn_geometry.hpp

@@ -0,0 +1,156 @@
+#pragma once
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+
+// pn_geometry.hpp
+//
+// Projective-metric geometry substrate — faithful port of the
+// `de.jreality.math.Pn` surface that the Java algorithmic core
+// uses pervasively:
+//   HyperbolicLayout, SphericalLayout, FundamentalPolygon,
+//   KoebePolyhedron, quasi-isothermic utilities, hyperelliptic θ.
+//
+// Porting this ~30-function surface once unblocks every downstream
+// geometric phase (9c / 10c' / 10e / 11c) instead of re-deriving the
+// metric ad hoc per phase.
+//
+// ┌──────────────────────────────────────────────────────────────────┐
+// │  Metric signature convention (matches jReality exactly)          │
+// │                                                                  │
+// │  EUCLIDEAN  = 0   ⟨u,v⟩ = Σ_{i<last} uᵢvᵢ         (spatial)  │
+// │  ELLIPTIC   = +1  ⟨u,v⟩ = Σ_{all} uᵢvᵢ             (S^n)     │
+// │  HYPERBOLIC = -1  ⟨u,v⟩ = u_last·v_last − Σ_{i<last} uᵢvᵢ    │
+// │                           (Minkowski; last coord = timelike)    │
+// │                                                                  │
+// │  HYPERBOLIC uses the timelike-positive ("upper sheet") convention│
+// │  so a normalised hyperbolic point satisfies ⟨p,p⟩ = +1 and     │
+// │  ⟨p,q⟩ ≥ 1 → d(p,q) = acosh(⟨p,q⟩) directly.                 │
+// │  This is identical to projective_math.hpp::hyperbolicDistance,  │
+// │  against which pn_distance_between(..., HYPERBOLIC) is anchored. │
+// └──────────────────────────────────────────────────────────────────┘
+//
+// Vectors are Eigen column vectors (homogeneous, length = dim+1).
+
+#include <Eigen/Core>
+#include <cmath>
+#include <algorithm>
+
+namespace conformallab {
+
+/// Metric signatures matching `de.jreality.math.Pn`.
+enum PnMetric : int {
+    PN_EUCLIDEAN  =  0,
+    PN_ELLIPTIC   = +1,
+    PN_HYPERBOLIC = -1
+};
+
+// ── Inner product ─────────────────────────────────────────────────────────────
+
+/// Bilinear form ⟨u,v⟩ for the given metric.
+/// Spatial part = all but the last coordinate; the last coord is the
+/// homogeneous/timelike one.
+inline double pn_inner_product(const Eigen::VectorXd& u,
+                               const Eigen::VectorXd& v,
+                               int metric)
+{
+    const int n       = static_cast<int>(u.size());
+    const double spat = u.head(n - 1).dot(v.head(n - 1));
+    const double last = u(n - 1) * v(n - 1);
+    switch (metric) {
+        case PN_HYPERBOLIC: return  last - spat;
+        case PN_ELLIPTIC:   return  last + spat;
+        case PN_EUCLIDEAN:
+        default:            return  spat;
+    }
+}
+
+// ── Norm / scale ──────────────────────────────────────────────────────────────
+
+/// Metric norm √|⟨p,p⟩|.  Absolute value guards against tiny negative
+/// round-off in the Euclidean / hyperbolic degenerate cases.
+inline double pn_norm(const Eigen::VectorXd& p, int metric)
+{
+    return std::sqrt(std::abs(pn_inner_product(p, p, metric)));
+}
+
+/// Dehomogenise: divide by the last component (jReality `Pn.dehomogenize`).
+inline Eigen::VectorXd pn_dehomogenize(const Eigen::VectorXd& p)
+{
+    return p / p(p.size() - 1);
+}
+
+/// Scale `p` to metric norm `length` (jReality `Pn.setToLength`).
+/// Returns `p` unchanged when its norm is numerically zero.
+inline Eigen::VectorXd pn_set_to_length(const Eigen::VectorXd& p,
+                                        double length, int metric)
+{
+    const double nrm = pn_norm(p, metric);
+    if (nrm < 1e-300) return p;
+    return p * (length / nrm);
+}
+
+/// Normalise to unit metric norm (jReality `Pn.normalize`).
+inline Eigen::VectorXd pn_normalize(const Eigen::VectorXd& p, int metric)
+{
+    return pn_set_to_length(p, 1.0, metric);
+}
+
+// ── Distance ──────────────────────────────────────────────────────────────────
+
+/// Geodesic distance between two points (jReality `Pn.distanceBetween`).
+///   EUCLIDEAN:  ‖p̂_spatial − q̂_spatial‖ (dehomogenised)
+///   ELLIPTIC:   acos(⟨p̂,q̂⟩)  (spherical angle, clamped to [−1,1])
+///   HYPERBOLIC: acosh(⟨p̂,q̂⟩)  (clamped ≥ 1)
+inline double pn_distance_between(const Eigen::VectorXd& p,
+                                  const Eigen::VectorXd& q,
+                                  int metric)
+{
+    if (metric == PN_EUCLIDEAN) {
+        const auto   pd = pn_dehomogenize(p);
+        const auto   qd = pn_dehomogenize(q);
+        const int    n  = static_cast<int>(pd.size());
+        return (pd.head(n - 1) - qd.head(n - 1)).norm();
+    }
+    const double np = pn_norm(p, metric);
+    const double nq = pn_norm(q, metric);
+    double c = pn_inner_product(p, q, metric) / (np * nq);
+    if (metric == PN_HYPERBOLIC)
+        return std::acosh(std::max(1.0, c));
+    // ELLIPTIC
+    c = std::clamp(c, -1.0, 1.0);
+    return std::acos(c);
+}
+
+// ── Geodesic interpolation ────────────────────────────────────────────────────
+
+/// Constant-speed geodesic interpolation at parameter t ∈ [0,1]
+/// (jReality `Pn.linearInterpolation`).
+///   EUCLIDEAN:  affine blend of the dehomogenised points.
+///   ELLIPTIC:   spherical slerp:
+///                 r(t) = (sin((1−t)d)·p̂ + sin(td)·q̂) / sin(d)
+///   HYPERBOLIC: hyperbolic slerp:
+///                 r(t) = (sinh((1−t)d)·p̂ + sinh(td)·q̂) / sinh(d)
+///   where d = pn_distance_between(p,q,metric) and p̂,q̂ are unit vectors.
+inline Eigen::VectorXd pn_linear_interpolation(const Eigen::VectorXd& p,
+                                               const Eigen::VectorXd& q,
+                                               double t, int metric)
+{
+    if (metric == PN_EUCLIDEAN) {
+        const auto pd = pn_dehomogenize(p);
+        const auto qd = pn_dehomogenize(q);
+        return (1.0 - t) * pd + t * qd;
+    }
+    const auto   ph = pn_normalize(p, metric);
+    const auto   qh = pn_normalize(q, metric);
+    const double d  = pn_distance_between(ph, qh, metric);
+    if (d < 1e-12) return ph;          // coincident — return either endpoint
+    if (metric == PN_ELLIPTIC) {
+        const double s = std::sin(d);
+        return (std::sin((1.0 - t) * d) * ph + std::sin(t * d) * qh) / s;
+    }
+    // HYPERBOLIC
+    const double s = std::sinh(d);
+    return (std::sinh((1.0 - t) * d) * ph + std::sinh(t * d) * qh) / s;
+}
+
+} // namespace conformallab

code/tests/cgal/CMakeLists.txt

@@ -86,6 +86,11 @@ add_executable(conformallab_cgal_tests
     #    newton_inversive_distance (FD Hessian).
     test_newton_phase9a.cpp
 
+    # ── Pn projective-metric substrate (de.jreality.math.Pn port) ───────────
+    #    Inner product / norm / distance / geodesic interpolation in
+    #    Euclidean, Elliptic, Hyperbolic signatures.
+    test_pn_geometry.cpp
+
     # ── Phase 8b-Lite: CGAL entry wrappers for the 4 non-Euclidean modes ─────
     #    Spherical, HyperIdeal, CircleP-Euclidean, Inversive-Distance via
     #    <CGAL/Discrete_*.h> public API + Conformal_layout.h wrapper.

code/tests/cgal/test_pn_geometry.cpp

@@ -0,0 +1,135 @@
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+
+// test_pn_geometry.cpp
+//
+// Verifies pn_geometry.hpp against:
+//   (a) closed-form analytic golden values, and
+//   (b) the already-verified projective_math.hpp::hyperbolicDistance
+//       (regression anchor for the HYPERBOLIC sign convention).
+//
+// Mirrored Java source: de.jreality.math.Pn
+
+#include "pn_geometry.hpp"
+#include "projective_math.hpp"
+#include <gtest/gtest.h>
+#include <Eigen/Core>
+#include <cmath>
+
+using namespace conformallab;
+
+namespace {
+Eigen::VectorXd V(std::initializer_list<double> xs) {
+    Eigen::VectorXd v(static_cast<int>(xs.size()));
+    int i = 0;
+    for (double x : xs) v(i++) = x;
+    return v;
+}
+constexpr double PN_PI = 3.14159265358979323846;
+} // namespace
+
+// ── Inner product per signature ──────────────────────────────────────────────
+TEST(PnGeometry, InnerProductSignatures)
+{
+    auto u = V({1.0, 2.0, 3.0});   // spatial=(1,2), last=3
+    auto v = V({4.0, 5.0, 6.0});   // spatial=(4,5), last=6
+    // spat = 1·4 + 2·5 = 14,  last = 3·6 = 18
+    EXPECT_NEAR(pn_inner_product(u, v, PN_EUCLIDEAN),   14.0,      1e-12);
+    EXPECT_NEAR(pn_inner_product(u, v, PN_ELLIPTIC),    14.0+18.0, 1e-12);
+    EXPECT_NEAR(pn_inner_product(u, v, PN_HYPERBOLIC),  18.0-14.0, 1e-12);
+}
+
+// ── Euclidean distance ────────────────────────────────────────────────────────
+TEST(PnGeometry, EuclideanDistance)
+{
+    // 3-4-5 right triangle in the plane (homogeneous w=1).
+    auto p  = V({0.0, 0.0, 1.0});
+    auto q  = V({3.0, 4.0, 1.0});
+    EXPECT_NEAR(pn_distance_between(p, q, PN_EUCLIDEAN), 5.0, 1e-12);
+
+    // Homogeneous scaling must not change the distance.
+    auto q2 = V({6.0, 8.0, 2.0});  // same point as q after dehomogenize
+    EXPECT_NEAR(pn_distance_between(p, q2, PN_EUCLIDEAN), 5.0, 1e-12);
+}
+
+// ── Elliptic distance = spherical angle ───────────────────────────────────────
+TEST(PnGeometry, EllipticDistanceIsAngle)
+{
+    auto e1 = V({1.0, 0.0, 0.0});
+    auto e2 = V({0.0, 1.0, 0.0});
+    auto d  = V({1.0, 1.0, 0.0});   // 45° between e1 and e2
+
+    EXPECT_NEAR(pn_distance_between(e1, e2, PN_ELLIPTIC), PN_PI / 2.0, 1e-12);
+    EXPECT_NEAR(pn_distance_between(e1, d,  PN_ELLIPTIC), PN_PI / 4.0, 1e-12);
+    EXPECT_NEAR(pn_distance_between(e1, e1, PN_ELLIPTIC), 0.0,      1e-12);
+}
+
+// ── Hyperbolic distance: closed form + projective_math.hpp anchor ─────────────
+TEST(PnGeometry, HyperbolicDistanceClosedFormAndAnchor)
+{
+    // Upper hyperboloid: p = apex (0,0,1); q = (sinh r, 0, cosh r) at distance r.
+    const double r = 0.873;
+    auto p  = V({0.0, 0.0, 1.0});
+    auto q  = V({std::sinh(r), 0.0, std::cosh(r)});
+
+    EXPECT_NEAR(pn_distance_between(p, q, PN_HYPERBOLIC), r, 1e-10);
+
+    // Regression anchor: identical to projective_math.hpp::hyperbolicDistance.
+    EXPECT_NEAR(pn_distance_between(p, q, PN_HYPERBOLIC),
+                hyperbolicDistance(p, q), 1e-12);
+
+    // Scale-invariance: homogeneous rescaling must not change distance.
+    auto q2 = (2.5 * q).eval();
+    EXPECT_NEAR(pn_distance_between(p, q2, PN_HYPERBOLIC), r, 1e-10);
+}
+
+// ── Norm / setToLength / normalize ───────────────────────────────────────────
+TEST(PnGeometry, NormAndScaling)
+{
+    auto p = V({3.0, 4.0, 0.0});
+    EXPECT_NEAR(pn_norm(p, PN_ELLIPTIC), 5.0, 1e-12);
+
+    auto s = pn_set_to_length(p, 2.0, PN_ELLIPTIC);
+    EXPECT_NEAR(pn_norm(s, PN_ELLIPTIC), 2.0, 1e-12);
+
+    auto u = pn_normalize(p, PN_ELLIPTIC);
+    EXPECT_NEAR(pn_norm(u, PN_ELLIPTIC), 1.0, 1e-12);
+
+    // A unit hyperboloid sheet point already has hyperbolic norm 1.
+    auto h = V({std::sinh(0.6), 0.0, std::cosh(0.6)});
+    EXPECT_NEAR(pn_norm(h, PN_HYPERBOLIC), 1.0, 1e-12);
+}
+
+// ── Geodesic interpolation ────────────────────────────────────────────────────
+TEST(PnGeometry, LinearInterpolationGeodesic)
+{
+    // EUCLIDEAN: affine midpoint.
+    {
+        auto p = V({0.0, 0.0, 1.0});
+        auto q = V({4.0, 0.0, 1.0});
+        auto m = pn_linear_interpolation(p, q, 0.5, PN_EUCLIDEAN);
+        EXPECT_NEAR(m(0), 2.0, 1e-12);
+    }
+
+    // ELLIPTIC: midpoint of e1,e2 is at half the π/2 arc = π/4 from each.
+    {
+        auto e1 = V({1.0, 0.0, 0.0});
+        auto e2 = V({0.0, 1.0, 0.0});
+        auto m  = pn_linear_interpolation(e1, e2, 0.5, PN_ELLIPTIC);
+        EXPECT_NEAR(pn_distance_between(e1, m, PN_ELLIPTIC), PN_PI / 4.0, 1e-10);
+        EXPECT_NEAR(pn_distance_between(e2, m, PN_ELLIPTIC), PN_PI / 4.0, 1e-10);
+        // Endpoint recovery at t=0.
+        auto m0 = pn_linear_interpolation(e1, e2, 0.0, PN_ELLIPTIC);
+        EXPECT_NEAR(pn_distance_between(e1, m0, PN_ELLIPTIC), 0.0, 1e-10);
+    }
+
+    // HYPERBOLIC: midpoint is at exactly half the geodesic distance from each end.
+    {
+        const double r = 1.2;
+        auto p = V({0.0, 0.0, 1.0});
+        auto q = V({std::sinh(r), 0.0, std::cosh(r)});
+        auto m = pn_linear_interpolation(p, q, 0.5, PN_HYPERBOLIC);
+        EXPECT_NEAR(pn_distance_between(p, m, PN_HYPERBOLIC), r / 2.0, 1e-9);
+        EXPECT_NEAR(pn_distance_between(q, m, PN_HYPERBOLIC), r / 2.0, 1e-9);
+    }
+}

doc/roadmap/java-parity.md

@@ -43,6 +43,20 @@ as the reference implementation for expected behaviour, edge cases, and test cas
 
 ---
 
+## Infrastructure / support layers
+
+These are not algorithmic phases but shared geometry substrates that multiple phases
+depend on.  Porting them once unblocks all downstream consumers.
+
+| Java library | What it provides | C++ status | Notes |
+|---|---|---|---|
+| **`de.jreality.math.Pn`** | Projective-metric primitives: inner product, norm, distance, geodesic interpolation in Euclidean / Elliptic / Hyperbolic signatures | ⚠️ **partial** → `pn_geometry.hpp` (2026-05-30) | 6 functions (~30 LOC) cover the entire surface used by the Java algorithmic core (27 files, ≈106 imports). Unblocks: HyperbolicLayout, SphericalLayout, FundamentalPolygon (9c), KoebePolyhedron (10c'), quasi-isothermic (10e), hyperelliptic θ, CircleDomain (11c). |
+| **`de.jreality.math.Rn`** | Real vector / matrix utilities (add, sub, scale, dot, norm, …) | ⚠️ partial → Eigen | The Rn surface used in the core (16 distinct methods) maps directly onto Eigen VectorXd / MatrixXd operations. No separate port needed; call-sites replace `Rn.foo(a,b)` with Eigen expressions inline. |
+| **`de.jreality.math.MatrixBuilder`** | Isometry construction: `rotateFromTo`, `translateFromTo` in Euclidean / Hyperbolic models | ❌ not yet | Used only by `HyperIdealHyperellipticUtility` (hyperelliptic θ). Port on demand when that utility is needed. |
+| **conformallab XML types** (`de.varylab.conformallab.types.*`) | JAXB-generated persistence format for Java GUI state (`HalfedgeEmbedding`, `UniformizationData`, …) | ❌ not planned | GUI / persistence layer, not an algorithmic substrate. conformallab++ has its own serialisation (`serialization.hpp`). Asset-specific parsers (e.g. for `lawson_curve_source.xml`) are added on demand, not as a general port. |
+
+---
+
 ## Java utility classes not yet ported
 
 These exist in `de.varylab.discreteconformal.util` in the Java library.