ConformalLabpp/code/include/period_matrix.hpp

#pragma once
// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// period_matrix.hpp
//
// Phase 7 — Period matrix for closed surfaces with Euclidean (flat) metric.
//
// For a closed genus-g surface with Euclidean conformal structure the holonomy
// group is generated by 2g translations ω_1, ..., ω_{2g} ∈ ℂ ≅ ℝ².
//
// ─── Genus-1 (flat torus) ────────────────────────────────────────────────────
//
//   The lattice Λ = ℤ·ω_1 ⊕ ℤ·ω_2 determines the conformal type.
//
//   Period ratio:  τ = ω_2 / ω_1   (as complex numbers)
//
//   By convention choose ω_1 such that Im(τ) > 0.
//   The conformal modulus / Teichmüller parameter is the SL(2,ℤ)-orbit of τ.
//
//   Reduction to fundamental domain {|τ| ≥ 1, −½ ≤ Re(τ) < ½, Im(τ) > 0}:
//     S:  τ ↦ −1/τ        (inversion)
//     T:  τ ↦  τ + 1      (translation)
//   Apply S and T repeatedly until τ is in the fundamental domain.
//
// ─── Genus g > 1 ─────────────────────────────────────────────────────────────
//
//   The full period matrix is a g×g complex symmetric matrix Ω with positive
//   definite imaginary part (Siegel upper half-space H_g).
//   Computing Ω from holonomy data requires integration of holomorphic
//   differentials — not implemented here.  For g > 1, this function returns
//   only the 2×2 block for the first pair of generators.
//
// ─── API ─────────────────────────────────────────────────────────────────────
//
//   PeriodData pd = compute_period_matrix(holonomy);
//   pd.tau          — complex period ratio τ (genus 1)
//   pd.omega        — holonomy generators as complex numbers (size = 2g)
//   pd.in_fundamental_domain — whether τ has been reduced
//
//   std::complex<double> reduce_to_fundamental_domain(τ)  — apply SL(2,ℤ)

#include "layout.hpp"
#include "discrete_elliptic_utility.hpp"   // normalizeModulus (Java-faithful)
#include <complex>
#include <cmath>
#include <vector>
#include <stdexcept>
#include <sstream>

namespace conformallab {

// ─────────────────────────────────────────────────────────────────────────────
// PeriodData
// ─────────────────────────────────────────────────────────────────────────────

/// Period-matrix data for a genus-g closed surface.  For genus 1 the
/// conformal type is fully captured by `τ = ω₂ / ω₁ ∈ ℍ`.
struct PeriodData {
    /// Lattice generators as complex numbers (one per cut edge).
    /// omega[i] = translations[i].x() + i·translations[i].y()
    std::vector<std::complex<double>> omega;

    /// Period ratio τ = omega[1] / omega[0]  (genus-1 only).
    /// Undefined (NaN) for genus != 1 or if holonomy has fewer than 2 generators.
    std::complex<double> tau = std::complex<double>(
        std::numeric_limits<double>::quiet_NaN(), 0.0);

    /// True if τ has been reduced to the standard fundamental domain.
    bool in_fundamental_domain = false;

