tests: port DiscreteEllipticUtility + P2 tests; stub HDS-blocked tests

Fully ported (pure math, no HDS required):
  test_discrete_elliptic_utility.cpp  – 2 tests
    normalizeModulus: move tau into SL(2,Z) fundamental domain
  test_p2_utility.cpp                 – 3 tests
    P2 projective geometry (perpendicularBisector, pointFromLines,
    makeDirectIsometryFromFrames double vs long double precision)

New headers:
  include/discrete_elliptic_utility.hpp  – normalizeModulus
  include/p2_utility.hpp                 – P2 Euclidean geometry (templated
    on scalar type so double and long double share one implementation)

Stubs (GTEST_SKIP, blocked until HDS port – Phase 4):
  test_hyper_ideal_functional.cpp          – 5 tests (1 @Ignore in Java)
  test_hyper_ideal_hyperelliptic_utility.cpp – 3 tests
  test_spherical_functional.cpp            – 5 tests
  All use CoHDS + HalfEdgeUtils which are not yet ported to C++.

Result: 34 tests total | 21 passed | 13 skipped | 0 failed

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

code/include/discrete_elliptic_utility.hpp

@@ -0,0 +1,43 @@
+#pragma once
+
+// Ported from de.varylab.discreteconformal.util.DiscreteEllipticUtility (Java).
+// Only the pure-math subset (no HDS required).
+
+#include <complex>
+#include <cmath>
+
+namespace conformallab {
+
+// Move tau into the fundamental domain of the modular group SL(2,Z):
+//   |Re(tau)| <= 0.5,  Im(tau) >= 0,  Re(tau) >= 0,  |tau| >= 1
+//
+// Algorithm: iteratively apply
+//   1. T-shift:   Re > 0.5 or Re < 0  →  Re -= sign(Re)
+//   2. Im-flip:   Im < 0              →  Im = -Im
+//   3. Re-flip:   Re < 0              →  Re = -Re
+//   4. S-invert:  |tau| < 1           →  tau = 1/tau
+//
+// Corresponds to Java DiscreteEllipticUtility.normalizeModulus(Complex).
+inline std::complex<double> normalizeModulus(std::complex<double> tau) {
+    int maxIter = 100;
+    while (--maxIter > 0) {
+        double re = tau.real();
+        double im = tau.imag();
+        // exit when all conditions satisfied
+        if (std::abs(re) <= 0.5 && im >= 0.0 && re >= 0.0 && std::abs(tau) >= 1.0)
+            break;
+
+        if (std::abs(re) > 0.5)
+            re -= (re > 0.0 ? 1.0 : -1.0);   // signum shift
+        if (im < 0.0)
+            im = -im;
+        if (re < 0.0)
+            re = -re;
+        tau = std::complex<double>(re, im);
+        if (std::abs(tau) < 1.0)
+            tau = 1.0 / tau;                  // S-transformation: invert
+    }
+    return tau;
+}
+
+} // namespace conformallab

code/include/p2_utility.hpp

