<repo_name>greenwald/RelHel<file_sep>/src/Contraction.cxx
#include "Contraction.h"

#include "ClebschGordan.h"

#include <algorithm>

#include <iostream>

namespace relhel {

//-------------------------
namespace contractions {
const std::vector<Contraction::Term::index> psi1_psi2 = {Contraction::Term::index::psi1, Contraction::Term::index::psi2};
const std::vector<Contraction::Term::index> psi1_chi  = {Contraction::Term::index::psi1, Contraction::Term::index::chi};
const std::vector<Contraction::Term::index> psi2_chi  = {Contraction::Term::index::psi2, Contraction::Term::index::chi};
const std::vector<Contraction::Term::index> psi1_phi  = {Contraction::Term::index::psi1, Contraction::Term::index::phi};
const std::vector<Contraction::Term::index> psi2_phi  = {Contraction::Term::index::psi2, Contraction::Term::index::phi};
const std::vector<Contraction::Term::index> chi_phi   = {Contraction::Term::index::chi,  Contraction::Term::index::phi};
}

//-------------------------
Contraction::Contraction(const CoupledWaveFunctions& psi, const OrbitalAngularMomentumWaveFunction& chi, const WaveFunction& phi)
{
    // if chi has no projection = zero wave functions
    if (chi.projections().at(0).empty())
        return;
    
    // loop over spin projections of first daughter
    for (const auto& m_wps : psi.phi()[0].projections()) {
        if (m_wps.second.empty())
            continue;

        // find appropriate spin projection of second daughter
        auto it = psi.phi()[1].projections().find(psi.delta() - m_wps.first);
        if (it == psi.phi()[1].projections().end() or it->second.empty())
            continue;

        // loop over WaveProduct's of first daughter
        for (const auto& wp1 : m_wps.second) {
            auto two_j1 = spin(wp1);
            auto wp1_coeff = squared_coefficient(wp1);
            
            // loop over WaveProduct's of second daughter
            for (const auto& wp2 : it->second) {
                auto two_j2 = spin(wp2);
                auto psi_coeff = wp1_coeff * squared_coefficient(wp2)
                    * ClebschGordan::squared_coefficient(two_j1, m_wps.first, two_j2, it->first, psi.twoS());
                
                // loop over WaveProduct's of orbital angular momentum
                for (const auto& wp_chi : chi.projections().at(0)) {
                    auto l = spin(wp_chi) / 2;

                    auto psi_chi_coeff = psi_coeff * squared_coefficient(wp_chi) * pow(decompose(2), l);
                                
                    // loop over projections of initial state
                    for (const auto& m_wps_phi : phi.projections())
                        // loop over WaveProduct's of initial state
                        for (const auto& wp_phi : m_wps_phi.second)
                    
                            Terms_.emplace_back(wp1, wp2, wp_chi, wp_phi, squared_coefficient(wp_phi) * psi_chi_coeff);
                                 
                }
            }
        }
    }
    prune();
}

//-------------------------
const std::vector<int> gamma_exponents(const Contraction::Term& t, const std::vector<Contraction::Term::index>& I)
{
     std::vector<int> g;
     g.reserve(I.size());
     std::transform(I.begin(), I.end(), std::back_inserter(g), [&t](Contraction::Term::index i){return t.GammaExponents[static_cast<size_t>(i)]; });
     return g;
}

//-------------------------
// helper function
GammaPolynomial& add_term(GammaPolynomial& P, const GammaPolynomial::key_type& ges, const GammaPolynomial::mapped_type::value_type& coeff)
{
    // look for term with same exponents
    auto it = P.find(ges);
    // if found
    if (it != P.end()) {
        // look for coefficient term with same argument under the square root
        auto it2 = std::find_if(it->second.begin(), it->second.end(),
                                [&coeff](const GammaPolynomial::mapped_type::value_type& r)
                                {return r.second == coeff.second;});
        // if one is found, add arguments outside of square root sign
        if (it2 != it->second.end()) {
            it2->first += coeff.first;
            if (it2->first == 0)
                it->second.erase(it2);
        }
        // else add new term completely
        else it->second.push_back(coeff);
    }
    // else add new term
    else {
        P[ges].push_back(coeff);
    }
    return P;
}

//-------------------------
// helper function
GammaPolynomial& add_term(GammaPolynomial& P, const GammaPolynomial::key_type& ges, const GammaPolynomial::mapped_type& coeffs)
{
    std::for_each(coeffs.begin(), coeffs.end(), [&](const GammaPolynomial::mapped_type::value_type& c){add_term(P, ges, c);});
    return P;
}

//-------------------------
// helper function
GammaPolynomial& add_term(GammaPolynomial& P, const GammaPolynomial::value_type& ges_coeffs)
{ return add_term(P, ges_coeffs.first, ges_coeffs.second); }

//-------------------------
const GammaPolynomial gamma_polynomial(const std::vector<Contraction::Term>& terms)
{
    if (std::any_of(terms.begin(), terms.end(), [](const Contraction::Term& t){return rank(t) > 0;}))
        throw std::invalid_argument("Terms are not scalar; gamma_polynomial(...)");
    
    GammaPolynomial P;

    for (const auto& t : terms) {
        auto c = factorize_sqrt(t.SquaredCoefficient);
        add_term(P, gamma_exponents(t, contractions::psi1_psi2), std::make_pair(static_cast<double>(sqrt(c[0])), c[1]));
    }
    
    return P;
}

//-------------------------
// helper function
std::string to_exp_string(const GammaPolynomial::key_type& ges)
{
    std::string S;
    for (size_t i = 0; i < ges.size(); ++i)
        S += exponential_string(" * g_" + std::to_string(i + 1), ges[i]);
    return (S.empty()) ? S : S.erase(0, 3);
}

//-------------------------
// helper function
std::string to_coef_string(const GammaPolynomial::mapped_type& coefs)
{
    return "("
        + std::accumulate(coefs.begin(), coefs.end(), std::string(""),
                          [](std::string& S, const GammaPolynomial::mapped_type::value_type& r)
                          {return S += " + " + to_sqrt_string(RationalNumber(r.first * std::abs(r.first)) * r.second);}).erase(0, 3)
        + ")";
}

//-------------------------
GammaPolynomial& operator+=(GammaPolynomial& lhs, const GammaPolynomial& rhs)
{
    for (const auto& term : rhs)
        add_term(lhs, term);
    return lhs;
}

//-------------------------
std::string to_string(const GammaPolynomial::value_type& ges_coefs)
{
    auto s_ge = to_exp_string(ges_coefs.first);
    return s_ge + (!s_ge.empty() ? " * " : "") + to_coef_string(ges_coefs.second);
}

//-------------------------
std::string to_string(const GammaPolynomial& gp)
{
    
    return std::accumulate(gp.begin(), gp.end(), std::string(),
                           [](std::string& s, const GammaPolynomial::value_type& t){return s += " + " + to_string(t);}).erase(0, 3);
}

//-------------------------
std::string to_string(const Contraction::Term& term)
{
    return std::accumulate(term.WaveProducts.begin(), term.WaveProducts.end(), std::string(),
                           [](std::string& s, const WaveProduct& w){return s += ", " + to_string(w);}).erase(0, 2)
        + " : Gamma(" + std::accumulate(term.GammaExponents.begin(), term.GammaExponents.end(), std::string(),
                                        [](std::string& s, int g){return s += ", " + std::to_string(g);}).erase(0, 2) + ")"
        + " * " + to_sqrt_string(term.SquaredCoefficient);
}

//-------------------------
std::string to_string(const std::vector<Contraction::Term>& terms)
{
    return std::accumulate(terms.begin(), terms.end(), std::string(),
                           [](std::string& s, const Contraction::Term& t){return s += "\n" + to_string(t);}).erase(0, 1);
}

//-------------------------
std::string to_string(Contraction::Term::index i)
{
    switch(i) {
    case Contraction::Term::index::psi1:
        return "psi1";
    case Contraction::Term::index::psi2:
        return "psi2";
    case Contraction::Term::index::chi:
        return "chi";
    case Contraction::Term::index::phi:
        return "phi";
    default:
        return std::to_string(static_cast<size_t>(i));
    }
}

//-------------------------
const std::vector<Contraction::Term::index> triple_contraction(const Contraction::Term& t)
{
    std::vector<Contraction::Term::index> I;
    if (rank(t) != 3 or std::count_if(t.WaveProducts.begin(), t.WaveProducts.end(), [](const WaveProduct& w){return rank(w) == 1;}) != 3)
        return I;
    I = {Contraction::Term::index::psi1, Contraction::Term::index::psi2, Contraction::Term::index::chi, Contraction::Term::index::phi};
    I.erase(I.begin() + std::distance(t.WaveProducts.begin(), std::find_if(t.WaveProducts.begin(), t.WaveProducts.end(), [](const WaveProduct& w){return rank(w) == 0;})));
    return I;
}

//-------------------------
std::string to_string(const std::vector<Contraction::Term::index>& I)
{ return std::accumulate(I.begin(), I.end(), std::string(""), [](std::string& s, Contraction::Term::index i){return s += ", " + to_string(i);}).erase(0, 2); }

//-------------------------
void Contraction::Term::contract(std::vector<index> I)
{
    // check size
    if (I.size() < 2 or I.size() > 3)
        throw std::invalid_argument("Invalid number of contraction indices; contract(" + to_string(I) + ")");

    // check no WaveFunction is empty
    if (std::any_of(I.begin(), I.end(), [&](index i){return WaveProducts[static_cast<size_t>(i)].empty();})) {
        std::vector<index> J;
        J.reserve(I.size());
        std::copy_if(I.begin(), I.end(), std::back_inserter(J), [&](index i){return WaveProducts[static_cast<size_t>(i)].empty();});
        throw std::invalid_argument("Rank-zero Wave function(s) " + to_string(J) + "; contract(" + to_string(I) + ")");
    }
    
    // sort I (aids checking in the 3-index situation below)
    std::sort(I.begin(), I.end(), [](index a, index b){return std::less<size_t>()(static_cast<size_t>(a), static_cast<size_t>(b));});
    
    // check no index is repeated
    if (std::adjacent_find(I.begin(), I.end()) != I.end())
        throw std::invalid_argument("Can't contract WaveFunction with itself; contract(" + to_string(I) + ")");

    // get back elements and pop them off
    std::vector<int> m;
    m.reserve(I.size());
    for (auto i : I) {
        m.push_back(WaveProducts[static_cast<size_t>(i)].back());
        WaveProducts[static_cast<size_t>(i)].pop_back();
    }

    // two WaveFunction's contracted
    if (I.size() == 2) {

        // if contracting with phi
        if (I[1] == index::phi) {
            // if not the same 
            if (m[0] != m[1])
                SquaredCoefficient = RationalNumber(0);
        }
        else {
            // if opposite signs
            if (m[0] * m[1] < 0)
                SquaredCoefficient.negate();
            
            // if contracting with chi, and not opposite signed (or both zero)
            if (I[1] == index::chi and m[0] != -m[1])
                SquaredCoefficient = RationalNumber(0);
        }
    }
    // three WaveFunction's contracted
    else {

        // if not one and only one spin projection is zero
        if (std::count(m.begin(), m.end(), 0) != 1)
            SquaredCoefficient = RationalNumber(0);
        // or the first two projections are equal
        else if (m[0] == m[1])
            SquaredCoefficient = RationalNumber(0);
        else if (m[2] != 0) {
            // if contracting with initial state
            if (I[2] == index::phi) {
                if (m[0] + m[1] + m[2] == 0)
                    SquaredCoefficient = RationalNumber(0);
            } else {
                if (m[2] != 0 and m[0] + m[1] + m[2] != 0)
                    SquaredCoefficient = RationalNumber(0);
            }
        }

        // if not yet zero'ed, check negation situations
        if (!is_zero(SquaredCoefficient)) {

            // if m[2] == 0, m[0] dictates sign
            if (m[2] == 0 and m[0] < 0)
                SquaredCoefficient.negate();
            
            // if m[1] == 0, m[2] dictates sign
            if (m[1] == 0 and m[2] < 0)
                SquaredCoefficient.negate();
            
            // if m[0] == 0, m[2] dictates sign
            if (m[0] == 0 and m[2] > 0)
                SquaredCoefficient.negate();
        }
    }

    // increase relevant gamma factors
    for (size_t i = 0; i < I.size(); ++i)
        if (m[i] == 0) ++GammaExponents[static_cast<size_t>(I[i])];
}

}
<file_sep>/CMakeLists.txt
cmake_minimum_required(VERSION 2.8.12)
project(RelHel)

