// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// coupling.hpp
// Angular momentum coupling coefficients (Wigner 3-j, 6-j, 9-j symbols)
//
// Provided functions:
//   wigner3j(j1,j2,j3, m1,m2,m3) — Wigner 3-j symbol
//   wigner6j(j1,j2,j3, j4,j5,j6) — Wigner 6-j symbol
//   wigner9j(j1,...,j9)           — Wigner 9-j symbol
//   clebschGordan(j1,m1, j2,m2, J,M) — Clebsch-Gordan coefficient
//
// j values are half-integers (managed as twice the int: two_j = 2j)
// Supported type: double (for exact computation, Rational/Int is recommended)

#ifndef SANGI_SPECIAL_COUPLING_HPP
#define SANGI_SPECIAL_COUPLING_HPP

#include <math/special/special_common.hpp>
#include <cmath>
#include <limits>
#include <algorithm>
#include <vector>

namespace sangi {
namespace special {

// ================================================================
// Internal utility: log factorial
// ================================================================

namespace detail {

inline double log_factorial(int n) {
    if (n < 0) return std::numeric_limits<double>::quiet_NaN();
    if (n <= 1) return 0.0;
    return std::lgamma(static_cast<double>(n + 1));
}

// Triangle condition: |j1-j2| <= j3 <= j1+j2 (tested in 2j units)
inline bool triangle(int two_j1, int two_j2, int two_j3) {
    if ((two_j1 + two_j2 + two_j3) % 2 != 0) return false;
    if (two_j3 < std::abs(two_j1 - two_j2)) return false;
    if (two_j3 > two_j1 + two_j2) return false;
    return true;
}

// Δ(abc) = sqrt[(a+b-c)!(a-b+c)!(-a+b+c)! / (a+b+c+1)!]  (2j input)
inline double triangle_coeff(int two_a, int two_b, int two_c) {
    int s = (two_a + two_b + two_c) / 2;
    int n1 = s - two_c;  // (a+b-c) in 2j units / 2? No.
    // In half-integer: a+b-c = (two_a+two_b-two_c)/2
    int abc = (two_a + two_b - two_c) / 2;
    int acb = (two_a - two_b + two_c) / 2;
    int bca = (-two_a + two_b + two_c) / 2;
    int sum = (two_a + two_b + two_c) / 2 + 1;

    return log_factorial(abc) + log_factorial(acb) + log_factorial(bca)
           - log_factorial(sum);
}

} // namespace detail

// ================================================================
// Wigner 3-j symbol
// (j1 j2 j3)     (-1)^{j1-j2-m3}
// (m1 m2 m3)  =  ───────────────── · CG(j1,m1,j2,m2; j3,-m3)
//                   sqrt(2j3+1)
//
// Input: two_j1, two_j2, two_j3, two_m1, two_m2, two_m3 (all in 2j, 2m units)
// ================================================================

inline double wigner3j(int two_j1, int two_j2, int two_j3,
                       int two_m1, int two_m2, int two_m3) {
    // Selection-rule check
    if (two_m1 + two_m2 + two_m3 != 0) return 0.0;
    if (std::abs(two_m1) > two_j1) return 0.0;
    if (std::abs(two_m2) > two_j2) return 0.0;
    if (std::abs(two_m3) > two_j3) return 0.0;
    if (!detail::triangle(two_j1, two_j2, two_j3)) return 0.0;

    // j+m must be integer → two_j + two_m must be even
    if ((two_j1 + two_m1) % 2 != 0) return 0.0;
    if ((two_j2 + two_m2) % 2 != 0) return 0.0;
    if ((two_j3 + two_m3) % 2 != 0) return 0.0;

    // Racah formula
    double log_prefactor = 0.5 * detail::triangle_coeff(two_j1, two_j2, two_j3);

    // sqrt factor from factorials of (j+m) and (j-m)
    log_prefactor += 0.5 * (
        detail::log_factorial((two_j1 + two_m1) / 2) +
        detail::log_factorial((two_j1 - two_m1) / 2) +
        detail::log_factorial((two_j2 + two_m2) / 2) +
        detail::log_factorial((two_j2 - two_m2) / 2) +
        detail::log_factorial((two_j3 + two_m3) / 2) +
        detail::log_factorial((two_j3 - two_m3) / 2)
    );

    // Sum over t
    int t_min = std::max({0,
        (two_j2 - two_j3 - two_m1) / 2,
        (two_j1 + two_m2 - two_j3) / 2});
    // Wait, need to be more careful. Use half-integer values.
    // t ranges where all factorial arguments are non-negative
    int a1 = (two_j1 + two_j2 - two_j3) / 2;    // j1+j2-j3
    int a2 = (two_j1 - two_m1) / 2;              // j1-m1
    int a3 = (two_j2 + two_m2) / 2;              // j2+m2
    int c1 = (-two_j2 + two_j3 + two_m1) / 2;   // j3-j2+m1
    int c2 = (-two_j1 + two_j3 - two_m2) / 2;   // j3-j1-m2

    t_min = std::max({0, -c1, -c2});
    int t_max = std::min({a1, a2, a3});

    double sum = 0.0;
    for (int t = t_min; t <= t_max; ++t) {
        double log_term = detail::log_factorial(t) +
                          detail::log_factorial(a1 - t) +
                          detail::log_factorial(a2 - t) +
                          detail::log_factorial(a3 - t) +
                          detail::log_factorial(t + c1) +
                          detail::log_factorial(t + c2);
        double sign = (t % 2 == 0) ? 1.0 : -1.0;
        sum += sign * std::exp(log_prefactor - log_term);
    }

    // Phase: (-1)^{(j1-j2-m3)/1} but in 2j units: (-1)^{(two_j1-two_j2-two_m3)/2}
    int phase_exp = (two_j1 - two_j2 - two_m3) / 2;
    double phase = (phase_exp % 2 == 0) ? 1.0 : -1.0;

    return phase * sum;
}

// ================================================================
// Clebsch-Gordan coefficient
// <j1 m1; j2 m2 | J M> = (-1)^{j1-j2+M} sqrt(2J+1) · (j1 j2  J )
//                                                       (m1 m2 -M)
// ================================================================

inline double clebschGordan(int two_j1, int two_m1,
                            int two_j2, int two_m2,
                            int two_J, int two_M) {
    if (two_m1 + two_m2 != two_M) return 0.0;

    int phase_exp = (two_j1 - two_j2 + two_M) / 2;
    double phase = (phase_exp % 2 == 0) ? 1.0 : -1.0;
    double sqrt_2J1 = std::sqrt(static_cast<double>(two_J + 1));

    return phase * sqrt_2J1 * wigner3j(two_j1, two_j2, two_J,
                                        two_m1, two_m2, -two_M);
}

// ================================================================
// Wigner 6-j symbol {j1 j2 j3}
//                  {j4 j5 j6}
// Racah formula
// ================================================================

inline double wigner6j(int two_j1, int two_j2, int two_j3,
                       int two_j4, int two_j5, int two_j6) {
    // Triangle conditions (4 triads)
    if (!detail::triangle(two_j1, two_j2, two_j3)) return 0.0;
    if (!detail::triangle(two_j1, two_j5, two_j6)) return 0.0;
    if (!detail::triangle(two_j4, two_j2, two_j6)) return 0.0;
    if (!detail::triangle(two_j4, two_j5, two_j3)) return 0.0;

    double log_delta = detail::triangle_coeff(two_j1, two_j2, two_j3) +
                       detail::triangle_coeff(two_j1, two_j5, two_j6) +
                       detail::triangle_coeff(two_j4, two_j2, two_j6) +
                       detail::triangle_coeff(two_j4, two_j5, two_j3);
    log_delta *= 0.5;

    // Range of t
    int t_min = std::max({
        (two_j1 + two_j2 + two_j3) / 2,
        (two_j1 + two_j5 + two_j6) / 2,
        (two_j4 + two_j2 + two_j6) / 2,
        (two_j4 + two_j5 + two_j3) / 2
    });
    int t_max = std::min({
        (two_j1 + two_j2 + two_j4 + two_j5) / 2,
        (two_j2 + two_j3 + two_j5 + two_j6) / 2,
        (two_j1 + two_j3 + two_j4 + two_j6) / 2
    });

    double sum = 0.0;
    for (int t = t_min; t <= t_max; ++t) {
        double log_num = detail::log_factorial(t + 1);
        double log_den =
            detail::log_factorial(t - (two_j1 + two_j2 + two_j3) / 2) +
            detail::log_factorial(t - (two_j1 + two_j5 + two_j6) / 2) +
            detail::log_factorial(t - (two_j4 + two_j2 + two_j6) / 2) +
            detail::log_factorial(t - (two_j4 + two_j5 + two_j3) / 2) +
            detail::log_factorial((two_j1 + two_j2 + two_j4 + two_j5) / 2 - t) +
            detail::log_factorial((two_j2 + two_j3 + two_j5 + two_j6) / 2 - t) +
            detail::log_factorial((two_j1 + two_j3 + two_j4 + two_j6) / 2 - t);
        double sign = (t % 2 == 0) ? 1.0 : -1.0;
        sum += sign * std::exp(log_delta + log_num - log_den);
    }
    return sum;
}

// ================================================================
// Wigner 9-j symbol: expressed as a sum of products of 3-j symbols
// {j1 j2 j3}
// {j4 j5 j6} = Σ_t (2t+1) {j1 j4 j7} {j2 j5 j8} {j3 j6 j9}
// {j7 j8 j9}                {j8 t  j6}  {j4 t  j3}  {j1 t  j5}
//                                         ... (product of Wigner 6-j symbols)
// ================================================================

inline double wigner9j(int two_j1, int two_j2, int two_j3,
                       int two_j4, int two_j5, int two_j6,
                       int two_j7, int two_j8, int two_j9) {
    // 9-j = Σ_t (-1)^{2t} (2t+1) {j1 j4 j7} {j2 j5 j8} {j3 j6 j9}
    //                               {j8 j9 t}  {j4 j6 t}  {j2 j7 t}
    // Wait, let me use the standard formula:
    // {j1 j2 j3}
    // {j4 j5 j6} = Σ_t (-1)^{2t} (2t+1)
    // {j7 j8 j9}     × {j1 j4 j7} {j2 j5 j8} {j3 j6 j9}
    //                   {j8 j9 t } {j4 j6 t } {j7 j2 t }
    // No, the standard definition:
    // 9j = Σ_{two_t} (-1)^{two_t} (two_t+1)
    //      × 6j(j1,j4,j7, j8,j9,t) × 6j(j2,j5,j8, j4,t,j6) × 6j(j3,j6,j9, t,j1,j2)
    // Wait, I need to look up the correct formula more carefully.
    // Standard: 9j{a b c; d e f; g h j} =
    //   Σ_x (-1)^{2x}(2x+1) 6j{a b c; f j x} 6j{d e f; b x h} 6j{g h j; x a d}

    int two_t_min = std::max({
        std::abs(two_j1 - two_j9),
        std::abs(two_j2 - two_j6),
        std::abs(two_j4 - two_j8)
    });
    int two_t_max = std::min({
        two_j1 + two_j9,
        two_j2 + two_j6,
        two_j4 + two_j8
    });

    double sum = 0.0;
    for (int two_t = two_t_min; two_t <= two_t_max; two_t += 2) {
        double sign = (two_t % 2 == 0) ? 1.0 : -1.0;  // (-1)^{2t} = (-1)^{two_t}
        // Actually (-1)^{2t} where t = two_t/2, so 2t = two_t. (-1)^{two_t}.
        double w6j_1 = wigner6j(two_j1, two_j2, two_j3,
                                two_j6, two_j9, two_t);
        double w6j_2 = wigner6j(two_j4, two_j5, two_j6,
                                two_j2, two_t, two_j8);
        double w6j_3 = wigner6j(two_j7, two_j8, two_j9,
                                two_t, two_j1, two_j4);

        sum += sign * static_cast<double>(two_t + 1) * w6j_1 * w6j_2 * w6j_3;
    }
    return sum;
}

} // namespace special
} // namespace sangi

#endif // SANGI_SPECIAL_COUPLING_HPP