@@ -0,0 +1,102 @@
+#pragma once
+
+// 2-D projective geometry utilities for the Euclidean signature.
+// Ported from de.jreality.math.P2 and de.varylab.discreteconformal.math.P2Big.
+//
+// Points and lines are represented as homogeneous 3-vectors (x, y, w).
+// In the Euclidean case a finite point (px, py) is stored as (px, py, 1).
+
+#include <Eigen/Dense>
+#include <cmath>
+
+namespace conformallab {
+
+// ── Point / line duality ──────────────────────────────────────────────────────
+
+// Intersection of two lines l1, l2  (or  line through two points p1, p2)
+// via the cross product.  Works for any P2 element.
+// Corresponds to Java P2.pointFromLines / P2.lineFromPoints.
+inline Eigen::Vector3d pointFromLines(const Eigen::Vector3d& l1,
+                                      const Eigen::Vector3d& l2) {
+    return l1.cross(l2);
+}
+
+// ── Euclidean perpendicular bisector ─────────────────────────────────────────
+
+// Returns the homogeneous line coordinates (a, b, c) of the perpendicular
+// bisector of the segment [p, q] in the Euclidean plane.
+// Coordinates: ax + by + c = 0  (after dehomogenizing p and q).
+//
+// Corresponds to Java P2.perpendicularBisector(p, q, Pn.EUCLIDEAN).
+inline Eigen::Vector3d perpendicularBisectorEuclidean(const Eigen::Vector3d& p_h,
+                                                      const Eigen::Vector3d& q_h) {
+    // Dehomogenize
+    Eigen::Vector2d p = p_h.head<2>() / p_h(2);
+    Eigen::Vector2d q = q_h.head<2>() / q_h(2);
+
+    // Direction vector (p → direction, matching jReality sign convention)
+    Eigen::Vector2d d = p - q;
+
+    // Midpoint
+    Eigen::Vector2d m = (p + q) * 0.5;
+
+    // Line: d[0]*(x - m[0]) + d[1]*(y - m[1]) = 0
+    //      = d[0]*x + d[1]*y - (d[0]*m[0] + d[1]*m[1])
+    double c = -(d(0) * m(0) + d(1) * m(1));
+    return {d(0), d(1), c};
+}
+
+// ── Euclidean distance between two P2 homogeneous points ─────────────────────
+
+inline double euclideanDistanceP2(const Eigen::Vector3d& p_h,
+                                  const Eigen::Vector3d& q_h) {
+    Eigen::Vector2d p = p_h.head<2>() / p_h(2);
+    Eigen::Vector2d q = q_h.head<2>() / q_h(2);
+    return (p - q).norm();
+}
+
+// ── Direct Euclidean isometry from two point-frames ──────────────────────────
+
+// Build the 3×3 projective matrix that represents the coordinate frame
+// anchored at p0 with p1 defining the positive x-direction.
+// Euclidean case: columns are [dehom(p0), unit_dir(p0→p1), perp_dir].
+//
+// Template parameter S allows float / double / long double.
+template <typename S>
+Eigen::Matrix<S, 3, 3> makeFrameMatrix(Eigen::Matrix<S, 3, 1> p0_h,
+                                        Eigen::Matrix<S, 3, 1> p1_h) {
+    // Dehomogenize
+    Eigen::Matrix<S, 3, 1> p0 = p0_h / p0_h(2);          // (px, py, 1)
+    Eigen::Matrix<S, 3, 1> p1_d = p1_h / p1_h(2);
+
+    // Unit direction p0 → p1
+    Eigen::Matrix<S, 2, 1> dir2 = (p1_d - p0).template head<2>();
+    dir2.normalize();
+    Eigen::Matrix<S, 3, 1> p1n(dir2(0), dir2(1), S(0));
+
+    // Perpendicular direction
+    Eigen::Matrix<S, 3, 1> p2(-dir2(1), dir2(0), S(0));
+
+    Eigen::Matrix<S, 3, 3> M;
+    M.col(0) = p0;
+    M.col(1) = p1n;
+    M.col(2) = p2;
+    return M;
+}
+
+// Find the 3×3 Euclidean isometry (as a projective matrix) that maps
+// the frame (s1, s2) to the frame (t1, t2).
+//
+// Corresponds to Java P2.makeDirectIsometryFromFrames(s1, s2, t1, t2, Pn.EUCLIDEAN)
+// and P2Big.makeDirectIsometryFromFrames(...) (the BigDecimal / high-precision variant).
+template <typename S>
+Eigen::Matrix<S, 3, 3> makeDirectIsometryFromFramesEuclidean(
+    Eigen::Matrix<S, 3, 1> s1, Eigen::Matrix<S, 3, 1> s2,
+    Eigen::Matrix<S, 3, 1> t1, Eigen::Matrix<S, 3, 1> t2)
+{
+    auto toS = makeFrameMatrix<S>(s1, s2);
+    auto toT = makeFrameMatrix<S>(t1, t2);
+    return toT * toS.inverse();
+}
+
+} // namespace conformallab

code/tests/CMakeLists.txt