set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmakeModules)

# Default C++ flags
include(CheckCXXCompilerFlag)
CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11)
if (COMPILER_SUPPORTS_CXX11)
  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
else()
	message(FATAL_ERROR "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
endif()

CHECK_CXX_COMPILER_FLAG("-fdiagnostics-color=always" COMPILER_SUPPORTS_DIAGNOSTIC_COLORS)
if (COMPILER_SUPPORTS_DIAGNOSTIC_COLORS)
  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fdiagnostics-color=always")
endif()

set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -pedantic -Wno-long-long -Wno-missing-braces -Werror=overloaded-virtual -pthread" )

# Default build type is 'Release' (cmake default one is '')
if(NOT CMAKE_BUILD_TYPE)
	set(CMAKE_BUILD_TYPE Release CACHE STRING
		"Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel."
		FORCE)
endif(NOT CMAKE_BUILD_TYPE)

# build with cmake -DCMAKE_BUILD_TYPE=Debug
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -g2" )

# build with cmake -DCMAKE_BUILD_TYPE=Release
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3 -D ELPP_DISABLE_DEBUG_LOGS" )

# Add subdirectories
add_subdirectory(src)
add_subdirectory(programs)

<file_sep>/src/MathUtils.cxx
#include "MathUtils.h"

#include <algorithm>
#include <cstdlib>
#include <functional>

namespace relhel {

//-------------------------
const std::vector<unsigned> triangle(unsigned two_s1, unsigned two_s2)
{
    unsigned two_smin = std::abs(static_cast<int>(two_s1) - static_cast<int>(two_s2));
    unsigned two_smax = two_s1 + two_s2;
    std::vector<unsigned> two_S;
    two_S.reserve((two_smax - two_smin) / 2 + 1);
    for (auto two_s = two_smin; two_s <= two_smax; two_s += 2)
        two_S.push_back(two_s);
    return two_S;
}

//-------------------------
const std::vector<int> projections(unsigned two_j)
{
    std::vector<int> two_M;
    for (int two_m = -two_j; two_m <= static_cast<int>(two_j); ++two_m)
        two_M.push_back(two_m);
    return two_M;
}
    
//-------------------------
const std::vector<std::vector<int> > projections(const std::vector<unsigned>& two_J)
{
    // initialize vector of spin projections to -two_j
    std::vector<int> two_M;
    two_M.reserve(two_J.size());
    std::transform(two_J.begin(), two_J.end(), std::back_inserter(two_M), std::negate<int>());
    
    std::vector<std::vector<int> > SPV;
    // fill SPV with "odometer"-style looping
    while (two_M.back() <= (int)two_J.back()) {
        SPV.push_back(two_M);
        two_M[0] += 2;
        for (size_t i = 0; (i < two_M.size() - 1) and (two_M[i] > (int)two_J[i]); ++i) {
            two_M[i] = -two_J[i];
            two_M[i + 1] += 2;
        }
    }
    return SPV;
}

}
<file_sep>/src/RationalNumber.cxx
#include "RationalNumber.h"

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <numeric>
#include <string>

namespace relhel {

namespace prime_cache {

    // initialize with first primes up to 1000
    static std::vector<unsigned> primes_ = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
                                            59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
                                            127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
                                            191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
                                            257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
                                            331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
                                            401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
                                            467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557,
                                            563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619,
                                            631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
                                            709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
                                            797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
                                            877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953,
                                            967, 971, 977, 983, 991, 997};

    //-------------------------
    bool is_prime(unsigned p)
    {
        if (primes_.empty() and p == 2)
            return true;

        // check for upper limit of search range: p / 2
        auto last = std::upper_bound(primes_.begin(), primes_.end(), 0.5 * p);

        // search if any prime divides p
        for (auto it = primes_.begin(); it != last; ++it)
            if (p % *it == 0)
                return false;

        // if couldn't search high enough range, through
        if (last == primes_.end())
            throw;

        // else p is prime
        return true;
    }

    //-------------------------
    std::vector<unsigned>::iterator next(std::vector<unsigned>::iterator it)
    {
        // increment to next prime
        ++it;
        // if prime available, return it
        if (it != primes_.end())
            return it;
        // else add to cache, brute force
        // initialize p to next value
        unsigned p = primes_.empty() ? 2 : primes_.back() + 1;
        // increase until prime is found
        while (!is_prime(p)) ++p;
        // add to cache
        primes_.push_back(p);
        return --primes_.end();
    }
}

//-------------------------
PrimeFactors::operator unsigned() const
{
    return std::accumulate(Factors_.begin(), Factors_.end(), 1u,
                           [](unsigned a, const map_type::value_type& p_e)
                           { return a * std::pow(p_e.first, p_e.second); });
}

//-------------------------
std::string to_string(const PrimeFactors& pf)
{
    return is_one(pf) ? "1" : 
        std::accumulate(pf.begin(), pf.end(), std::string(),
                        [](std::string& s, const PrimeFactors::map_type::value_type& p_e)
                        { return s += " * " + exponential_string(std::to_string(p_e.first), p_e.second); }
            ).erase(0, 3);
}

//-------------------------
PrimeFactors decompose(const std::vector<unsigned>& N)
{
    PrimeFactors PF;
    for (unsigned n : N) {
        if (n == 0)
            PF[0] = 1;
        else {
            // check if in prime list
            auto pit = std::lower_bound(prime_cache::primes_.begin(), prime_cache::primes_.end(), n);
            if (pit != prime_cache::primes_.end() and *pit == n)
                PF[*pit] = 1;
            else {
                for (auto it = prime_cache::primes_.begin(); n > 1; it = prime_cache::next(it)) {
                    unsigned m = 0;
                    while (n % *it == 0) {
                        ++m;
                        n /= *it;
                    }
                    if (m > 0)
                        PF[*it] = m;
                }
            }
        }
    }
    return PF;
}

//-------------------------
PrimeFactors decompose_factorial(const std::vector<unsigned>& F)
{
    std::vector<unsigned> N;
    N.reserve(std::max(std::accumulate(F.begin(), F.end(), size_t(0)), F.size()) - F.size());
    for (unsigned f : F)
        while (f > 1)
            N.push_back(f--);

    return decompose(N);
}

//-------------------------
PrimeFactors& remove_zeroes(PrimeFactors& PF)
{
    for (auto it = PF.begin(); it != PF.end();)
        if (it->second == 0)
            it = PF.erase(it);
        else
            ++it;
    return PF;
}

//-------------------------
void remove_common(PrimeFactors& A, PrimeFactors& B)
{    
    for (auto& p_e : A) {
        auto it = B.find(p_e.first);
        if (it != B.end()) {
            unsigned c = std::min(p_e.second, it->second);
            p_e.second -= c;
            it->second -= c;
        }
    }
}

//-------------------------
PrimeFactors common(PrimeFactors& A, PrimeFactors& B)
{
    PrimeFactors C;
    for (auto& p_e : A) {
        auto it = B.find(p_e.first);
        if (it != B.end()) {
            unsigned c = std::min(p_e.second, it->second);
            C[p_e.first] = c;
            p_e.second -= c;
            it->second -= c;
        }
    }
    return C;
}

//-------------------------
PrimeFactors sqrt(PrimeFactors pf)
{
    for (auto& p_e : pf)
        if (is_even(p_e.second))
            p_e.second /= 2;
        else
            throw std::runtime_error("sqrt is not integral");
    return pf;
}

//-------------------------
std::array<PrimeFactors, 2> factorize_sqrt(PrimeFactors pf)
{
    PrimeFactors pf2;
    for (auto& p_e : pf)
        if (is_odd(p_e.second)) {
            pf2[p_e.first] = 1;
            --p_e.second;
        }
    return {pf, pf2};
}

//-------------------------
std::string to_string(const std::vector<long>& V)
{
    return std::accumulate(V.begin(), V.end(), std::string(),
                           [](std::string& s, long v)
                           {return s += " " + std::to_string(v);}).erase(0, 1);
}

//-------------------------
double sign_factor(int n, int d)
{
    if (n != 0 and d != 0)
        return sign_of(n * d);
    if (d != 0)
        return 0;
    if (n == 0)
        return std::numeric_limits<double>::quiet_NaN();
    return sign_of(n) * std::numeric_limits<double>::infinity();
}

//-------------------------
void remove_common(RationalNumber& A, RationalNumber& B)
{
    auto nA = A.numerator();
    auto nB = B.numerator();
    auto dA = A.denominator();
    auto dB = B.denominator();
    remove_common(nA, nB);
    remove_common(dA, dB);
    A = RationalNumber(nA, dA, A.sign());
    B = RationalNumber(nB, dB, B.sign());
}

//-------------------------
RationalNumber common(RationalNumber& A, RationalNumber& B)
{
    auto nA = A.numerator();
    auto nB = B.numerator();
    auto dA = A.denominator();
    auto dB = B.denominator();
    auto nC = common(nA, nB);
    auto dC = common(dA, dB);
    A = RationalNumber(nA, dA, A.sign());
    B = RationalNumber(nB, dB, B.sign());
    return RationalNumber(nC, dC);
}

//-------------------------
const bool operator<(const RationalNumber& lhs, const RationalNumber& rhs)
{
    // check by sign first
    if (lhs.sign() < rhs.sign())
        return true;

    if (lhs.sign() > rhs.sign())
        return false;

    // signs must be the same now:
    if (is_zero(lhs)) // therefore if both are zero
        return false;

    // remove common factors and compare uncommon factors
    auto lhs_temp = lhs;
    auto rhs_temp = rhs;
    remove_common(lhs_temp, rhs_temp);
    return static_cast<double>(lhs_temp) < static_cast<double>(rhs_temp);
}

//-------------------------
RationalNumber& operator+=(RationalNumber& lhs, const RationalNumber& rhs)
{
    if (rhs.sign() == 0 or !std::isfinite(lhs.sign()))
        return lhs;
    if (lhs.sign() == 0 or !std::isfinite(rhs.sign()))
        return lhs = rhs;

    // form new numerator
    auto LR = lhs.numerator() * rhs.denominator();
    auto RL = rhs.numerator() * lhs.denominator();
    // remove common factors
    auto C = common(LR, RL);

    // calculate sum of noncommon factors, including signs
    int num = lhs.sign() * static_cast<unsigned>(LR) + rhs.sign() * static_cast<unsigned>(RL);

    return lhs = RationalNumber(C * decompose(std::abs(num)), lhs.denominator() * rhs.denominator(), static_cast<double>(num));
}

//-------------------------
RationalNumber& operator*=(RationalNumber& lhs, const RationalNumber& rhs)
{
    return lhs = RationalNumber(lhs.numerator() * rhs.numerator(),
                          lhs.denominator() * rhs.denominator(),
                          multiply_sign_factors(lhs.sign(), rhs.sign()));
}

//-------------------------
double multiply_sign_factors(double s1, double s2)
{
    if ((s1 == 0 and std::isinf(s2)) or (s2 == 0 and std::isinf(s1)))
        return std::numeric_limits<double>::quiet_NaN();
    return s1 * s2;
}

//-------------------------
double invert_sign_factor(double s)
{
    if (s == 0)
        return std::numeric_limits<double>::infinity();
    if (std::isinf(s))
        return 0;
    return s;
}

//-------------------------
std::string to_string(const RationalNumber& f)
{
    if (f.sign() == 0)
        return "0";
    if (!std::isfinite(f.sign()))
        return std::to_string(f.sign());
    auto N = static_cast<unsigned>(f.numerator());
    auto D = static_cast<unsigned>(f.denominator());
    return (f.sign() > 0 ? "" : "-") + std::to_string(N) + (D == 1 ? "" : " / " + std::to_string(D));
}

std::string to_header_string(const RationalNumber& f)
{
    if (f.sign() == 0)
        return "{0,1}";
    if (std::isinf(f.sign()))
        return "{1,0}";
    if (std::isnan(f.sign()))
        return "{0,0}";
    return "{" + std::to_string(static_cast<unsigned>(f.numerator()) * f.sign()) 
        + "," + std::to_string(static_cast<unsigned>(f.denominator())) + "}";
}

//-------------------------
std::string to_detailed_string(const RationalNumber& f)
{
    return (std::isfinite(f.sign()) ? (f.sign() > 0 ? "" : "-") : (std::to_string(f.sign()) + " * "))
                    + "(" + to_string(f.numerator()) + ") / (" + to_string(f.denominator()) + ")";
}

//-------------------------
std::array<RationalNumber, 2> factorize_sqrt(const RationalNumber& f)
{
    // factorize numerator and denominator
    auto N = factorize_sqrt(f.numerator());
    auto D = factorize_sqrt(f.denominator());
    // and create new rational numbers
    return {RationalNumber(N[0], D[0], f.sign()), RationalNumber(N[1], D[1], 1)};
}

//-------------------------
std::string to_sqrt_string(const std::array<RationalNumber, 2>& f)
{
    if (is_zero(f[0]) or is_zero(f[1]))
        return "0";

    if (is_one(f[1]))
        return to_string(sqrt(f[0]));

    std::string s;
    if (is_one(f[0].numerator()) and is_one(f[0].denominator()))
        s = static_cast<double>(sqrt(f[0])) < 0 ? "-" : "";
    else
        s = to_string(sqrt(f[0])) + " * ";
    
    return s + "sqrt("
        + std::to_string(static_cast<unsigned>(f[1].numerator())) + " / "
        + std::to_string(static_cast<unsigned>(f[1].denominator())) + ")";
}


}
<file_sep>/programs/test.cxx
#include "ClebschGordan.h"
#include "MathUtils.h"
#include "RationalNumber.h"

