ConformalLabpp/code/include/clausen.hpp

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


// 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().
//
// Note on return value: the loop decrements n one extra time after the
// stopping condition is met, so the returned value is one less than the
// last index checked.  Callers pass the result directly to csevl() as
// the term count, which evaluates terms [0, n-1] — this is intentional
// and matches the Java Clausen.inits() / csevl() contract exactly.
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 `Cl₂(x) = −∫₀ˣ log|2 sin(t/2)| dt`,
/// computed via a high-precision Chebyshev expansion.
/// Same as 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) = Cl₂(2x) / 2`.
/// Same as 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(Li₂(z))` for complex `z`.
/// Same as 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