@@ -1,8 +1,17 @@
 add_executable(conformallab_tests
+    # ── Fully ported (pure math, no HDS) ────────────────────────────────────
     test_clausen.cpp
     test_hyper_ideal_utility.cpp
     test_matrix_utility.cpp
     test_surface_curve_utility.cpp
+    test_discrete_elliptic_utility.cpp
+    test_p2_utility.cpp
+
+    # ── Stubs: blocked until HDS port (Phase 4) ──────────────────────────────
+    # All tests call GTEST_SKIP() with a clear explanation.
+    test_hyper_ideal_functional.cpp
+    test_hyper_ideal_hyperelliptic_utility.cpp
+    test_spherical_functional.cpp
 )
 
 target_include_directories(conformallab_tests SYSTEM PRIVATE

code/tests/test_discrete_elliptic_utility.cpp

@@ -0,0 +1,39 @@
+// Port of de.varylab.discreteconformal.util.DiscreteEllipticUtilityTest (Java/JUnit).
+// Tests the normalizeModulus function that moves a complex number tau into the
+// fundamental domain of the modular group SL(2,Z).
+
+#include "discrete_elliptic_utility.hpp"
+#include <gtest/gtest.h>
+#include <complex>
+#include <cmath>
+
+using namespace conformallab;
+
+// Corresponds to Java testNormalizeModulus()
+TEST(DiscreteEllipticUtilityTest, NormalizeModulus) {
+    // tau already in fundamental domain → should be returned unchanged
+    std::complex<double> tau(0.45, 1.1);
+    auto tauNorm = normalizeModulus(tau);
+    EXPECT_NEAR(0.45, tauNorm.real(), 1E-12);
+    EXPECT_NEAR(1.1,  tauNorm.imag(), 1E-12);
+
+    // tau = i/3  (|tau| < 1) → inversion gives 3i
+    tau = std::complex<double>(0.0, 1.0 / 3.0);
+    tauNorm = normalizeModulus(tau);
+    EXPECT_NEAR(3.0, tauNorm.imag(), 1E-12);
+    EXPECT_NEAR(0.0, tauNorm.real(), 1E-12);
+}
+
+// Corresponds to Java testNormalizeModulusPeriodShift()
+// Two tau values that differ by a T-shift (integer shift of Re) must normalize
+// to the same point in the fundamental domain.
+TEST(DiscreteEllipticUtilityTest, NormalizeModulusPeriodShift) {
+    std::complex<double> tau1(0.3,  1.0);
+    std::complex<double> tau2(-0.7, 1.0);   // tau2 = tau1 - 1
+
+    auto n1 = normalizeModulus(tau1);
+    auto n2 = normalizeModulus(tau2);
+
+    EXPECT_NEAR(n1.real(), n2.real(), 1E-12) << "real parts should be equal";
+    EXPECT_NEAR(n1.imag(), n2.imag(), 1E-12) << "imag parts should be equal";
+}

code/tests/test_hyper_ideal_functional.cpp

@@ -0,0 +1,37 @@
+// Stub for de.varylab.discreteconformal.functional.HyperIdealFunctionalTest (Java/JUnit).
+//
+// STATUS: BLOCKED – requires HDS port (Phase 4).
+//
+// These tests evaluate gradient and Hessian of the HyperIdealFunctional on
+// actual mesh data (CoHDS + HyperIdealGenerator).  They cannot be ported
+// until the HalfEdge data structure (CoHDS), the functional evaluation
+// framework, and the mesh generators are available in C++.
+//
+// Java tests and their status:
+//   testHessian()                           – @Ignore in Java (skipped here too)
+//   testGradientWithHyperIdealAndIdealPoints – blocked: needs HDS
+//   testGradientInTheExtendedDomain          – blocked: needs HDS
+//   testGradientWithHyperellipticCurve       – blocked: needs HDS
+//   testFunctionalAtNaNValue                 – blocked: needs HDS
+
+#include <gtest/gtest.h>
+
+TEST(HyperIdealFunctionalTest, TestHessian_IgnoredInJava) {
+    GTEST_SKIP() << "@Ignore in Java – skipped here too";
+}
+
+TEST(HyperIdealFunctionalTest, GradientWithHyperIdealAndIdealPoints) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)";
+}
+
+TEST(HyperIdealFunctionalTest, GradientInTheExtendedDomain) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)";
+}
+
+TEST(HyperIdealFunctionalTest, GradientWithHyperellipticCurve) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)";
+}
+
+TEST(HyperIdealFunctionalTest, FunctionalAtNaNValue) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealFunctional)";
+}