#include <iostream>
#include <map>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>

using namespace relhel;

//-------------------------
using f_map = std::map<unsigned, RationalNumber>;

//-------------------------
f_map f(unsigned two_j, int two_lambda)
{
    if (is_odd(two_j + two_lambda))
        throw std::invalid_argument("j and lambda must be either both integral or half-integral");
    
    unsigned j_p_lambda = (two_j + two_lambda) / 2;
    unsigned j_m_lambda = (two_j - two_lambda) / 2;
    
    if (abs(two_lambda) > two_j)
        throw std::runtime_error("abs(lambda) > j");

    RationalNumber a(decompose_factorial(j_p_lambda) * decompose_factorial(j_m_lambda),
                     decompose_factorial(two_j));
    
    f_map F;
    for (unsigned n = is_even(j_p_lambda) ? 0 : 1; n <= (two_j - abs(two_lambda)) / 2; n += 2)  {
        F.emplace(n, a * RationalNumber(decompose_factorial(two_j / 2) * pow(decompose(2), n),
                                        decompose_factorial(n)
                                        * decompose_factorial((j_p_lambda - n) / 2)
                                        * decompose_factorial((j_m_lambda - n) / 2)));
    }
    return F;
}

//-------------------------
using ff_map = std::map<std::array<unsigned, 2>, RationalNumber>;