    /// Genus of the surface = `|omega| / 2`.
    int genus() const { return static_cast<int>(omega.size()) / 2; }
};

// ─────────────────────────────────────────────────────────────────────────────
// reduce_to_fundamental_domain
//
// Applies SL(2,ℤ) generators S: τ↦−1/τ and T: τ↦τ+1 to bring τ into
//   F = { τ ∈ ℍ : |τ| ≥ 1,  −½ ≤ Re(τ) < ½ }
//
// Returns the reduced τ.  Throws if Im(τ) ≤ 0 (not in upper half-plane).
// ─────────────────────────────────────────────────────────────────────────────
/// Reduce `τ ∈ ℍ` to the standard SL(2,ℤ) fundamental domain
/// `F = { τ ∈ ℍ : |τ| ≥ 1, −½ ≤ Re τ < ½ }` via the generators
/// `S: τ↦−1/τ` and `T: τ↦τ+1`.  Throws if `Im τ ≤ 0`.
inline std::complex<double> reduce_to_fundamental_domain(std::complex<double> tau)
{
    if (tau.imag() <= 0.0) {
        std::ostringstream msg;
        msg << "period_matrix: τ = " << tau.real() << " + " << tau.imag()
            << "i is not in the upper half-plane (Im(τ) must be > 0).";
        throw std::domain_error(msg.str());
    }

    // Iterate at most 200 times (convergence is rapid for well-conditioned τ)
    for (int k = 0; k < 200; ++k) {
        // T step: shift Re(τ) into [−½, ½)
        double re = tau.real();
        long   n  = static_cast<long>(std::floor(re + 0.5));
        tau -= std::complex<double>(static_cast<double>(n), 0.0);

        // S step: if |τ| < 1, apply τ ← −1/τ
        if (std::abs(tau) < 1.0 - 1e-12) {
            tau = -1.0 / tau;
        } else {
            break;
        }
    }
    return tau;
}

// ─────────────────────────────────────────────────────────────────────────────
// is_in_fundamental_domain — check membership in F with tolerance tol.
// ─────────────────────────────────────────────────────────────────────────────
/// `true` iff `τ` lies inside the standard SL(2,ℤ) fundamental domain
/// with tolerance `tol`.
inline bool is_in_fundamental_domain(std::complex<double> tau, double tol = 1e-9)
{
    if (tau.imag() <= 0.0) return false;
    if (std::abs(tau.real()) > 0.5 + tol) return false;
    if (std::abs(tau) < 1.0 - tol) return false;
    return true;
}

// ─────────────────────────────────────────────────────────────────────────────
// compute_period_matrix
//
// Computes the period data from the Euclidean holonomy translations.
// For genus-1 surfaces, also reduces τ to the fundamental domain.
// ─────────────────────────────────────────────────────────────────────────────
/// Compute the period data from the Euclidean holonomy translations.
/// For genus 1, also normalises `τ` when `reduce` is `true` (default)
/// using `normalizeModulus` — the Java-faithful reduction
/// (`DiscreteEllipticUtility.normalizeModulus`), which folds τ into
/// `0 ≤ Re(τ) ≤ ½`, `Im(τ) ≥ 0`, `|τ| ≥ 1` (the extra `Re ≥ 0` fold
/// uses the mirror symmetry `τ ≅ −τ̄`).  This matches the upstream Java
/// output exactly (Finding 6).  For the canonical SL(2,ℤ) domain
/// (`−½ ≤ Re τ < ½`, no mirror fold) call `reduce_to_fundamental_domain`
/// on `pd.tau` instead.
inline PeriodData compute_period_matrix(const HolonomyData& hol, bool reduce = true)
{
    PeriodData pd;
    pd.omega.reserve(hol.translations.size());
    for (auto& t : hol.translations)
        pd.omega.push_back(std::complex<double>(t.x(), t.y()));

    if (pd.omega.size() < 2) return pd;   // need at least 2 generators

    // τ = ω_2 / ω_1  — choose ω_1 such that Im(τ) > 0
    std::complex<double> w1 = pd.omega[0];
    std::complex<double> w2 = pd.omega[1];
    if (std::abs(w1) < 1e-14) return pd;

    std::complex<double> tau = w2 / w1;
    if (tau.imag() < 0.0) {
        tau = std::conj(tau);              // swap orientation
        w1  = std::conj(w1);
        w2  = std::conj(w2);
        pd.omega[0] = w1;
        pd.omega[1] = w2;
    }
    if (tau.imag() < 0.0) return pd;      // degenerate

    if (reduce) {
        // Java-faithful normalisation (Finding 6): folds τ into
        // 0 ≤ Re ≤ ½, Im ≥ 0, |τ| ≥ 1 via DiscreteEllipticUtility.normalizeModulus.
        tau = normalizeModulus(tau);
        pd.in_fundamental_domain = true;
    }
    pd.tau = tau;
    return pd;
}

} // namespace conformallab