code/tests/test_hyper_ideal_hyperelliptic_utility.cpp

@@ -0,0 +1,26 @@
+// Stub for de.varylab.discreteconformal.functional.HyperIdealHyperellipticUtilityTest.
+//
+// STATUS: BLOCKED – requires HDS port (Phase 4).
+//
+// Tests compute intersection angles of circles associated with hyper-ideal
+// vertices using CoHDS + HalfEdgeUtils.  All three tests operate on mesh
+// data structures that are not yet available in C++.
+//
+// Java tests and their status:
+//   testCalculateCircleIntersections         – blocked: needs CoHDS + HalfEdgeUtils
+//   testCalculateCircleIntersectionsInfinite – blocked: needs CoHDS + HalfEdgeUtils
+//   testLawsonHyperellipticAngles            – blocked: needs CoHDS + HyperIdealGenerator
+
+#include <gtest/gtest.h>
+
+TEST(HyperIdealHyperellipticUtilityTest, CalculateCircleIntersections) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HalfEdgeUtils)";
+}
+
+TEST(HyperIdealHyperellipticUtilityTest, CalculateCircleIntersectionsInfinite) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HalfEdgeUtils)";
+}
+
+TEST(HyperIdealHyperellipticUtilityTest, LawsonHyperellipticAngles) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + HyperIdealGenerator)";
+}

code/tests/test_p2_utility.cpp