//-------------------------
ff_map ff(unsigned two_j1, int two_lambda1, unsigned two_j2, int two_lambda2)
{
    auto F1 = f(two_j1, two_lambda1);
    auto F2 = f(two_j2, two_lambda2);
    ff_map FF;
    for (const auto& f1 : F1)
        for (const auto& f2 : F2)
            FF.emplace(ff_map::key_type({f1.first, f2.first}), f1.second * f2.second);
    return FF;
}

//-------------------------
const bool is_one(const ff_map& FF)
{ return FF.empty() or (FF.size() == 1 and FF.begin()->first[0] == 0 and FF.begin()->first[1] == 0); }

//-------------------------
std::string to_string(const ff_map::value_type& ff)
{
    return (!is_one(ff.second) or (ff.first[0] == 0 and ff.first[1] == 0) ? to_string(ff.second) : "")
        + (!is_one(ff.second) and (ff.first[0] != 0 or ff.first[1] != 0) ? " * " : "")
        + exponential_string("g_1", ff.first[0])
        + (ff.first[0] != 0 and ff.first[1] != 0 ? " * " : "")
        + exponential_string("g_2", ff.first[1]);
}

//-------------------------
std::string to_string(const ff_map& ff)
{
    return std::accumulate(ff.begin(), ff.end(), std::string(),
                           [](std::string& s, const ff_map::value_type& f)
                           { return s += " + " + to_string(f); }).erase(0, 3);
}

//-------------------------
std::string prefix_string(int j)
{ if (j > 0) return "+"; return (j == 0) ? " " : ""; }



//-------------------------
int main()
{

    unsigned two_J = 0;
    std::vector<unsigned> two_j = {2 , 2};
    
    for (auto two_S : triangle(two_j[0], two_j[1]))
        for (auto two_L : triangle(two_J, two_S)) {
            RationalNumber LJ2(two_L + 1, two_J + 1);
            for (auto two_lambda : projections(two_j)) {

                auto CG2 = ClebschGordan::squared_coefficient(two_L, 0, two_S, two_lambda[0] - two_lambda[1], two_J)
                    * ClebschGordan::squared_coefficient(two_j[0], two_lambda[0], two_j[1], -two_lambda[1], two_S);

                if (is_zero(CG2))
                    continue;
                
                std::cout << spin_to_string(two_J) << " -> "
                          << spin_to_string(two_j[0]) << " [" << prefix_string(two_lambda[0]) << spin_to_string(two_lambda[0]) << "]"
                          << " + "
                          << spin_to_string(two_j[1]) << " [" << prefix_string(two_lambda[1]) << spin_to_string(two_lambda[1]) << "]"
                          << ", L = " << spin_to_string(two_L)
                          << ", S = " << spin_to_string(two_S)
                          << "\t" << std::flush;

                if (is_zero(CG2)) {
                    std::cout << "0" << std::endl;
                    continue;
                }

                auto term = ff(two_j[0], two_lambda[0], two_j[1], two_lambda[1]);

                if (is_one(CG2) and is_one(LJ2)) {
                    std::cout << to_string(term) << std::endl;
                    continue;
                }

                std::cout << to_sqrt_string(LJ2 * CG2) << std::flush;

                if (!is_one(term))
                    std::cout << " * "
                              << (term.size() > 1 ? "(" : "")
                              << to_string(term)
                              << (term.size() > 1 ? ")" : "") << std::flush;

                std::cout << std::endl;
            }
        }
    
    return 0;
}
<file_sep>/README.md
# RelHel

Relativistic spin-amplitude calculator for integer spins.
<file_sep>/include/RationalNumber.h
/// \file
/// \author <NAME>

#ifndef RationalNumber_h
#define RationalNumber_h

#include "MathUtils.h"

#include <array>
#include <cmath>
#include <limits>
#include <map>
#include <string>
#include <vector>

namespace relhel {

/// \class PrimeFactors
/// \brief Prime factor decomposition of an unsigned integer
/// \author <NAME>
class PrimeFactors
{
 public:

    /// maps prime factor to exponent
    using map_type = std::map<unsigned, unsigned>;
    
    /// access operator
    map_type::mapped_type& operator[](const map_type::key_type& key)
    { return Factors_[key]; }

    /// access begin
    map_type::iterator begin()
    { return Factors_.begin(); }

    /// access begin
    map_type::const_iterator begin() const
    { return Factors_.begin(); }

    /// access end
    map_type::iterator end()
    { return Factors_.end(); }

    /// access end
    map_type::const_iterator end() const
    { return Factors_.end(); }

    const bool empty() const
    { return Factors_.empty(); }
    
    /// erase
    map_type::iterator erase(map_type::const_iterator pos)
    { return Factors_.erase(pos); }

    map_type::iterator find(const map_type::key_type& key)
    { return Factors_.find(key); }
    
    /// explicit conversion to unsigned
    explicit operator unsigned() const;

    /// const access to factors map
    const map_type& factors() const
    { return Factors_; }
    
 private:

    map_type Factors_;
};

/// check if equal to unity
inline const bool is_one(const PrimeFactors& pf)
{ return pf.empty(); }

/// convert to string
std::string to_string(const PrimeFactors& pf);

/// \return prime number decomposition of product of arguments
PrimeFactors decompose(const std::vector<unsigned>& N);

/// \return prime number decomposition of product of arguments
template <typename ... Types>
PrimeFactors decompose(unsigned n, Types ... other_n)
{ std::vector<unsigned> V({n, other_n...}); return decompose(V); }

/// \return prime number decomposition of product of factorials of arguments
PrimeFactors decompose_factorial(const std::vector<unsigned>& N);

/// \return prime number decomposition of product of factorials of arguments
template <typename ... Types>
PrimeFactors decompose_factorial(unsigned n, Types ... other_n)
{ std::vector<unsigned> V({n, other_n...}); return decompose_factorial(V); }

/// equality operator
inline const bool operator==(const PrimeFactors& lhs, const PrimeFactors& rhs)
{ return lhs.factors() == rhs.factors(); }

/// remove factors with exponent = 0
PrimeFactors& remove_zeroes(PrimeFactors& PF);

/// remove common factors from arguments;
/// remove_zeroes is not called.
void remove_common(PrimeFactors& A, PrimeFactors& B);

/// remove common factors from arguments, storing into return value;
/// remove_zeroes is not called.
PrimeFactors common(PrimeFactors& A, PrimeFactors& B);

/// multiplication assignment
inline PrimeFactors& operator*=(PrimeFactors& lhs, const PrimeFactors& rhs)
{ for (const auto& p_e : rhs) lhs[p_e.first] += p_e.second; return lhs; }

/// multiplication
inline PrimeFactors operator*(PrimeFactors A, const PrimeFactors& B)
{ return A *= B; }

/// \return exponentiation of factorized unsigned
inline PrimeFactors pow(PrimeFactors pf, unsigned n)
{ for (auto& p_e : pf) p_e.second *= n; return pf; }

/// \return square root of factorized unsigned, throws if sqrt is not unsigned.
PrimeFactors sqrt(PrimeFactors pf);

/// factorize PrimeFactors into sqrt-able and non-sqrt-able portions;
/// note, no sqrt is applied
std::array<PrimeFactors, 2> factorize_sqrt(PrimeFactors pf);

/// \return sign factor
double sign_factor(int n, int d);

/// \class RationalNumber
/// \brief Ractional number represented by prime decompositions of numerator and denominator
/// \author <NAME>
class RationalNumber
{
public:

    /// constructor
    /// \param N exponents of the numerator's prime numbers
    /// \param D exponents of the denominator's prime numbers
    /// \param s sign of number (converted to +-1, 0, inf, or nan)
    RationalNumber(const PrimeFactors& N, const PrimeFactors& D, double s = 1)
        : Numerator_(N), Denominator_(D), Sign_(std::isfinite(s) ? sign_of(s) : s)
    {
        remove_common(Numerator_, Denominator_);
        remove_zeroes(Numerator_);
        remove_zeroes(Denominator_);
    }

    /// int constructor
    /// \param N numerator
    /// \param D denominator
    explicit RationalNumber(int N, int D = 1)
        : RationalNumber(decompose(abs(N)), decompose(abs(D)), sign_factor(N, D)) {}

    /// \return Numerator_
    const PrimeFactors& numerator() const
    { return Numerator_; }

    /// \return Denominator_
    const PrimeFactors& denominator() const
    { return Denominator_; }

    /// \return Sign_
    const double sign() const
    { return Sign_; }

    /// convert to double
    explicit operator double() const
    { return Sign_ * static_cast<unsigned>(Numerator_) / static_cast<unsigned>(Denominator_); }

    /// flip the sign
    void negate()
    { Sign_ *= -1; }
    
private:
    
    /// prime factorization of numerator
    PrimeFactors Numerator_;

    /// prime factorization of denominator
    PrimeFactors Denominator_;
    
