add GTest infrastructure and port first mathematical tests from Java
- Add Google Test via CMake FetchContent (v1.14.0) - Add clausen.hpp: Clausen integral, Lobachevsky function, Im(Li2) - Add hyper_ideal_utility.hpp: generalized and ideal-vertex hyperbolic tetrahedron volume formulas using Eigen for the 4x4 Gram determinant - Port ClausenTest (5 tests) and HyperIdealUtilityTest (8 tests) from Java/JUnit — all 13 pass with same tolerances as the Java originals - Fix pre-existing VIEWER/viewer case mismatch in CMakeLists.txt Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
179
code/include/clausen.hpp
Normal file
179
code/include/clausen.hpp
Normal file
@@ -0,0 +1,179 @@
|
||||
#pragma once
|
||||
|
||||
// Clausen integral, Lobachevsky function, and Im(Li2).
|
||||
// Ported from de.varylab.discreteconformal.functional.Clausen (Java).
|
||||
// Original algorithm: Boris Springborn, Stefan Sechelmann (TU Berlin).
|
||||
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <cstdlib>
|
||||
|
||||
namespace conformallab {
|
||||
|
||||
namespace detail {
|
||||
|
||||
// Evaluate a Chebyshev series at x using n terms.
|
||||
// Corresponds to Java Clausen.csevl().
|
||||
inline double csevl(double x, const double* cs, int n) noexcept {
|
||||
double b2 = 0.0, b1 = 0.0, b0 = 0.0;
|
||||
const double twox = 2.0 * x;
|
||||
while (n-- > 0) {
|
||||
b2 = b1;
|
||||
b1 = b0;
|
||||
b0 = twox * b1 - b2 + cs[n];
|
||||
}
|
||||
return 0.5 * (b0 - b2);
|
||||
}
|
||||
|
||||
// Count Chebyshev terms needed so truncation error <= eta.
|
||||
// Corresponds to Java Clausen.inits().
|
||||
inline int inits(const double* series, int n, double eta) noexcept {
|
||||
double err = 0.0;
|
||||
while (err <= eta && n-- != 0) {
|
||||
err += std::abs(series[n]);
|
||||
}
|
||||
return n;
|
||||
}
|
||||
|
||||
// Chebyshev expansion of `cl(t)/t + log(t)` around Pi/6.
|
||||
inline const double* clpi6() noexcept {
|
||||
static const double a[] = {
|
||||
2 * 1.0057346496467363858,
|
||||
.0076523796971586786263,
|
||||
.0019223823523180480014,
|
||||
.53333368801173950429e-5,
|
||||
.68684944849366102659e-6,
|
||||
.63769755654413855855e-8,
|
||||
.57069363812137970721e-9,
|
||||
.87936343137236194448e-11,
|
||||
.62365831120408524691e-12,
|
||||
.12996625954032513221e-13,
|
||||
.78762044080566097484e-15,
|
||||
.20080243561666612900e-16,
|
||||
.10916495826127475499e-17,
|
||||
.32027217200949691956e-19,
|
||||
};
|
||||
return a;
|
||||
}
|
||||
|
||||
// Chebyshev expansion around Pi/2.
|
||||
inline const double* clpi2() noexcept {
|
||||
static const double a[] = {
|
||||
2 * .017492908851746863924 + 2 * 1.0057346496467363858,
|
||||
.023421240075284860656 + .0076523796971586786263,
|
||||
.0060025281630108248332 + .0019223823523180480014,
|
||||
.000085934211448718844330 + .53333368801173950429e-5,
|
||||
.000012155033501044820317 + .68684944849366102659e-6,
|
||||
.46587486310623464413e-6 + .63769755654413855855e-8,
|
||||
.50732554559130493329e-7 + .57069363812137970721e-9,
|
||||
.28794458754760053792e-8 + .87936343137236194448e-11,
|
||||
.27792370776596244150e-9 + .62365831120408524691e-12,
|
||||
.19340423475636663004e-10 + .12996625954032513221e-13,
|
||||
.17726134256574610202e-11 + .78762044080566097484e-15,
|
||||
.13811355237660945692e-12 + .20080243561666612900e-16,
|
||||
.12433074161771699487e-13 + .10916495826127475499e-17,
|
||||
.10342683357723940535e-14 + .32027217200949691956e-19,
|
||||
.92910354101990447850e-16,
|
||||
.80428334724548559541e-17,
|
||||
.72598441354406482972e-18,
|
||||
.64475701884829384587e-19,
|
||||
.58630185185185185187e-20,
|
||||
};
|
||||
return a;
|
||||
}
|
||||
|
||||
// Chebyshev expansion of `-cl(Pi-t)/(Pi-t) + log(2)` around 5*Pi/6.
|
||||
inline const double* cl5pi6() noexcept {
|
||||
static const double a[] = {
|
||||
2 * .017492908851746863924,
|
||||
.023421240075284860656,
|
||||
.0060025281630108248332,
|
||||
.000085934211448718844330,
|
||||
.000012155033501044820317,
|
||||
.46587486310623464413e-6,
|
||||
.50732554559130493329e-7,
|
||||
.28794458754760053792e-8,
|
||||
.27792370776596244150e-9,
|
||||
.19340423475636663004e-10,
|
||||
.17726134256574610202e-11,
|
||||
.13811355237660945692e-12,
|
||||
.12433074161771699487e-13,
|
||||
.10342683357723940535e-14,
|
||||
.92910354101990447850e-16,
|
||||
.80428334724548559541e-17,
|
||||
.72598441354406482972e-18,
|
||||
.64475701884829384587e-19,
|
||||
.58630185185185185187e-20,
|
||||
};
|
||||
return a;
|
||||
}
|
||||
|
||||
// Machine epsilon for double (IEEE 754).
|
||||
constexpr double kMachEps = 1.1102230246251565e-16;
|
||||
constexpr double kLn2 = 0.693147180559945309417232121458;
|
||||
|
||||
// Number of Chebyshev terms for each expansion (computed like Java static init).
|
||||
inline int nclpi6() noexcept {
|
||||
static const int n = inits(clpi6(), 14, kMachEps / 10.0);
|
||||
return n;
|
||||
}
|
||||
inline int nclpi2() noexcept {
|
||||
static const int n = inits(clpi2(), 19, kMachEps / 10.0);
|
||||
return n;
|
||||
}
|
||||
inline int ncl5pi6() noexcept {
|
||||
static const int n = inits(cl5pi6(), 19, kMachEps / 10.0);
|
||||
return n;
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
|
||||
// Clausen's integral Cl2(x) = -integral_0^x log|2 sin(t/2)| dt.
|
||||
// High-precision Chebyshev implementation.
|
||||
// Corresponds to Java Clausen.clausen2().
|
||||
inline double clausen2(double x) noexcept {
|
||||
using namespace detail;
|
||||
constexpr double pi = 3.14159265358979323846264338328;
|
||||
constexpr double two_pi = 2.0 * pi;
|
||||
|
||||
int rh = 0;
|
||||
x = std::fmod(x, two_pi);
|
||||
if (x < 0.0) x += two_pi;
|
||||
if (x > pi) { rh = 1; x = two_pi - x; }
|
||||
|
||||
double f;
|
||||
if (x == 0.0) {
|
||||
f = 0.0;
|
||||
} else if (x <= pi / 3.0) {
|
||||
f = csevl(x * (6.0 / pi) - 1.0, clpi6(), nclpi6()) * x - x * std::log(x);
|
||||
} else if (x <= 2.0 * pi / 3.0) {
|
||||
f = csevl(x * (3.0 / pi) - 1.0, clpi2(), nclpi2()) * x - x * std::log(x);
|
||||
} else {
|
||||
f = (kLn2 - csevl(5.0 - x * (6.0 / pi), cl5pi6(), ncl5pi6())) * (pi - x);
|
||||
}
|
||||
return rh ? -f : f;
|
||||
}
|
||||
|
||||
// Milnor's Lobachevsky function Л(x) = Cl2(2x)/2.
|
||||
// Corresponds to Java Clausen.Л().
|
||||
inline double Lobachevsky(double x) noexcept {
|
||||
constexpr double pi = 3.14159265358979323846264338328;
|
||||
x = std::fmod(x, pi);
|
||||
if (x <= 0.0) x += pi;
|
||||
return clausen2(2.0 * x) / 2.0;
|
||||
}
|
||||
|
||||
// Imaginary part of the dilogarithm Im(Li2(z)).
|
||||
// Corresponds to Java Clausen.ImLi2().
|
||||
inline double ImLi2(std::complex<double> z) noexcept {
|
||||
auto a = std::log(1.0 - std::conj(z)); // log(1 - conj(z))
|
||||
auto b = std::log(1.0 - z); // log(1 - z)
|
||||
double x = std::log(std::abs(z));
|
||||
double y = 0.5 * (a - b).imag();
|
||||
double phi = std::arg(z);
|
||||
return y * x + 0.5 * (clausen2(2.0 * y)
|
||||
- clausen2(2.0 * (y + phi))
|
||||
+ clausen2(2.0 * phi));
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
108
code/include/hyper_ideal_utility.hpp
Normal file
108
code/include/hyper_ideal_utility.hpp
Normal file
@@ -0,0 +1,108 @@
|
||||
#pragma once
|
||||
|
||||
// Hyperbolic tetrahedron volume formulas.
|
||||
// Ported from de.varylab.discreteconformal.functional.HyperIdealUtility (Java).
|
||||
|
||||
#include "clausen.hpp"
|
||||
|
||||
#include <Eigen/Dense>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
|
||||
namespace conformallab {
|
||||
|
||||
// Volume of a generalized hyperbolic tetrahedron with dihedral angles A..F.
|
||||
// Formula: Meyerhoff / Ushijima (Springer 2006).
|
||||
// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolume().
|
||||
inline double calculateTetrahedronVolume(double A, double B, double C,
|
||||
double D, double E, double F) {
|
||||
constexpr double pi = 3.14159265358979323846264338328;
|
||||
|
||||
// Degenerate if any angle equals pi.
|
||||
if (A == pi || B == pi || C == pi || D == pi || E == pi || F == pi)
|
||||
return 0.0;
|
||||
|
||||
const double sA = std::sin(A), sB = std::sin(B), sC = std::sin(C);
|
||||
const double sD = std::sin(D), sE = std::sin(E), sF = std::sin(F);
|
||||
const double cA = std::cos(A), cB = std::cos(B), cC = std::cos(C);
|
||||
const double cD = std::cos(D), cE = std::cos(E), cF = std::cos(F);
|
||||
|
||||
// Unit complex numbers e^(i*angle).
|
||||
using Cx = std::complex<double>;
|
||||
auto polar = [](double angle) { return std::polar(1.0, angle); };
|
||||
|
||||
Cx ad = polar(A + D), be = polar(B + E), cf = polar(C + F);
|
||||
Cx abc = polar(A + B + C), abf = polar(A + B + F);
|
||||
Cx ace = polar(A + C + E), aef = polar(A + E + F);
|
||||
Cx bcd = polar(B + C + D), bdf = polar(B + D + F);
|
||||
Cx def = polar(D + E + F), cde = polar(C + D + E);
|
||||
Cx abde = ad * be, acdf = ad * cf, bcef = be * cf;
|
||||
Cx abcdef = abc * def;
|
||||
|
||||
Cx z = ad + be + cf + abf + ace + bcd + def + abcdef;
|
||||
|
||||
// Gram matrix of the tetrahedron.
|
||||
Eigen::Matrix4d G;
|
||||
G << 1.0, -cA, -cB, -cF,
|
||||
-cA, 1.0, -cC, -cE,
|
||||
-cB, -cC, 1.0, -cD,
|
||||
-cF, -cE, -cD, 1.0;
|
||||
|
||||
Cx sqrtG = std::sqrt(Cx(G.determinant(), 0.0));
|
||||
|
||||
Cx f = Cx(sA*sD + sB*sE + sC*sF, 0.0);
|
||||
Cx f1 = f - sqrtG;
|
||||
Cx f2 = f + sqrtG;
|
||||
Cx z1 = -2.0 * f1 / z;
|
||||
Cx z2 = -2.0 * f2 / z;
|
||||
|
||||
auto U = [&](Cx zi) {
|
||||
return 0.5 * (
|
||||
+ ImLi2(zi)
|
||||
+ ImLi2(abde * zi)
|
||||
+ ImLi2(acdf * zi)
|
||||
+ ImLi2(bcef * zi)
|
||||
- ImLi2(-abc * zi)
|
||||
- ImLi2(-aef * zi)
|
||||
- ImLi2(-bdf * zi)
|
||||
- ImLi2(-cde * zi)
|
||||
);
|
||||
};
|
||||
|
||||
return (U(z1) - U(z2)) / 2.0;
|
||||
}
|
||||
|
||||
// Volume of a hyperideal tetrahedron with one ideal vertex (at gamma).
|
||||
// Dihedral angles at the ideal vertex: gamma1, gamma2, gamma3.
|
||||
// Dihedral angles at opposite edges: alpha23, alpha31, alpha12.
|
||||
// Formula: Kolpakov–Mednykh (arxiv math/0603097).
|
||||
// Corresponds to Java HyperIdealUtility.calculateTetrahedronVolumeWithIdealVertexAtGamma().
|
||||
inline double calculateTetrahedronVolumeWithIdealVertexAtGamma(
|
||||
double gamma1, double gamma2, double gamma3,
|
||||
double alpha23, double alpha31, double alpha12)
|
||||
{
|
||||
constexpr double pi = 3.14159265358979323846264338328;
|
||||
auto L = [](double x) { return Lobachevsky(x); };
|
||||
|
||||
double result = L(gamma1) + L(gamma2) + L(gamma3);
|
||||
|
||||
result += L((pi + alpha31 - alpha12 - gamma1) / 2.0);
|
||||
result += L((pi + alpha12 - alpha23 - gamma2) / 2.0);
|
||||
result += L((pi + alpha23 - alpha31 - gamma3) / 2.0);
|
||||
|
||||
result += L((pi - alpha31 + alpha12 - gamma1) / 2.0);
|
||||
result += L((pi - alpha12 + alpha23 - gamma2) / 2.0);
|
||||
result += L((pi - alpha23 + alpha31 - gamma3) / 2.0);
|
||||
|
||||
result += L((pi + alpha31 + alpha12 - gamma1) / 2.0);
|
||||
result += L((pi + alpha12 + alpha23 - gamma2) / 2.0);
|
||||
result += L((pi + alpha23 + alpha31 - gamma3) / 2.0);
|
||||
|
||||
result += L((pi - alpha31 - alpha12 - gamma1) / 2.0);
|
||||
result += L((pi - alpha12 - alpha23 - gamma2) / 2.0);
|
||||
result += L((pi - alpha23 - alpha31 - gamma3) / 2.0);
|
||||
|
||||
return result / 2.0;
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
Reference in New Issue
Block a user