@@ -0,0 +1,87 @@
+// Port of de.varylab.discreteconformal.math.P2BigTest (Java/JUnit).
+// Tests 2-D projective geometry utilities: perpendicular bisectors,
+// point-from-lines, and direct isometries in the Euclidean plane.
+//
+// The Java test compared double precision (P2) against BigDecimal precision
+// (P2Big) to 1E-10.  Here we compare double against long double to the
+// same tolerance.
+
+#include "p2_utility.hpp"
+#include <gtest/gtest.h>
+#include <Eigen/Dense>
+#include <cmath>
+
+using namespace conformallab;
+
+// Corresponds to Java P2BigTest.testMakeDirectIsometryFromFramesEuclidean()
+//
+// Computes the Euclidean isometry mapping frame (s1,s2) to frame (t1,t2)
+// with both double and long-double precision, and checks:
+//   1. The two precisions agree to 1E-10 (precision stability).
+//   2. The matrix actually maps s1→t1 and s2→t2.
+TEST(P2UtilityTest, MakeDirectIsometryFromFramesEuclidean) {
+    using V3d  = Eigen::Vector3d;
+    using V3ld = Eigen::Matrix<long double, 3, 1>;
+
+    V3d s1(-1.4142135623730963, 0.0, 1.0);
+    V3d s2( 1.4142135623730951, 0.0, 1.0);
+    V3d t1(-2.828427124746189,  2.4494897427831805, 1.0);
+    V3d t2( 0.0,               2.4494897427831783, 1.0);
+
+    // double precision
+    auto T = makeDirectIsometryFromFramesEuclidean<double>(s1, s2, t1, t2);
+
+    // long double precision (analogous to Java's BigDecimal P2Big)
+    V3ld s1l = s1.cast<long double>();
+    V3ld s2l = s2.cast<long double>();
+    V3ld t1l = t1.cast<long double>();
+    V3ld t2l = t2.cast<long double>();
+    auto Tl = makeDirectIsometryFromFramesEuclidean<long double>(s1l, s2l, t1l, t2l);
+
+    // 1. double vs long double must agree to 1E-10
+    for (int i = 0; i < 3; ++i)
+        for (int j = 0; j < 3; ++j)
+            EXPECT_NEAR((double)Tl(i,j), T(i,j), 1E-10)
+                << "element (" << i << "," << j << ") differs between precisions";
+
+    // 2. T must map s1 → t1 and s2 → t2  (verify isometry correctness)
+    auto map_s1 = T * s1;
+    auto map_s2 = T * s2;
+    EXPECT_NEAR(euclideanDistanceP2(map_s1, t1), 0.0, 1E-9) << "T*s1 should equal t1";
+    EXPECT_NEAR(euclideanDistanceP2(map_s2, t2), 0.0, 1E-9) << "T*s2 should equal t2";
+}
+
+// Corresponds to Java P2BigTest.testPerpendicularBisector()
+TEST(P2UtilityTest, PerpendicularBisector) {
+    Eigen::Vector3d p1(0.5, 0.0, 1.0);
+    Eigen::Vector3d q1(0.0, 0.5, 1.0);
+
+    auto bisector = perpendicularBisectorEuclidean(p1, q1);
+
+    EXPECT_NEAR( 0.5, bisector(0), 1E-10);
+    EXPECT_NEAR(-0.5, bisector(1), 1E-10);
+    EXPECT_NEAR( 0.0, bisector(2), 1E-10);
+}
+
+// Corresponds to Java P2BigTest.testPerpendicularBisectorIntersection()
+//
+// The intersection of the perpendicular bisectors of two edges must be
+// equidistant from the endpoints of each edge (circumcenter property).
+TEST(P2UtilityTest, PerpendicularBisectorIntersection) {
+    Eigen::Vector3d p1(0.5, 0.0, 1.0);
+    Eigen::Vector3d q1(0.0, 1.0, 1.0);
+    Eigen::Vector3d p2(1.0, 0.0, 1.0);
+    Eigen::Vector3d q2(0.0, 1.5, 1.0);
+
+    auto l1 = perpendicularBisectorEuclidean(p1, q1);
+    auto l2 = perpendicularBisectorEuclidean(p2, q2);
+    auto o  = pointFromLines(l1, l2);    // circumcenter
+
+    // o must be equidistant from p1 and q1
+    EXPECT_NEAR(euclideanDistanceP2(p1, o),
+                euclideanDistanceP2(q1, o), 1E-10);
+
+    // o must be equidistant from p2 and q2
+    EXPECT_NEAR(euclideanDistanceP2(p2, o),
+                euclideanDistanceP2(q2, o), 1E-10);
+}

code/tests/test_spherical_functional.cpp

@@ -0,0 +1,37 @@
+// Stub for de.varylab.discreteconformal.functional.SphericalFunctionalTest (Java/JUnit).
+//
+// STATUS: BLOCKED – requires HDS port (Phase 4).
+//
+// Tests evaluate gradient and Hessian of the SphericalFunctional on meshes
+// built via CoHDS + ConvexHull, and check that a regular spherical metric
+// is a critical point of the functional.  All tests require the HalfEdge
+// data structure and the functional evaluation framework in C++.
+//
+// Java tests and their status:
+//   testReducedGradient   – blocked: needs CoHDS + SphericalFunctional
+//   testReducedHessian    – blocked: needs CoHDS + SphericalFunctional
+//   testGradient          – blocked: needs CoHDS + SphericalFunctional
+//   testHessian           – blocked: needs CoHDS + SphericalFunctional
+//   testCriticalPoint     – blocked: needs CoHDS + ConvexHull + SphericalFunctional
+
+#include <gtest/gtest.h>
+
+TEST(SphericalFunctionalTest, ReducedGradient) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)";
+}
+
+TEST(SphericalFunctionalTest, ReducedHessian) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)";
+}
+
+TEST(SphericalFunctionalTest, Gradient) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)";
+}
+
+TEST(SphericalFunctionalTest, Hessian) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + SphericalFunctional)";
+}
+
+TEST(SphericalFunctionalTest, CriticalPoint) {
+    GTEST_SKIP() << "Blocked: requires HDS port (CoHDS + ConvexHull + SphericalFunctional)";
+}