    // Prefactor, including sign
    double Sign_{1};

};

/// convert to string as [sign][N]/[D]
std::string to_string(const RationalNumber& f);

/// convert to string as {N, D}
std::string to_header_string(const RationalNumber& f);

/// convert to detailed string
std::string to_detailed_string(const RationalNumber& f);

/// equality operator
inline const bool operator==(const RationalNumber& lhs, const RationalNumber& rhs)
{ return lhs.sign() == rhs.sign() and lhs.numerator() == rhs.numerator() and lhs.denominator() == rhs.denominator(); }

/// inequality operator
inline const bool operator!=(const RationalNumber& lhs, const RationalNumber& rhs)
{ return !(lhs == rhs); }

/// less than operator
const bool operator<(const RationalNumber& lhs, const RationalNumber& rhs);

/// \return whether RationalNumber is zero
inline const bool is_zero(const RationalNumber& f)
{ return f.sign() == 0; }

/// \return whether RationalNumber is unity
inline const bool is_one(const RationalNumber& f)
{ return is_one(f.numerator()) and is_one(f.denominator()) and f.sign() == 1; }

/// remove common factors from arguments;
void remove_common(RationalNumber& A, RationalNumber& B);

/// remove common factors from arguments, storing into return value;
RationalNumber common(RationalNumber& A, RationalNumber& B);

/// unary minus
inline RationalNumber operator-(const RationalNumber& f)
{ return RationalNumber(f.numerator(), f.denominator(), -f.sign()); }

/// addition assignment
RationalNumber& operator+=(RationalNumber& lhs, const RationalNumber& rhs);

/// addition
inline RationalNumber operator+(RationalNumber lhs, const RationalNumber& rhs)
{ return lhs += rhs; }

/// multiplication assignment
RationalNumber& operator*=(RationalNumber& lhs, const RationalNumber& rhs);

/// multiplication
inline RationalNumber operator*(RationalNumber lhs, const RationalNumber& rhs)
{ return lhs *= rhs; }

/// multiplication assignment
inline RationalNumber& operator*=(RationalNumber& lhs, const PrimeFactors& rhs)
{ return lhs = RationalNumber(lhs.numerator() * rhs, lhs.denominator(), lhs.sign()); }

/// multiplication
inline RationalNumber operator*(RationalNumber lhs, const PrimeFactors& rhs)
{ return lhs *= rhs; }

/// multiplication
inline RationalNumber operator*(const PrimeFactors& lhs, RationalNumber rhs)
{ return rhs *= lhs; }

/// division of two PrimeFactors
inline RationalNumber operator/(const PrimeFactors& N, const PrimeFactors& D)
{ return RationalNumber(N, D); }

/// multiple sign factors
double multiply_sign_factors(double s1, double s2);

/// invert sign
double invert_sign_factor(double s);

/// invert RationalNumber
inline RationalNumber invert(const RationalNumber& f)
{ return RationalNumber(f.denominator(), f.numerator(), invert_sign_factor(f.sign())); }

/// \return square root of absolute value of number with sign
/// preserved, potentially throws if sqrt is not int
inline RationalNumber sqrt(const RationalNumber& f)
{ return RationalNumber(sqrt(f.numerator()), sqrt(f.denominator()), f.sign()); }

/// factorize rational number into cleanly sqrt-able and
/// non-cleanly-sqrt-able portions; note the sqrt is not taken of
/// either and the sign is transfered to the sqrt-able portion.
/// \return array := [sqrt-able, non-sqrt-able]
std::array<RationalNumber, 2> factorize_sqrt(const RationalNumber& f);

/// convert to sqrt string := [sign][N_0]/[D_0] sqrt([N_1][D_1]);
std::string to_sqrt_string(const std::array<RationalNumber, 2>& f);

/// convert to sqrt string := [sign][N]/[D] sqrt([n][d]);
/// where clean sqrt is pulled out
inline std::string to_sqrt_string(const RationalNumber& f)
{ return to_sqrt_string(factorize_sqrt(f)); }


/// \return exponentiated rational number with sign preserved
inline RationalNumber pow(const RationalNumber& f, unsigned n)
{ return RationalNumber(pow(f.numerator(), n), pow(f.denominator(), n), f.sign()); }

/// \return absolute value of number
inline RationalNumber abs(const RationalNumber& f)
{ return RationalNumber(f.numerator(), f.denominator(), std::abs(f.sign())); }

}
 
#endif
<file_sep>/programs/CalcAmpl.cxx
// #include "TFhh.h"
// #include "TJSS.h"

#include "Contraction.h"
#include "WaveFunction.h"

#include <iostream>
#include <stdexcept>
#include <vector>

using namespace relhel;

//-------------------------
// printing function
std::string to_string(unsigned count, const std::vector<Contraction::Term::index>& I)
// { if (count == 0) return ""; return ((count > 1) ? (std::to_string(count) + " x ") : "") + "(" + to_string(I) + "); "; }
{ return std::to_string(count) + " x (" + to_string(I) + "); "; }

int main(int narg, char** args)
{

    if (narg < 4) {
        std::cout << "Usage: " << args[0] << " (Mother 2J) (Daughter 2J) (Daughter 2J) [L] [2S] [2M]" << std::endl;
        return 1;
    }

    WaveFunction phi(std::atoi(args[1]));
    WaveFunction psi1(std::atoi(args[2]));
    WaveFunction psi2(std::atoi(args[3]));

    if (narg < 7 and spin(phi) != 0)
        throw std::invalid_argument("M is ambiguous, you must state it explicitly; main");

    unsigned L = 0;
    if (narg < 5) {
        // can fix S?
        if (abs(spin(psi1) - spin(psi2)) == (spin(psi1) + spin(psi2))
            and abs(spin(phi) - spin(psi1) - spin(psi2)) == (spin(phi) + spin(psi1) + spin(psi2)))
            L = (spin(phi) + spin(psi1) + spin(psi2)) / 2;
        else
            throw std::invalid_argument("L is ambiguous, you must state it explicitly; main");
    } else
        L = std::atoi(args[4]);

    unsigned two_S = 0;
    if (narg < 6) {
        if (abs(spin(phi) - 2 * L) == (spin(phi) + 2 * L))
            two_S = spin(phi) + 2 * L;
        else if (abs(spin(psi1) - spin(psi2)) == (spin(psi1) + spin(psi2)))
            two_S = spin(psi1) + spin(psi2);
        else
            throw std::invalid_argument("S is ambiguous, you must state it explicitly; main");
    } else
        two_S = std::atoi(args[5]);

    int two_M = (narg < 7) ? 0 : std::atoi(args[6]);


    std::cout << "(" + spin_to_string(spin(phi)) + ") -> (" + spin_to_string(spin(psi1)) + ") + (" + spin_to_string(spin(psi2)) + ")"
              << ", L = " << L
              << ", S = " << spin_to_string(two_S)
              << ", M = " << spin_to_string(two_M)
              << std::endl;

    if (!triangle(spin(psi1), spin(psi2), two_S))
        throw std::invalid_argument("j1:j2:S triangle broken; main");

    if (!triangle(spin(phi), 2 * L, two_S))
        throw std::invalid_argument("j0:L:S triangle broken; main");
    
    CoupledWaveFunctions psi(psi1, psi2, two_S, two_M);

    OrbitalAngularMomentumWaveFunction chi(L);

    Contraction C(psi, chi, phi);

    GammaPolynomial P;
    
    for (unsigned psi_internal = 0; psi_internal <= std::min(rank(psi1), rank(psi2)); ++psi_internal) {

        auto K1 = C;
        
        K1.contract(contractions::psi1_psi2, psi_internal);
        if (K1.terms().empty()) continue;
        
        for (unsigned psi_chi = 0; psi_chi <= std::min(rank(psi), rank(chi)); ++psi_chi) {
            for (unsigned psi1_chi = 0; psi1_chi <= (rank(psi1) - psi_internal) and psi1_chi <= psi_chi; ++psi1_chi) {
                auto K2 = K1;
                
                K2.contract(contractions::psi1_chi, psi1_chi);
                if (K2.terms().empty()) continue;
                
                // limit psi2_chi to remaining rank of psi2 after contraction with psi1
                unsigned psi2_chi = std::min(psi_chi - psi1_chi, rank(psi2) - psi_internal);
                K2.contract(contractions::psi2_chi, psi2_chi);
                if (K2.terms().empty()) continue;

                for (unsigned psi_phi = 0; psi_phi <= std::min(rank(psi), rank(phi)); ++psi_phi) {
                    for (unsigned psi1_phi = 0; psi1_phi <= (rank(psi1) - psi_internal - psi1_chi) and psi1_phi <= psi_phi; ++psi1_phi) {
                        auto K3 = K2;
                        
                        K3.contract(contractions::psi1_phi, psi1_phi);
                        if (K3.terms().empty()) continue;

                        // limit psi2_phi to remaining rank of psi2 after contraction with psi1 and chi
                        unsigned psi2_phi = std::min(psi_phi - psi1_phi, rank(psi2) - psi_internal - psi2_chi);
                        K3.contract(contractions::psi2_phi, psi2_phi);
                        if (K3.terms().empty()) continue;

                        for (unsigned chi_phi = 0; chi_phi <= std::min(rank(phi) - psi1_phi - psi2_phi, rank(chi) - psi1_chi - psi2_chi); ++chi_phi) {
                            auto K4 = K3;

                            K4.contract(contractions::chi_phi, chi_phi);
                            if (K4.terms().empty()) continue;

                            auto I3 = triple_contraction(K4);
                            if (I3.size() == 3) {
                                K4.contract(I3);
                                if (K4.terms().empty()) continue;
                            }

                            if (rank(K4) == 0) {

                                // std::cout << "contractions: "
                                //           << to_string(psi_internal, contractions::psi1_psi2)
                                //           << to_string(psi1_chi,     contractions::psi1_chi)
                                //           << to_string(psi2_chi,     contractions::psi2_chi)
                                //           << to_string(psi1_phi,     contractions::psi1_phi)
                                //           << to_string(psi2_phi,     contractions::psi2_phi)
                                //           << to_string(chi_phi,      contractions::chi_phi)
                                //           << to_string(I3.size() == 3, I3)
                                //           << "produce " << K4.terms().size() << " terms" << std::endl;

                                P += gamma_polynomial(K4.terms());
                            }
                        }
                    }
                }
            }
        }
    }

    std::cout << to_string(P) << std::endl;
    
    return 0;
}

<file_sep>/include/MathUtils.h
#ifndef relhel__MathUtils_h
#define relhel__MathUtils_h

#include <cmath>
#include <cstdlib>
#include <stdexcept>
#include <string>
#include <vector>

namespace relhel {

/// \return whether value is odd
constexpr const bool is_odd(int v) noexcept
{ return v & 0x1; }

/// \return whether value is even
constexpr const bool is_even(int v) noexcept
{ return !is_odd(v); }

/// (-1)^n
constexpr int pow_negative_one(int n)
{ return is_odd(n) ? -1 : +1; }

/// \return whether three spins fulfil triangle relationship
inline const bool triangle(unsigned two_j1, unsigned two_j2, unsigned two_j3)
{ return two_j3 >= std::abs((int)two_j1 - (int)two_j2) and two_j3 <= (two_j1 + two_j2); }

/// \return vector of (twice) spins from |two_s1-two_s2| to (two_s1 + two_s2)
const std::vector<unsigned> triangle(unsigned two_s1, unsigned two_s2);

/// \return vector of (twice) spin projections from -two_j to two_j
const std::vector<int> projections(unsigned two_j);

/// \return vector of spin projections of a collection of particles
const std::vector<std::vector<int> > projections(const std::vector<unsigned>& two_J);

/// \return sign of argument
template <typename T>
constexpr const T sign_of(T t) noexcept
{ return ((T(0) < t) - (t < T(0))); }

/// \return sign_of(val) * sqrt(abs(val))
inline const double signed_sqrt(double val) noexcept
{ return sign_of(val) * std::sqrt(std::abs(val)); }

/// \return exponent string := "" (n == 0), "s" (n == 1), "s^n" (otherwise)
inline std::string exponential_string(std::string s, int n)
{ return n == 0 ? "" : (s + (n == 1 ? "" : "^" + std::to_string(n))); }

/// \return spin as string
inline std::string spin_to_string(int two_s)
{ return is_even(two_s) ? std::to_string(two_s / 2) : (std::to_string(two_s) + "/2"); }

}

#endif
<file_sep>/include/WaveFunction.h
#ifndef relhel__WaveFunction_h
#define relhel__WaveFunction_h

#include "RationalNumber.h"

#include <algorithm>
#include <cstdlib>
#include <numeric>
#include <string>
#include <vector>

namespace relhel {

/// vector of (twice) the spin projection of unit-spin wave functions
using WaveProduct = std::vector<int>;

/// \return whether all projections are +1, 0, or -1
inline const bool valid(const WaveProduct& v)
{ return std::all_of(v.begin(), v.end(), [](int m){return std::abs(m) <= 1;}); }

/// \return projection of WaveProduct
inline const int projection(const WaveProduct& v)
{ return std::accumulate(v.begin(), v.end(), 0); }

/// \return number of projection-zero waves in WaveProduct
inline const unsigned zeroes(const WaveProduct& v)
{ return std::count(v.begin(), v.end(), 0); }

/// \return rank of WaveProduct
inline const unsigned rank(const WaveProduct& w)
{ return w.size(); }

/// \return total spin of WaveProduct
inline const unsigned spin(const WaveProduct& v)
{ return 2 * v.size(); }

/// \return coefficient of WaveProduct
const RationalNumber squared_coefficient(const WaveProduct& v);

/// convert to string
std::string to_string(const WaveProduct& wp);

/// terms in sum of WaveProduct's (without coefficients)
using WaveProductSum = std::vector<WaveProduct>;

/// \return rank of WaveProductSum
inline const unsigned rank(const WaveProductSum& wps)
{ return wps.empty() ? 0 : rank(wps[0]); }

/// \return spin of WaveProductSum
inline const unsigned spin(const WaveProductSum& wps)
{ return wps.empty() ? 0 : spin(wps[0]); }

/// convert to string
std::string to_string(const WaveProductSum& wps);

/// holds wave functions
class WaveFunction
{
public:
    /// map of spin projection to WaveProductSum
    using map_type = std::map<int, WaveProductSum>;
    
    /// Constructor
    /// \param two_j (twice) the spin of wave function
    explicit WaveFunction(unsigned two_j);

    /// \return Projections_
    const map_type& projections() const
    { return Projections_; }
    
protected:

    /// map of spin projection to WaveProductSum
    map_type Projections_;
};

/// \return spin of WaveFunction
const unsigned spin(const WaveFunction& wf);

/// \return rank of WaveFunction
const unsigned rank(const WaveFunction& wf);

/// convert to string
std::string to_string(const WaveFunction& wf);

/// holds orbital angular momentum wave function
struct OrbitalAngularMomentumWaveFunction : public WaveFunction {
    /// Constructor
    /// \param l orbital angular momentum
    OrbitalAngularMomentumWaveFunction(unsigned l) : WaveFunction(2 * l)
    {
        for (auto it = Projections_.begin(); it != Projections_.end(); )
            if (it->first != 0)
                it = Projections_.erase(it);
            else
                ++it;
    }
};

/// Coupling of two WaveFunction's
class CoupledWaveFunctions
{
public:
    
    /// Constructor
    CoupledWaveFunctions(const WaveFunction& phi1, const WaveFunction& phi2, unsigned two_S, int delta);

    /// \return array of coupled WaveFunction's
    const std::array<WaveFunction, 2>& phi() const
    { return Phi_; }

    /// \return total coupled spin
    const unsigned twoS() const
    { return TwoS_; }

    /// \return spin projection
    const int delta() const
    { return Delta_; }
    
private:

    /// First WaveFunction to couple
    std::array<WaveFunction, 2> Phi_;

    /// total coupled spin
    unsigned TwoS_;

    /// projection of total coupled spin
    int Delta_;
    
};

/// \return rank of coupled wave functions
inline unsigned rank(const CoupledWaveFunctions& cwf)
{ return std::accumulate(cwf.phi().begin(), cwf.phi().end(), 0u, [](unsigned r, const WaveFunction& w){return r + rank(w);}); }

}

#endif
<file_sep>/programs/CMakeLists.txt
include_directories(${RelHel_SOURCE_DIR}/include)

set(PROGRAMS
  test
  CalcAmpl
  )

foreach(program ${PROGRAMS})
	add_executable(${program} ${program}.cxx)
	target_link_libraries(${program} RelHel)
endforeach(program)
<file_sep>/src/ClebschGordan.cxx
#include "ClebschGordan.h"

#include "MathUtils.h"

#include <algorithm>
#include <cmath>

namespace relhel {

//-------------------------
std::string ClebschGordan::to_string(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J, int two_M)
{
    return std::string("(") + spin_to_string(two_j1) + " " + spin_to_string(two_m1)
           + ", " + spin_to_string(two_j2) + " " + spin_to_string(two_m2)
           + " | " + spin_to_string(two_J) + " " + spin_to_string(two_M) + ")";
}

//-------------------------
const bool ClebschGordan::consistent(unsigned two_J, int two_M)
{
    return (std::abs(two_M) <= (int)two_J) and is_even(two_J + two_M);
}
    
//-------------------------
const bool ClebschGordan::nonzero(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J, int two_M)
{
    // and that (j1+j2) and J are consistent
    if (is_odd(two_J + two_j1 + two_j2))
        return false;

    // check input spin-projection compatibilities
    if (!consistent(two_j1, two_m1) or !consistent(two_j2, two_m2) or !consistent(two_J,  two_M ))
        return false;

    // check input spin-projections
    if (two_M != two_m1 + two_m2)
        return false;

    // check whether J lies between |j1 - j2| and (j1 + j2)
    if ((int)two_J < std::abs((int)two_j1 - (int)two_j2) or two_J > (two_j1 + two_j2))
        return false;

    // when either daughter spin is zero
    if (two_j1 == 0 and two_j2 != two_J)
        return false;
    if (two_j2 == 0 and two_j1 != two_J)
        return false;

    // check for (j1 0 j2 0 | J 0), j1 + j2 + J is odd
    if (two_m1 == 0 and two_m2 == 0 and is_odd((two_j1 + two_j2 + two_J) / 2))
        return false;

    // (3/2 +-1/2, 3/2 +-1/2 | 2 +-1) == 0
    if (two_j1 == 3 and std::abs(two_m1) == 1 and two_j2 == 3 and two_m2 == two_m1 and two_J == 4)
        return false;

    // (2 +-1, 3/2 -+1/2 | 3/2 +-1/2) == 0
    if (two_j1 == 4 and std::abs(two_m1) == 2 and two_j2 == 3 and two_m2 == -two_m1 / 2 and two_J == 3)
        return false;

    // (2 +-1, 2 +-1 | 3 +-2) == 0
    if (two_j1 == 4 and std::abs(two_m1) == 2 and two_j2 == 4 and two_m2 == two_m1 and two_J == 6)
        return false;

    return true;
}
    
//-------------------------
const RationalNumber ClebschGordan::squared_coefficient(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J, int two_M)
{
    if (!nonzero(two_j1, two_m1, two_j2, two_m2, two_J, two_M))
        return RationalNumber(0);

    // simple case of spin zero (must be 1 since check above already would have found 0)
    if (two_j1 == 0 or two_j2 == 0)
        return RationalNumber(1);

    // z range dictated by factorials in denominator ( 1/n! = 0 when n < 0)
    unsigned z_min = std::max<int>({0, (int)two_j2 - two_m1 - (int)two_J, (int)two_j1 + two_m2 - (int)two_J}) / 2;
    unsigned z_max = std::min<int>({(int)two_j1 + (int)two_j2 - (int)two_J, (int)two_j1 - two_m1, (int)two_j2 + two_m2}) / 2;

    RationalNumber z_sum(0);
    for (unsigned z = z_min; z <= z_max; ++z)
        // z'th term := (-)^z / z! / (j1+j2-J-z)! / (j1-m1-z)! / (j2+m2-z)! / (J-j2+m1+z)! / (J-j1-m2+z)!
        z_sum += RationalNumber(decompose(1),
                                decompose_factorial(z) *
                                decompose_factorial(((int)two_j1 + (int)two_j2 - (int)two_J) / 2 - z) *
                                decompose_factorial(((int)two_j1 - two_m1) / 2 - z) *
                                decompose_factorial(((int)two_j2 + two_m2) / 2 - z) *
                                decompose_factorial(((int)two_J - (int)two_j2 + two_m1) / 2 + z) *
                                decompose_factorial(((int)two_J - (int)two_j1 - two_m2) / 2 + z),
                                static_cast<double>(pow_negative_one(z)));
    
    // C-G coef = sqrt( (2J+1) (j1+j2-J)! (j1-j2+J)! (j2-j1+J)! / (J+j1+j2+1)! )
    //          * sqrt( (j1+m1)! (j1-m1)! (j2+m2)! (j2-m2)! (J+M)! (J-M)! )
    //          * z_sum
    return pow(z_sum, 2)
        * RationalNumber(decompose(two_J + 1)
                         * decompose_factorial(((int)two_j1 + (int)two_j2 - (int)two_J) / 2)
                         * decompose_factorial(((int)two_j1 - (int)two_j2 + (int)two_J) / 2)
                         * decompose_factorial(((int)two_j2 - (int)two_j1 + (int)two_J) / 2)
                         * decompose_factorial(((int)two_j1 + two_m1) / 2)
                         * decompose_factorial(((int)two_j1 - two_m1) / 2)
                         * decompose_factorial(((int)two_j2 + two_m2) / 2)
                         * decompose_factorial(((int)two_j2 - two_m2) / 2)
                         * decompose_factorial(((int)two_J + two_M) / 2)
                         * decompose_factorial(((int)two_J - two_M) / 2),
                         decompose_factorial(((int)two_j1 + (int)two_j2 + (int)two_J) / 2 + 1));
}
}
<file_sep>/src/WaveFunction.cxx
#include "WaveFunction.h"

#include <stdexcept>

namespace relhel {

//-------------------------
const RationalNumber squared_coefficient(const WaveProduct& v)
{
    unsigned j = spin(v) / 2;
    int two_m = projection(v);
    if (is_odd(two_m))
        throw std::runtime_error("projection is not integer; squared_coefficient(...)");
    int m = two_m / 2;
    return decompose_factorial(j + m) * decompose_factorial(j - m) * pow(decompose(2), zeroes(v)) / decompose_factorial(2 * j);
}

//-------------------------
std::string to_string(const WaveProduct& wp)
{
    return std::accumulate(wp.begin(), wp.end(), std::string("("), [](std::string& s, int m){return s += (m > 0) ? "+" : ((m < 0) ? "-" : "0");}) + ")";
}

//-------------------------
std::string to_string(const WaveProductSum& wps)
{
    return std::accumulate(wps.begin(), wps.end(), std::string(), [](std::string& s, const WaveProduct& w){return s += ", " + to_string(w);}).erase(0, 2);
}

//-------------------------
WaveFunction::WaveFunction(unsigned two_j)
{
    if (is_odd(two_j))
        throw std::invalid_argument("currently only supporting integer spins; WaveFunction::WaveFunction");

    WaveProduct WP(two_j / 2, -2);

    if (two_j == 0)
        Projections_[0] = WaveProductSum(1, WP);
    else {
         // loop through permutations using odometer method
        while (WP.back() <= 2) {
            // add product into sum for corresponding projection
            Projections_[projection(WP)].push_back(WP);
            // tick over the first projection
            WP[0] += 2;
            // check what projections must be reset or ticked over
            for (size_t i = 0; (i < WP.size() - 1) and (WP[i] > 2); ++i) {
                WP[i] = - 2;
                WP[i + 1] += 2;
            }
        }
    }
}

//-------------------------
const unsigned spin(const WaveFunction& wf)
{
    // look for non-empty spin projection in WaveFunction
    auto it = std::find_if(wf.projections().begin(), wf.projections().end(), [](const WaveFunction::map_type::value_type& m_wps){return !m_wps.second.empty();});
    // if none, return 0
    if (it == wf.projections().end())
        return 0;
    // return spin of WaveProductSum
    return spin(it->second);
}

//-------------------------
const unsigned rank(const WaveFunction& wf)
{
    // look for non-empty spin projection in WaveFunction
    auto it = std::find_if(wf.projections().begin(), wf.projections().end(), [](const WaveFunction::map_type::value_type& m_wps){return !m_wps.second.empty();});
    // if none, return 0
    if (it == wf.projections().end())
        return 0;
    // return tank of WaveProductSum
    return rank(it->second);
}

//-------------------------
// helper function
std::string spin_projection_string(unsigned two_j, int two_m)
{ return std::string("|") + spin_to_string(two_j) + ", " + spin_to_string(two_m) + ">"; }

//-------------------------
std::string to_string(const WaveFunction& wf)
{
    auto two_j = spin(wf);
    return std::accumulate(wf.projections().begin(), wf.projections().end(), std::string(),
                           [&](std::string& s, const WaveFunction::map_type::value_type& m_wps)
                           {return s += "\n" + spin_projection_string(two_j, m_wps.first) + " = " + to_string(m_wps.second);}).erase(0, 1);
}

//-------------------------
CoupledWaveFunctions::CoupledWaveFunctions(const WaveFunction& phi1, const WaveFunction& phi2, unsigned two_s, int delta)
    : Phi_({phi1, phi2}), TwoS_(two_s), Delta_(delta)
{
    auto two_j0 = spin(Phi_[0]);
    auto two_j1 = spin(Phi_[1]);
    
    if (!triangle(two_j0, two_j1, TwoS_))
        throw std::invalid_argument("invalid spins for coupling: triangle("
                                    + spin_to_string(two_j0) + ", "
                                    + spin_to_string(two_j1) + ", "
                                    + spin_to_string(TwoS_) + ") = false; get_spin_coupled_tensor_sum(...)");
    
    if (is_odd(two_j0 + two_j1 + Delta_) or std::abs(Delta_) > (two_j0 + two_j1))
        throw std::invalid_argument("invalid spin projection for coupled state; get_spin_coupled_tensor_sum(...)");
}

}
<file_sep>/include/ClebschGordan.h
/// \file
/// \brief code adapted from YAP (www.github.com/YAP/YAP)
/// \author <NAME>

#ifndef __ClebschGordan__h
#define __ClebschGordan__h

#include "RationalNumber.h"

#include <string>

namespace relhel {

/// \namespace ClebschGordan
/// \brief Clebsch-Gordan coefficient calculation
/// \author <NAME>
namespace ClebschGordan {

/// \return Clebsch-Gordan coefficient string
/// \param two_j1 2*spin of first particle
/// \param two_m1 2*spin-projection of first particle
/// \param two_j2 2*spin of second particle
/// \param two_m2 2*spin-projection of second particle
/// \param two_J  2*spin of composite system
/// \param two_M  2*spin-projection of composite system
std::string to_string(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J, int two_M);

/// \return Clebsch-Gordan coefficient string
/// \param two_j1 2*spin of first particle
/// \param two_m1 2*spin-projection of first particle
/// \param two_j2 2*spin of second particle
/// \param two_m2 2*spin-projection of second particle
/// \param two_J  2*spin of composite system
inline std::string to_string(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J)
{ return to_string(two_j1, two_m1, two_j2, two_m2, two_J, two_m1 + two_m2); }

/// \return consistency of spin and spin projection
/// \param two_J 2*spin
/// \param two_M 2*spin-projection
const bool consistent(unsigned two_J, int two_M);

/// \return Whether Clebsch-Gordan coefficient is nonzero
/// \param two_j1 2*spin of first particle
/// \param two_m1 2*spin-projection of first particle
/// \param two_j2 2*spin of second particle
/// \param two_m2 2*spin-projection of second particle
/// \param two_J  2*spin of composite system
/// \param two_M  2*spin-projection of composite system
const bool nonzero(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J, int two_M);

/// \return Whether Clebsch-Gordan coefficient is nonzero, with M := m1 + m2.
/// \param two_j1 2*spin of first particle
/// \param two_m1 2*spin-projection of first particle
/// \param two_j2 2*spin of second particle
/// \param two_m2 2*spin-projection of second particle
/// \param two_J  2*spin of composite system
inline const bool nonzero(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J)
{ return nonzero(two_j1, two_m1, two_j2, two_m2, two_J, two_m1 + two_m2); }

/// Calculates _signed_ _squared_ C-G coefficienty. For return value
/// C, the Clebsch-Gordan coefficient is sign(C)sqrt(abs(C))
/// \return signed squared Clebsch-Gordan coefficient (j1 m1 j2 m2 | J M)
/// Implemented from Eq. (16) from <NAME>, "Theory of Complex Spectra. II", Phys. Rev. 62, 438 (1942)
/// \param two_j1 2*spin of first particle
/// \param two_m1 2*spin-projection of first particle
/// \param two_j2 2*spin of second particle
/// \param two_m2 2*spin-projection of second particle
/// \param two_J  2*spin of composite system
/// \param two_M  2*spin-projection of composite system
const RationalNumber squared_coefficient(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J, int two_M);

/// Calculates _signed_ _squared_ C-G coefficienty. For return value
/// C, the Clebsch-Gordan coefficient is sign(C)sqrt(abs(C))
/// \return signed squared Clebsch-Gordan coefficient (j1 m1 j2 m2 | J M)
/// Implemented from Eq. (16) from <NAME>, "Theory of Complex Spectra. II", Phys. Rev. 62, 438 (1942)
/// \param two_j1 2*spin of first particle
/// \param two_m1 2*spin-projection of first particle
/// \param two_j2 2*spin of second particle
/// \param two_m2 2*spin-projection of second particle
/// \param two_J  2*spin of composite system
inline const RationalNumber squared_coefficient(unsigned two_j1, int two_m1, unsigned two_j2, int two_m2, unsigned two_J)
{ return squared_coefficient(two_j1, two_m1, two_j2, two_m2, two_J, two_m1 + two_m2); }

}

}

#endif
<file_sep>/include/Contraction.h
#ifndef relhel__Contraction_h
#define relhel__Contraction_h

#include "RationalNumber.h"
#include "WaveFunction.h"

#include <array>
#include <string>
#include <vector>

namespace relhel {

/// performs contractions over tensors
class Contraction
{
public:

    /// Contraction Term
    struct Term {

        /// index for WaveProducts
        enum index : size_t {
            /// first daughter in final state
            psi1 = 0,
            /// second daughter in final state
            psi2 = 1,
            /// orbital angular momentum
            chi = 2,
            /// (complex conjugated) initial state
            phi = 3
        };

        /// Daughter WaveProducts
        std::vector<WaveProduct> WaveProducts;
        
        /// Square of coefficients
        RationalNumber SquaredCoefficient;

        /// exponents of gamma factors of contributing components
        std::vector<int> GammaExponents;
        
        /// Constructor
        Term(const WaveProduct& wp1, const WaveProduct& wp2, const WaveProduct& wp_chi, const WaveProduct& phi, const RationalNumber& c2)
        : WaveProducts({wp1, wp2, wp_chi, phi}), SquaredCoefficient(c2), GammaExponents(WaveProducts.size(), 0) {}

        /// contract components
        /// \param I vector of components to contract.
        /// if I is empty and only 3 components remain to be contracted, they are chosen
        void contract(std::vector<index> I = {});
    };

    /// Constructor
    /// \param psi coupled daughter state WaveFunction's
    /// \param chi orbital angular momentum WaveFunction
    Contraction(const CoupledWaveFunctions& psi, const OrbitalAngularMomentumWaveFunction& chi, const WaveFunction& phi);

    /// \return Terms_
    const std::vector<Term>& terms() const
    { return Terms_; }

    // prune off zero terms
    void prune()
    { Terms_.erase(std::remove_if(Terms_.begin(), Terms_.end(), [](const Term& t){return is_zero(t.SquaredCoefficient);}), Terms_.end()); }

    /// contract daughters
    void contract(const std::vector<Term::index>& I = {})
    { std::for_each(Terms_.begin(), Terms_.end(), [&I](Term& t){t.contract(I);}); prune(); }

    /// contract daughters
    void contract(const std::vector<Term::index>& I, unsigned n)
    { for (unsigned i = 0; i < n and !Terms_.empty(); ++i) contract(I); }

private:
    /// Terms in contraction
    std::vector<Term> Terms_;
};

namespace contractions {
    extern const std::vector<Contraction::Term::index> psi1_psi2;
    extern const std::vector<Contraction::Term::index> psi1_chi;
    extern const std::vector<Contraction::Term::index> psi2_chi;
    extern const std::vector<Contraction::Term::index> psi1_phi;
    extern const std::vector<Contraction::Term::index> psi2_phi;
    extern const std::vector<Contraction::Term::index> chi_phi;
}

/// \return vector of only the requested gamma exponents
const std::vector<int> gamma_exponents(const Contraction::Term& t, const std::vector<Contraction::Term::index>& I);

/// \typedef GammaPolynomial
 using GammaPolynomial = std::map<std::vector<int>, std::vector<std::pair<int, RationalNumber> > >;

/// add polynomials
GammaPolynomial& operator+=(GammaPolynomial& lhs, const GammaPolynomial& rhs);
 
/// add polynomials
const GammaPolynomial operator+(GammaPolynomial lhs, const GammaPolynomial& rhs)
{ return lhs += rhs; }
 
/// \return string of polynomial term
std::string to_string(const GammaPolynomial::value_type& ges_coefs);

/// \return string of gamma polynomial
std::string to_string(const GammaPolynomial& gp);
 
/// \return GammaPolynomial
const GammaPolynomial gamma_polynomial(const std::vector<Contraction::Term>& terms);
 
/// \return rank of contraction term
inline unsigned rank(const Contraction::Term& t)
{ return std::accumulate(t.WaveProducts.begin(), t.WaveProducts.end(), 0u, [](unsigned r, const WaveProduct& w){return r + rank(w);}); }

/// \return rank of contraction
inline unsigned rank(const Contraction& C)
{ return C.terms().empty() ? 0 : rank(C.terms()[0]); }

/// \return indices for triple contraction; empty if none possible
const std::vector<Contraction::Term::index> triple_contraction(const Contraction::Term& t);

/// \return indices for triple contraction; empty if none possible
inline const std::vector<Contraction::Term::index> triple_contraction(const Contraction& C)
{ if(C.terms().empty()) return std::vector<Contraction::Term::index>(); return triple_contraction(C.terms()[0]); }

/// convert to string
std::string to_string(const Contraction::Term::index i);

/// convert to string
std::string to_string(const std::vector<Contraction::Term::index>& I);
 
/// convert to string
std::string to_string(const Contraction::Term& term);

/// convert to string
std::string to_string(const std::vector<Contraction::Term>& terms);

}

#endif
<file_sep>/src/CMakeLists.txt
set(INCLUDE_DIRECTORIES
	${RelHel_SOURCE_DIR}/include)
include_directories(${INCLUDE_DIRECTORIES})

set(RelHel_SOURCES
  ClebschGordan.cxx
  Contraction.cxx
  MathUtils.cxx
  RationalNumber.cxx
  WaveFunction.cxx
)

add_library(RelHel SHARED ${RelHel_SOURCES})

# install destinations can be passed via the command line:
# cmake -DLIBRARY_OUTPUT_DIRECTORY:PATH=<lib_path> <path_to_CMakeLists.tex>
# otherwise, default LD_LIBRARY_PATH
if(NOT DEFINED LIBRARY_OUTPUT_DIRECTORY)
	set(LIBRARY_OUTPUT_DIRECTORY ${RelHel_SOURCE_DIR}/lib/${CMAKE_BUILD_TYPE})
endif()

if(NOT DEFINED INCLUDE_OUTPUT_DIRECTORY)
	set(INCLUDE_OUTPUT_DIRECTORY ${RelHel_SOURCE_DIR}/include/RelHel)
endif()

install(TARGETS RelHel LIBRARY DESTINATION ${LIBRARY_OUTPUT_DIRECTORY})

# Matches all the headers in ${RelHelDID}/include and its subdirs
file(GLOB_RECURSE
	 INSTALL_INCLUDES ${RelHel_SOURCE_DIR}/include/*.h)

# message(STATUS "${INSTALL_INCLUDES}")

# install(FILE       ${INSTALL_INCLUDES}
#	    DESTINATION ${INCLUDE_OUTPUT_DIRECTORY}
#	)
