initial commit.

diff --git a/.gitignore b/.gitignore
new file mode 100644 (file)
index 0000000..7893817
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+CMakeCache.txt
+CMakeFiles
+*~
+*.a
+cmake_install.cmake
+Makefile
+clebschgordan

diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644 (file)
index 0000000..70efbe2
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,20 @@
+project(SU2CLEBSCH)
+
+cmake_minimum_required(VERSION 3.0)
+
+set(CMAKE_BUILD_TYPE Release)
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+
+#set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR})
+
+find_package(LAPACK REQUIRED)
+
+add_library(clebsch clebsch.cpp)
+target_link_libraries(clebsch lapack)
+
+add_executable(clebschgordan clebschgordan.cpp)
+target_link_libraries(clebschgordan clebsch)
+
+add_library(su2clebsch su2clebsch.cpp)
+target_link_libraries(su2clebsch clebsch)

diff --git a/clebsch.cpp b/clebsch.cpp
new file mode 100644 (file)
index 0000000..f5f852d
--- /dev/null
+++ b/clebsch.cpp
@@ -0,0 +1,816 @@
+#include "clebsch.h"
+
+// implementation of "binomial_t" starts here
+
+clebsch::binomial_t binomial;
+
+int clebsch::binomial_t::operator()(int n, int k) {
+    if (N <= n) {
+        for (cache.resize((n + 1) * (n + 2) / 2); N <= n; ++N) {
+            cache[N * (N + 1) / 2] = cache[N * (N + 1) / 2 + N] = 1;
+            for (int k = 1; k < N; ++k) {
+                cache[N * (N + 1) / 2 + k] = cache[(N - 1) * N / 2 + k]
+                                           + cache[(N - 1) * N / 2 + k - 1];
+            }
+        }
+    }
+
+    return cache[n * (n + 1) / 2 + k];
+}
+
+// implementation of "weight" starts here
+
+clebsch::weight::weight(int N) : elem(N), N(N) {}
+
+clebsch::weight::weight(int N, int index) : elem(N, 0), N(N) {
+    for (int i = 0; index > 0 && i < N; ++i) {
+        for (int j = 1; binomial(N - i - 1 + j, N - i - 1) <= index; j <<= 1) {
+            elem[i] = j;
+        }
+
+        for (int j = elem[i] >> 1; j > 0; j >>= 1) {
+            if (binomial(N - i - 1 + (elem[i] | j), N - i - 1) <= index) {
+                elem[i] |= j;
+            }
+        }
+
+        index -= binomial(N - i - 1 + elem[i]++, N - i - 1);
+    }
+}
+
+clebsch::weight &clebsch::weight::operator=(const clebsch::weight &w) {
+    int &n = const_cast<int &>(N);
+    elem = w.elem;
+    n = w.N;
+    return *this;
+}
+
+int &clebsch::weight::operator()(int k) {
+    assert(1 <= k && k <= N);
+    return elem[k - 1];
+}
+
+const int &clebsch::weight::operator()(int k) const {
+    assert(1 <= k && k <= N);
+    return elem[k - 1];
+}
+
+bool clebsch::weight::operator<(const weight &w) const {
+    assert(w.N == N);
+    for (int i = 0; i < N; ++i) {
+        if (elem[i] - elem[N - 1] != w.elem[i] - w.elem[N - 1]) {
+            return elem[i] - elem[N - 1] < w.elem[i] - w.elem[N - 1];
+        }
+    }
+    return false;
+}
+
+bool clebsch::weight::operator==(const weight &w) const {
+    assert(w.N == N);
+
+    for (int i = 1; i < N; ++i) {
+        if (w.elem[i] - w.elem[i - 1] != elem[i] - elem[i - 1]) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+clebsch::weight clebsch::weight::operator+(const weight &w) const {
+    weight result(N);
+
+    transform(elem.begin(), elem.end(), w.elem.begin(), result.elem.begin(), std::plus<int>());
+
+    return result;
+}
+
+int clebsch::weight::index() const {
+    int result = 0;
+
+    for (int i = 0; elem[i] > elem[N - 1]; ++i) {
+        result += binomial(N - i - 1 + elem[i] - elem[N - 1] - 1, N - i - 1);
+    }
+
+    return result;
+}
+
+long long clebsch::weight::dimension() const {
+    long long numerator = 1, denominator = 1;
+
+    for (int i = 1; i < N; ++i) {
+        for (int j = 0; i + j < N; ++j) {
+            numerator *= elem[j] - elem[i + j] + i;
+            denominator *= i;
+        }
+    }
+
+    return numerator / denominator;
+}
+
+// implementation of "pattern" starts here
+
+clebsch::pattern::pattern(const pattern &p) : elem(p.elem), N(p.N) {}
+
+clebsch::pattern::pattern(const weight &irrep, int index) :
+        elem((irrep.N * (irrep.N + 1)) / 2), N(irrep.N) {
+    for (int i = 1; i <= N; ++i) {
+        (*this)(i, N) = irrep(i);
+    }
+
+    for (int l = N - 1; l >= 1; --l) {
+        for (int k = 1; k <= l; ++k) {
+            (*this)(k, l) = (*this)(k + 1, l + 1);
+        }
+    }
+
+    while (index-- > 0) {
+        bool b = ++(*this);
+
+        assert(b);
+    }
+}
+
+int &clebsch::pattern::operator()(int k, int l) {
+    return elem[(N * (N + 1) - l * (l + 1)) / 2 + k - 1];
+}
+
+const int &clebsch::pattern::operator()(int k, int l) const {
+    return elem[(N * (N + 1) - l * (l + 1)) / 2 + k - 1];
+}
+
+bool clebsch::pattern::operator++() {
+    int k = 1, l = 1;
+
+    while (l < N && (*this)(k, l) == (*this)(k, l + 1)) {
+        if (--k == 0) {
+            k = ++l;
+        }
+    }
+
+    if (l == N) {
+        return false;
+    }
+
+    ++(*this)(k, l);
+
+    while (k != 1 || l != 1) {
+        if (++k > l) {
+            k = 1;
+            --l;
+        }
+
+        (*this)(k, l) = (*this)(k + 1, l + 1);
+    }
+
+    return true;
+}
+
+bool clebsch::pattern::operator--() {
+    int k = 1, l = 1;
+
+    while (l < N && (*this)(k, l) == (*this)(k + 1, l + 1)) {
+        if (--k == 0) {
+            k = ++l;
+        }
+    }
+
+    if (l == N) {
+        return false;
+    }
+
+    --(*this)(k, l);
+
+    while (k != 1 || l != 1) {
+        if (++k > l) {
+            k = 1;
+            --l;
+        }
+
+        (*this)(k, l) = (*this)(k, l + 1);
+    }
+
+    return true;
+}
+
+int clebsch::pattern::index() const {
+    int result = 0;
+
+    for (pattern p(*this); --p; ++result) {}
+
+    return result;
+}
+
+clebsch::weight clebsch::pattern::get_weight() const {
+    clebsch::weight result(N);
+
+    for (int prev = 0, l = 1; l <= N; ++l) {
+        int now = 0;
+
+        for (int k = 1; k <= l; ++k) {
+            now += (*this)(k, l);
+        }
+
+        result(l) = now - prev;
+        prev = now;
+    }
+
+    return result;
+}
+
+double clebsch::pattern::lowering_coeff(int k, int l) const {
+    double result = 1.0;
+
+    for (int i = 1; i <= l + 1; ++i) {
+        result *= (*this)(i, l + 1) - (*this)(k, l) + k - i + 1;
+    }
+    
+    for (int i = 1; i <= l - 1; ++i) {
+        result *= (*this)(i, l - 1) - (*this)(k, l) + k - i;
+    }
+
+    for (int i = 1; i <= l; ++i) {
+        if (i == k) continue;
+        result /= (*this)(i, l) - (*this)(k, l) + k - i + 1;
+        result /= (*this)(i, l) - (*this)(k, l) + k - i;
+    }
+
+    return std::sqrt(-result);
+}
+
+double clebsch::pattern::raising_coeff(int k, int l) const {
+    double result = 1.0;
+
+    for (int i = 1; i <= l + 1; ++i) {
+        result *= (*this)(i, l + 1) - (*this)(k, l) + k - i;
+    }
+
+    for (int i = 1; i <= l - 1; ++i) {
+        result *= (*this)(i, l - 1) - (*this)(k, l) + k - i - 1;
+    }
+
+    for (int i = 1; i <= l; ++i) {
+        if (i == k) continue;
+        result /= (*this)(i, l) - (*this)(k, l) + k - i;
+        result /= (*this)(i, l) - (*this)(k, l) + k - i  - 1;
+    }
+
+    return std::sqrt(-result);
+}
+
+// implementation of "decomposition" starts here
+
+clebsch::decomposition::decomposition(const weight &factor1, const weight &factor2) :
+        N(factor1.N), factor1(factor1), factor2(factor2) {
+    assert(factor1.N == factor2.N);
+    std::vector<clebsch::weight> result;
+    pattern low(factor1), high(factor1);
+    weight trial(factor2);
+    int k = 1, l = N;
+
+    do {
+        while (k <= N) {
+            --l;
+            if (k <= l) {
+                low(k, l) = std::max(high(k + N - l, N), high(k, l + 1) + trial(l + 1) - trial(l));
+                high(k, l) = high(k, l + 1);
+                if (k > 1 && high(k, l) > high(k - 1, l - 1)) {
+                    high(k, l) = high(k - 1, l - 1);
+                }
+                if (l > 1 && k == l && high(k, l) > trial(l - 1) - trial(l)) {
+                    high(k, l) = trial(l - 1) - trial(l);
+                }
+                if (low(k, l) > high(k, l)) {
+                    break;
+                }
+                trial(l + 1) += high(k, l + 1) - high(k, l);
+            } else {
+                trial(l + 1) += high(k, l + 1);
+                ++k;
+                l = N;
+            }
+        }
+
+        if (k > N) {
+            result.push_back(trial);
+            for (int i = 1; i <= N; ++i) {
+                result.back()(i) -= result.back()(N);
+            }
+        } else {
+            ++l;
+        }
+
+        while (k != 1 || l != N) {
+            if (l == N) {
+                l = --k - 1;
+                trial(l + 1) -= high(k, l + 1);
+            } else if (low(k, l) < high(k, l)) {
+                --high(k, l);
+                ++trial(l + 1);
+                break;
+            } else {
+                trial(l + 1) -= high(k, l + 1) - high(k, l);
+            }
+            ++l;
+        }
+    } while (k != 1 || l != N);
+
+    sort(result.begin(), result.end());
+    for (std::vector<clebsch::weight>::iterator it = result.begin(); it != result.end(); ++it) {
+        if (it != result.begin() && *it == weights.back()) {
+            ++multiplicities.back();
+        } else {
+            weights.push_back(*it);
+            multiplicities.push_back(1);
+        }
+    }
+}
+
+int clebsch::decomposition::size() const {
+    return weights.size();
+}
+
+const clebsch::weight &clebsch::decomposition::operator()(int j) const {
+    return weights[j];
+}
+
+int clebsch::decomposition::multiplicity(const weight &irrep) const {
+    assert(irrep.N == N);
+    std::vector<clebsch::weight>::const_iterator it
+        = std::lower_bound(weights.begin(), weights.end(), irrep);
+
+    return it != weights.end() && *it == irrep ? multiplicities[it - weights.begin()] : 0;
+}
+
+// implementation of "index_adapter" starts here
+
+clebsch::index_adapter::index_adapter(const clebsch::decomposition &decomp) :
+        N(decomp.N),
+        factor1(decomp.factor1.index()),
+        factor2(decomp.factor2.index()) {
+    for (int i = 0, s = decomp.size(); i < s; ++i) {
+        indices.push_back(decomp(i).index());
+        multiplicities.push_back(decomp.multiplicity(decomp(i)));
+    }
+}
+
+int clebsch::index_adapter::size() const {
+    return indices.size();
+}
+
+int clebsch::index_adapter::operator()(int j) const {
+    return indices[j];
+}
+
+int clebsch::index_adapter::multiplicity(int irrep) const {
+    std::vector<int>::const_iterator it = std::lower_bound(indices.begin(), indices.end(), irrep);
+
+    return it != indices.end() && *it == irrep ? multiplicities[it - indices.begin()] : 0;
+}
+
+// implementation of "clebsch" starts here
+
+void clebsch::coefficients::set(int factor1_state,
+                                int factor2_state,
+                                int multiplicity_index,
+                                int irrep_state,
+                                double value) {
+    assert(0 <= factor1_state && factor1_state < factor1_dimension);
+    assert(0 <= factor2_state && factor2_state < factor2_dimension);
+    assert(0 <= multiplicity_index && multiplicity_index < multiplicity);
+    assert(0 <= irrep_state && irrep_state < irrep_dimension);
+
+    int coefficient_label[] = { factor1_state,
+                                factor2_state,
+                                multiplicity_index,
+                                irrep_state };
+    clzx[std::vector<int>(coefficient_label, coefficient_label
+            + sizeof coefficient_label / sizeof coefficient_label[0])] = value;
+}
+
+void clebsch::coefficients::highest_weight_normal_form() {
+    int hws = irrep_dimension - 1;
+
+    // bring CGCs into reduced row echelon form
+    for (int h = 0, i = 0; h < multiplicity - 1 && i < factor1_dimension; ++i) {
+        for (int j = 0; h < multiplicity - 1 && j < factor2_dimension; ++j) {
+            int k0 = h;
+
+            for (int k = h + 1; k < multiplicity; ++k) {
+                if (fabs((*this)(i, j, k, hws)) > fabs((*this)(i, j, k0, hws))) {
+                    k0 = k;
+                }
+            }
+
+            if ((*this)(i, j, k0, hws) < -EPS) {
+                for (int i2 = i; i2 < factor1_dimension; ++i2) {
+                    for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) {
+                        set(i2, j2, k0, hws, -(*this)(i2, j2, k0, hws));
+                    }
+                }
+            } else if ((*this)(i, j, k0, hws) < EPS) {
+                continue;
+            }
+
+            if (k0 != h) {
+                for (int i2 = i; i2 < factor1_dimension; ++i2) {
+                    for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) {
+                        double x = (*this)(i2, j2, k0, hws);
+                        set(i2, j2, k0, hws, (*this)(i2, j2, h, hws));
+                        set(i2, j2, h, hws, x);
+                    }
+                }
+            }
+
+            for (int k = h + 1; k < multiplicity; ++k) {
+                for (int i2 = i; i2 < factor1_dimension; ++i2) {
+                    for (int j2 = i2 == i ? j : 0; j2 < factor2_dimension; ++j2) {
+                        set(i2, j2, k, hws, (*this)(i2, j2, k, hws) - (*this)(i2, j2, h, hws)
+                                * (*this)(i, j, k, hws) / (*this)(i, j, h, hws));
+                    }
+                }
+            }
+
+            // next 3 lines not strictly necessary, might improve numerical stability
+            for (int k = h + 1; k < multiplicity; ++k) {
+                set(i, j, k, hws, 0.0);
+            }
+
+            ++h;
+        }
+    }
+
+    // Gram-Schmidt orthonormalization
+    for (int h = 0; h < multiplicity; ++h) {
+        for (int k = 0; k < h; ++k) {
+            double overlap = 0.0;
+            for (int i = 0; i < factor1_dimension; ++i) {
+                for (int j = 0; j < factor2_dimension; ++j) {
+                    overlap += (*this)(i, j, h, hws) * (*this)(i, j, k, hws);
+                }
+            }
+
+            for (int i = 0; i < factor1_dimension; ++i) {
+                for (int j = 0; j < factor2_dimension; ++j) {
+                    set(i, j, h, hws, (*this)(i, j, h, hws) - overlap * (*this)(i, j, k, hws));
+                }
+            }
+        }
+
+        double norm = 0.0;
+        for (int i = 0; i < factor1_dimension; ++i) {
+            for (int j = 0; j < factor2_dimension; ++j) {
+                norm += (*this)(i, j, h, hws) * (*this)(i, j, h, hws);
+            }
+        }
+        norm = std::sqrt(norm);
+
+        for (int i = 0; i < factor1_dimension; ++i) {
+            for (int j = 0; j < factor2_dimension; ++j) {
+                set(i, j, h, hws, (*this)(i, j, h, hws) / norm);
+            }
+        }
+    }
+}
+
+void clebsch::coefficients::compute_highest_weight_coeffs() {
+    if (multiplicity == 0) {
+        return;
+    }
+
+    std::vector<std::vector<int> > map_coeff(factor1_dimension,
+                                             std::vector<int>(factor2_dimension, -1));
+    std::vector<std::vector<int> > map_states(factor1_dimension,
+                                              std::vector<int>(factor2_dimension, -1));
+    int n_coeff = 0, n_states = 0;
+    pattern p(factor1, 0);
+
+    for (int i = 0; i < factor1_dimension; ++i, ++p) {
+        weight pw(p.get_weight());
+        pattern q(factor2, 0);
+        for (int j = 0; j < factor2_dimension; ++j, ++q) {
+            if (pw + q.get_weight() == irrep) {
+                map_coeff[i][j] = n_coeff++;
+            }
+        }
+    }
+
+    if (n_coeff == 1) {
+        for (int i = 0; i < factor1_dimension; ++i) {
+            for (int j = 0; j < factor2_dimension; ++j) {
+                if (map_coeff[i][j] >= 0) {
+                    set(i, j, 0, irrep_dimension - 1, 1.0);
+                    return;
+                }
+            }
+        }
+    }
+
+    double *hw_system = new double[n_coeff * (factor1_dimension * factor2_dimension)];
+    assert(hw_system != NULL);
+    memset(hw_system, 0, n_coeff * (factor1_dimension * factor2_dimension) * sizeof (double));
+
+    pattern r(factor1, 0);
+    for (int i = 0; i < factor1_dimension; ++i, ++r) {
+        pattern q(factor2, 0);
+
+        for (int j = 0; j < factor2_dimension; ++j, ++q) {
+            if (map_coeff[i][j] >= 0) {
+                for (int l = 1; l <= N - 1; ++l) {
+                    for (int k = 1; k <= l; ++k) {
+                        if ((k == 1 || r(k, l) + 1 <= r(k - 1, l - 1)) && r(k, l) + 1 <= r(k, l + 1)) {
+                            ++r(k, l);
+                            int h = r.index();
+                            --r(k, l);
+
+                            if (map_states[h][j] < 0) {
+                                map_states[h][j] = n_states++;
+                            }
+
+                            hw_system[n_coeff * map_states[h][j] + map_coeff[i][j]]
+                                += r.raising_coeff(k, l);
+                        }
+
+                        if ((k == 1 || q(k, l) + 1 <= q(k - 1, l - 1)) && q(k, l) + 1 <= q(k, l + 1)) {
+                            ++q(k, l);
+                            int h = q.index();
+                            --q(k, l);
+
+                            if (map_states[i][h] < 0) {
+                                map_states[i][h] = n_states++;
+                            }
+
+
+                            hw_system[n_coeff * map_states[i][h] + map_coeff[i][j]]
+                                += q.raising_coeff(k, l);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    int lwork = -1, info;
+    double worksize;
+
+    double *singval = new double[std::min(n_coeff, n_states)];
+    assert(singval != NULL);
+    double *singvec = new double[n_coeff * n_coeff];
+    assert(singvec != NULL);
+
+    dgesvd_("A",
+            "N",
+            &n_coeff,
+            &n_states,
+            hw_system,
+            &n_coeff,
+            singval,
+            singvec,
+            &n_coeff,
+            NULL,
+            &n_states,
+            &worksize,
+            &lwork,
+            &info);
+    assert(info == 0);
+
+    lwork = worksize;
+    double *work = new double[lwork];
+    assert(work != NULL);
+
+    dgesvd_("A",
+            "N",
+            &n_coeff,
+            &n_states,
+            hw_system,
+            &n_coeff,
+            singval,
+            singvec,
+            &n_coeff,
+            NULL,
+            &n_states,
+            work,
+            &lwork,
+            &info);
+    assert(info == 0);
+
+    for (int i = 0; i < multiplicity; ++i) {
+        for (int j = 0; j < factor1_dimension; ++j) {
+            for (int k = 0; k < factor2_dimension; ++k) {
+                if (map_coeff[j][k] >= 0) {
+                    double x = singvec[n_coeff * (n_coeff - 1 - i) + map_coeff[j][k]];
+
+                    if (fabs(x) > EPS) {
+                        set(j, k, i, irrep_dimension - 1, x);
+                    }
+                }
+            }
+        }
+    }
+
+    // uncomment next line to bring highest-weight coefficients into "normal form"
+    // highest_weight_normal_form();
+
+    delete[] work;
+    delete[] singvec;
+    delete[] singval;
+    delete[] hw_system;
+}
+
+void clebsch::coefficients::compute_lower_weight_coeffs(int multip_index,
+                                                        int state,
+                                                        std::vector<char> &done) {
+    weight statew(pattern(irrep, state).get_weight());
+    pattern p(irrep, 0);
+    std::vector<int> map_parent(irrep_dimension, -1),
+                     map_multi(irrep_dimension, -1),
+                     which_l(irrep_dimension, -1);
+    int n_parent = 0, n_multi = 0;
+
+    for (int i = 0; i < irrep_dimension; ++i, ++p) {
+        weight v(p.get_weight());
+
+        if (v == statew) {
+            map_multi[i] = n_multi++;
+        } else for (int l = 1; l < N; ++l) {
+            --v(l);
+            ++v(l + 1);
+            if (v == statew) {
+                map_parent[i] = n_parent++;
+                which_l[i] = l;
+                if (!done[i]) {
+                    compute_lower_weight_coeffs(multip_index, i, done);
+                }
+                break;
+            }
+            --v(l + 1);
+            ++v(l);
+        }
+    }
+
+    double *irrep_coeffs = new double[n_parent * n_multi];
+    assert(irrep_coeffs != NULL);
+    memset(irrep_coeffs, 0, n_parent * n_multi * sizeof (double));
+
+    double *prod_coeffs = new double[n_parent * factor1_dimension * factor2_dimension];
+    assert(prod_coeffs != NULL);
+    memset(prod_coeffs, 0, n_parent * factor1_dimension * factor2_dimension * sizeof (double));
+
+    std::vector<std::vector<int> > map_prodstat(factor1_dimension,
+                                                std::vector<int>(factor2_dimension, -1));
+    int n_prodstat = 0;
+
+    pattern r(irrep, 0);
+    for (int i = 0; i < irrep_dimension; ++i, ++r) {
+        if (map_parent[i] >= 0) {
+            for (int k = 1, l = which_l[i]; k <= l; ++k) {
+                if (r(k, l) > r(k + 1, l + 1) && (k == l || r(k, l) > r(k, l - 1))) {
+                    --r(k, l);
+                    int h = r.index();
+                    ++r(k, l);
+
+                    irrep_coeffs[n_parent * map_multi[h] + map_parent[i]] += r.lowering_coeff(k, l);
+                }
+            }
+
+            pattern q1(factor1, 0);
+            for (int j1 = 0; j1 < factor1_dimension; ++j1, ++q1) {
+                pattern q2(factor2, 0);
+
+                for (int j2 = 0; j2 < factor2_dimension; ++j2, ++q2) {
+                    if (std::fabs((*this)(j1, j2, multip_index, i)) > EPS) {
+                        for (int k = 1, l = which_l[i]; k <= l; ++k) {
+                            if (q1(k, l) > q1(k + 1, l + 1) && (k == l || q1(k, l) > q1(k, l - 1))) {
+                                --q1(k, l);
+                                int h = q1.index();
+                                ++q1(k, l);
+
+                                if (map_prodstat[h][j2] < 0) {
+                                    map_prodstat[h][j2] = n_prodstat++;
+                                }
+
+                                prod_coeffs[n_parent * map_prodstat[h][j2] + map_parent[i]] +=
+                                        (*this)(j1, j2, multip_index, i) * q1.lowering_coeff(k, l);
+                            }
+
+                            if (q2(k, l) > q2(k + 1, l + 1) && (k == l || q2(k, l) > q2(k, l - 1))) {
+                                --q2(k, l);
+                                int h = q2.index();
+                                ++q2(k, l);
+
+                                if (map_prodstat[j1][h] < 0) {
+                                    map_prodstat[j1][h] = n_prodstat++;
+                                }
+
+                                prod_coeffs[n_parent * map_prodstat[j1][h] + map_parent[i]] +=
+                                        (*this)(j1, j2, multip_index, i) * q2.lowering_coeff(k, l);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    double worksize;
+    int lwork = -1, info;
+
+    dgels_("N",
+           &n_parent,
+           &n_multi,
+           &n_prodstat,
+           irrep_coeffs,
+           &n_parent,
+           prod_coeffs,
+           &n_parent,
+           &worksize,
+           &lwork,
+           &info);
+    assert(info == 0);
+
+    lwork = worksize;
+    double *work = new double[lwork];
+    assert(work != NULL);
+
+    dgels_("N",
+           &n_parent,
+           &n_multi,
+           &n_prodstat,
+           irrep_coeffs,
+           &n_parent,
+           prod_coeffs,
+           &n_parent,
+           work,
+           &lwork,
+           &info);
+    assert(info == 0);
+
+    for (int i = 0; i < irrep_dimension; ++i) {
+        if (map_multi[i] >= 0) {
+            for (int j = 0; j < factor1_dimension; ++j) {
+                for (int k = 0; k < factor2_dimension; ++k) {
+                    if (map_prodstat[j][k] >= 0) {
+                        double x = prod_coeffs[n_parent * map_prodstat[j][k] + map_multi[i]];
+
+                        if (fabs(x) > EPS) {
+                            set(j, k, multip_index, i, x);
+                        }
+                    }
+                }
+            }
+
+            done[i] = true;
+        }
+    }
+
+    delete[] work;
+    delete[] prod_coeffs;
+    delete[] irrep_coeffs;
+}
+
+clebsch::coefficients::coefficients(const weight &irrep, const weight &factor1, const weight &factor2) :
+        N(irrep.N),
+        factor1(factor1),
+        factor2(factor2),
+        irrep(irrep),
+        factor1_dimension(factor1.dimension()),
+        factor2_dimension(factor2.dimension()),
+        irrep_dimension(irrep.dimension()),
+        multiplicity(decomposition(factor1, factor2).multiplicity(irrep)) {
+    assert(factor1.N == irrep.N);
+    assert(factor2.N == irrep.N);
+
+    compute_highest_weight_coeffs();
+
+    for (int i = 0; i < multiplicity; ++i) {
+        std::vector<char> done(irrep_dimension, 0);
+        done[irrep_dimension - 1] = true;
+        for (int j = irrep_dimension - 1; j >= 0; --j) {
+            if (!done[j]) {
+                compute_lower_weight_coeffs(i, j, done);
+            }
+        }
+    }
+}
+
+double clebsch::coefficients::operator()(int factor1_state,
+                                         int factor2_state,
+                                         int multiplicity_index,
+                                         int irrep_state) const {
+
+  assert(0 <= factor1_state && factor1_state < factor1_dimension);
+  assert(0 <= factor2_state && factor2_state < factor2_dimension);
+  assert(0 <= multiplicity_index && multiplicity_index < multiplicity);
+  assert(0 <= irrep_state && irrep_state < irrep_dimension);
+   
+  int coefficient_label[] = { factor1_state,
+                             factor2_state,
+                             multiplicity_index,
+                             irrep_state };
+  std::map<std::vector<int>, double>::const_iterator it(
+            clzx.find(std::vector<int>(coefficient_label, coefficient_label
+                                      + sizeof coefficient_label / sizeof coefficient_label[0])));
+  
+  return it != clzx.end() ? it->second : 0.0;
+}
+

diff --git a/clebsch.h b/clebsch.h
new file mode 100644 (file)
index 0000000..21fe4ce
--- /dev/null
+++ b/clebsch.h
@@ -0,0 +1,236 @@
+#ifndef CLEBSCH_H
+#define CLEBSCH_H
+
+//#define NDEBUG
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <functional>
+#include <fstream>
+#include <iostream>
+#include <map>
+#include <numeric>
+#include <vector>
+
+// Declaration of LAPACK subroutines
+// Make sure the data types match your version of LAPACK
+
+extern "C" void dgesvd_(char const* JOBU,
+                        char const* JOBVT,
+                        int const* M,
+                        int const* N,
+                        double* A,
+                        int const* LDA,
+                        double* S,
+                        double* U,
+                        int const* LDU,
+                        double* VT,
+                        int const* LDVT,
+                        double* WORK,
+                        int const* LWORK,
+                        int *INFO);
+
+extern "C" void dgels_(char const* TRANS,
+                       int const* M,
+                       int const* N,
+                       int const* NRHS,
+                       double* A,
+                       int const* LDA,
+                       double* B,
+                       int const* LDB,
+                       double* WORK,
+                       int const* LWORK,
+                       int *INFO);
+
+namespace clebsch {
+    const double EPS = 1e-12;
+
+    // binomial coefficients
+    class binomial_t {
+        std::vector<int> cache;
+        int N;
+
+    public:
+        int operator()(int n, int k);
+    };
+
+    // Eq. (19) and (25)
+    class weight {
+        std::vector<int> elem;
+
+    public:
+        // the N in "SU(N)"
+        const int N;
+
+        // create a non-initialized weight
+        weight(int N);
+
+        // create irrep weight of given index
+        // Eq. (C2)
+        weight(int N, int index);
+
+        // assign from another instance
+        clebsch::weight &operator=(const clebsch::weight &w);
+
+        // access elements of this weight (k = 1, ..., N)
+        int &operator()(int k);
+        const int &operator()(int k) const;
+
+        // compare weights
+        // Eq. (C1)
+        bool operator<(const weight &w) const;
+        bool operator==(const weight &w) const;
+
+        // element-wise sum of weights
+        clebsch::weight operator+(const weight &w) const;
+
+        // returns the index of this irrep weight (index = 0, 1, ...)
+        // Eq. (C2)
+        int index() const;
+
+        // returns the dimension of this irrep weight
+        // Eq. (22)
+        long long dimension() const;
+    };
+
+    // Eq. (20)
+    class pattern {
+        std::vector<int> elem;
+
+    public:
+        // the N in "SU(N)"
+        const int N;
+
+        // copy constructor
+        pattern(const pattern &pat);
+
+        // create pattern of given index from irrep weight
+        // Eq. (C7)
+        pattern(const weight &irrep, int index = 0);
+
+        // access elements of this pattern (l = 1, ..., N; k = 1, ..., l)
+        int &operator()(int k, int l);
+        const int &operator()(int k, int l) const;
+
+        // find succeeding/preceding pattern, return false if not possible
+        // Eq. (C9)
+        bool operator++();
+        bool operator--();
+
+        // returns the pattern index (index = 0, ..., dimension - 1)
+        // Eq. (C7)
+        int index() const;
+
+        // returns the pattern weight
+        // Eq. (25)
+        clebsch::weight get_weight() const;
+
+        // returns matrix element of lowering operator J^(l)_-
+        // between this pattern minus M^(k,l) and this pattern
+        // (l = 1, ..., N; k = 1, ..., l)
+        // Eq. (28)
+        double lowering_coeff(int k, int l) const;
+
+        // returns matrix element of raising operator J^(l)_+
+        // between this pattern plus M^(k,l) and this pattern
+        // (l = 1, ..., N; k = 1, ..., l)
+        // Eq. (29)
+        double raising_coeff(int k, int l) const;
+    };
+
+    class decomposition {
+        std::vector<clebsch::weight> weights;
+        std::vector<int> multiplicities;
+
+    public:
+        // the N in "SU(N)"
+        const int N;
+
+        // save given irreps for later use
+        const weight factor1, factor2;
+
+        // construct the decomposition of factor1 times factor2 into irreps
+        // Eq. (31)
+        decomposition(const weight &factor1, const weight &factor2);
+
+        // return the number of occurring irreps
+        int size() const;
+
+        // access the occurring irreps
+        // j = 0, ..., size() - 1
+        const clebsch::weight &operator()(int j) const;
+
+        // return the outer multiplicity of irrep in this decomposition
+        int multiplicity(const weight &irrep) const;
+    };
+
+    class index_adapter {
+        std::vector<int> indices;
+        std::vector<int> multiplicities;
+
+    public:
+        // the N in "SU(N)"
+        const int N;
+
+        // save given irreps for later use
+        const int factor1, factor2;
+
+        // construct this index_adapter from a given decomposition
+        index_adapter(const clebsch::decomposition &decomp);
+
+        // return the number of occurring irreps
+        int size() const;
+
+        // access the occurring irreps
+        int operator()(int j) const;
+
+        // return the outer multiplicity of irrep in this decomposition
+        int multiplicity(int irrep) const;
+    };
+
+    class coefficients {
+        std::map<std::vector<int>, double> clzx;
+
+        // access Clebsch-Gordan coefficients in convenient manner
+        void set(int factor1_state,
+                 int factor2_state,
+                 int multiplicity_index,
+                 int irrep_state,
+                 double value);
+
+        // internal functions, doing most of the work
+        void highest_weight_normal_form(); // Eq. (37)
+        void compute_highest_weight_coeffs(); // Eq. (36)
+        void compute_lower_weight_coeffs(int multip_index, int state, std::vector<char> &done); // Eq. (40)
+
+    public:
+        // the N in "SU(N)"
+        const int N;
+
+        // save irreps and their dimensions for later use
+        const weight factor1, factor2, irrep;
+        const int factor1_dimension, factor2_dimension, irrep_dimension;
+
+        // outer multiplicity of irrep in this decomposition
+        const int multiplicity;
+
+        // construct all Clebsch-Gordan coefficients of this decomposition
+        coefficients(const weight &irrep, const weight &factor1, const weight &factor2);
+
+        // access Clebsch-Gordan coefficients (read-only)
+        // multiplicity_index = 0, ..., multiplicity - 1
+        // factor1_state = 0, ..., factor1_dimension - 1
+        // factor2_state = 0, ..., factor2_dimension - 1
+        // irrep_state = 0, ..., irrep_dimension
+        double operator()(int factor1_state,
+                          int factor2_state,
+                          int multiplicity_index,
+                          int irrep_state) const;
+    };
+};
+
+#endif

diff --git a/clebschgordan.cpp b/clebschgordan.cpp
new file mode 100644 (file)
index 0000000..81a37bd
--- /dev/null
+++ b/clebschgordan.cpp
@@ -0,0 +1,251 @@
+#include <iostream>
+
+#include "clebsch.h"
+
+using namespace std;
+
+int main() {
+    while (true) {
+        int choice, N;
+
+        cout << "What would you like to do?" << endl;
+        cout << "1) Translate an i-weight S to its index P(S)" << endl;
+        cout << "2) Recover an i-weight S from its index P(S)" << endl;
+        cout << "3) Translate a pattern M to its index Q(M)" << endl;
+        cout << "4) Recover a pattern M from its index Q(M)" << endl;
+        cout << "5) Calculate Clebsch-Gordan coefficients for S x S' -> S''" << endl;
+        cout << "6) Calculate all Glebsch-Gordan coefficients for S x S'" << endl;
+        cout << "0) Quit" << endl;
+
+        do {
+            cin >> choice;
+        } while (choice < 0 || choice > 6);
+
+        if (choice == 0) {
+            break;
+        }
+
+        cout << "N (e.g. 3): ";
+        cin >> N;
+
+        switch (choice) {
+            case 1: {
+                clebsch::weight S(N);
+                cout << "Irrep S: ";
+                for (int k = 1; k <= N; ++k) {
+                    cin >> S(k);
+                }
+                cout << S.index() << endl;
+                break;
+            }
+            case 2: {
+                int P;
+                cout << "Index: ";
+                cin >> P;
+                clebsch::weight S(N, P);
+                cout << "I-weight:";
+                for (int k = 1; k <= N; ++k) {
+                    cout << ' ' << S(k);
+                }
+                cout << endl;
+                break;
+            }
+            case 3: {
+                clebsch::pattern M(N);
+                for (int l = N; l >= 1; --l) {
+                    cout << "Row l = " << l << ": ";
+                    for (int k = 1; k <= l; ++k) {
+                        cin >> M(k, l);
+                    }
+                }
+                cout << "Index: " << M.index() + 1 << endl;
+                break;
+            }
+            case 4: {
+                clebsch::weight S(N);
+                cout << "Irrep S: ";
+                for (int i = 1; i <= N; ++i) {
+                    cin >> S(i);
+                }
+
+                int Q;
+                cout << "Index (1..dim(S)): ";
+                cin >> Q;
+                clebsch::pattern M(S, Q - 1);
+                for (int l = N; l >= 1; --l) {
+                    for (int k = 1; k <= l; ++k) {
+                        cout << M(k, l) << '\t';
+                    }
+                    cout << endl;
+                }
+                break;
+            }
+            case 5: {
+                clebsch::weight S(N);
+                cout << "Irrep S (e.g.";
+                for (int k = N - 1; k >= 0; --k) {
+                    cout << ' ' << k;
+                }
+                cout << "): ";
+                for (int k = 1; k <= N; ++k) {
+                    cin >> S(k);
+                }
+
+                clebsch::weight Sprime(N);
+                cout << "Irrep S' (e.g.";
+                for (int k = N - 1; k >= 0; --k) {
+                    cout << ' ' << k;
+                }
+                cout << "): ";
+                for (int k = 1; k <= N; ++k) {
+                    cin >> Sprime(k);
+                }
+
+                clebsch::decomposition decomp(S, Sprime);
+                cout << "Littlewood-Richardson decomposition S \\otimes S' = \\oplus S'':" << endl;
+                cout << "[irrep index] S'' (outer multiplicity) {dimension d_S}" << endl;
+                for (int i = 0; i < decomp.size(); ++i) {
+                    cout << "[" << decomp(i).index() << "] ";
+                    for (int k = 1; k <= N; ++k) {
+                        cout << decomp(i)(k) << ' ';
+                    }
+                    cout << '(' << decomp.multiplicity(decomp(i)) << ") {"
+                         << decomp(i).dimension() << "}" << endl;;
+                }
+
+                clebsch::weight Sdoubleprime(N);
+                for (bool b = true; b; ) {
+                    cout << "Irrep S'': ";
+                    for (int k = 1; k <= N; ++k) {
+                        cin >> Sdoubleprime(k);
+                    }
+                    for (int i = 0; i < decomp.size(); ++i) {
+                        if (decomp(i) == Sdoubleprime) {
+                            b = false;
+                            break;
+                        }
+                    }
+                    if (b) {
+                        cout << "S'' does not occur in the decomposition" << endl;
+                    }
+                }
+
+                int alpha;
+                while (true) {
+                    cout << "Outer multiplicity index: ";
+                    cin >> alpha;
+                    if (1 <= alpha && alpha <= decomp.multiplicity(Sdoubleprime)) {
+                        break;
+                    }
+                    cout << "S'' does not have this multiplicity" << endl;
+                }
+
+                string file_name;
+                cout << "Enter file name to write to file (leave blank for screen output): ";
+                cin.ignore(1234, '\n');
+                getline(cin, file_name);
+
+                const clebsch::coefficients C(Sdoubleprime, S, Sprime);
+                int dimS = S.dimension(),
+                    dimSprime = Sprime.dimension(),
+                    dimSdoubleprime = Sdoubleprime.dimension();
+
+                ofstream os(file_name.c_str());
+                (file_name.empty() ? cout : os).setf(ios::fixed);
+                (file_name.empty() ? cout : os).precision(15);
+                (file_name.empty() ? cout : os) << "List of nonzero CGCs for S x S' => S'', alpha" << endl;
+                (file_name.empty() ? cout : os) << "Q(M)\tQ(M')\tQ(M'')\tCGC" << endl;
+                for (int i = 0; i < dimSdoubleprime; ++i) {
+                    for (int j = 0; j < dimS; ++j) {
+                        for (int k = 0; k < dimSprime; ++k) {
+                            double x = double(C(j, k, alpha - 1, i));
+
+                            if (fabs(x) > clebsch::EPS) {
+                                (file_name.empty() ? cout : os) << j + 1 << '\t'
+                                    << k + 1 << '\t' << i + 1 << '\t' << x << endl;
+                            }
+                        }
+                    }
+                }
+
+                break;
+            }
+            case 6: {
+                clebsch::weight S(N);
+                cout << "Irrep S (e.g.";
+                for (int k = N - 1; k >= 0; --k) {
+                    cout << ' ' << k;
+                }
+                cout << "): ";
+                for (int k = 1; k <= N; ++k) {
+                    cin >> S(k);
+                }
+
+                clebsch::weight Sprime(N);
+                cout << "Irrep S' (e.g.";
+                for (int k = N - 1; k >= 0; --k) {
+                    cout << ' ' << k;
+                }
+                cout << "): ";
+                for (int k = 1; k <= N; ++k) {
+                    cin >> Sprime(k);
+                }
+
+                string file_name;
+                cout << "Enter file name to write to file (leave blank for screen output): ";
+                cin.ignore(1234, '\n');
+                getline(cin, file_name);
+
+                ofstream os(file_name.c_str());
+                (file_name.empty() ? cout : os).setf(ios::fixed);
+                (file_name.empty() ? cout : os).precision(15);
+
+                clebsch::decomposition decomp(S, Sprime);
+                (file_name.empty() ? cout : os) <<
+                    "Littlewood-Richardson decomposition S \\otimes S' = \\oplus S'':" << endl;
+                (file_name.empty() ? cout : os) <<
+                    "[irrep index] S'' (outer multiplicity) {dimension d_S}" << endl;
+                for (int i = 0; i < decomp.size(); ++i) {
+                    (file_name.empty() ? cout : os) << "[" << decomp(i).index() << "] ";
+                    for (int k = 1; k <= N; ++k) {
+                        (file_name.empty() ? cout : os) << decomp(i)(k) << ' ';
+                    }
+                    (file_name.empty() ? cout : os) << '(' << decomp.multiplicity(decomp(i)) << ") {"
+                         << decomp(i).dimension() << "}" << endl;;
+                }
+
+                for (int i = 0; i < decomp.size(); ++i) {
+                    const clebsch::coefficients C(decomp(i),S, Sprime);
+                    int dimS = S.dimension(),
+                        dimSprime = Sprime.dimension(),
+                        dimSdoubleprime = decomp(i).dimension();
+
+                    for (int m = 0; m < C.multiplicity; ++m) {
+                        (file_name.empty() ? cout : os) << "List of nonzero CGCs for S x S' => S'' = (";
+                        for (int j = 1; j <= N; ++j) cout << decomp(i)(j) << (j < N ? ' ' : ')');
+                        (file_name.empty() ? cout : os) << ", alpha = " << m + 1 << endl;
+                        (file_name.empty() ? cout : os) << "Q(M)\tQ(M')\tQ(M'')\tCGC" << endl;
+                        for (int i = 0; i < dimSdoubleprime; ++i) {
+                            for (int j = 0; j < dimS; ++j) {
+                                for (int k = 0; k < dimSprime; ++k) {
+                                    double x = double(C(j, k, m, i));
+
+                                    if (fabs(x) > clebsch::EPS) {
+                                        (file_name.empty() ? cout : os) << j  + 1<< '\t'
+                                            << k + 1 << '\t' << i + 1 << '\t' << x << endl;
+                                    }
+                                }
+                            }
+                        }
+
+                        (file_name.empty() ? cout : os) << endl;
+                    }
+                }
+
+                break;
+            }
+        }
+    }
+
+    return 0;
+}

diff --git a/su2clebsch.cpp b/su2clebsch.cpp
new file mode 100644 (file)
index 0000000..49a8047
--- /dev/null
+++ b/su2clebsch.cpp
@@ -0,0 +1,60 @@
+#include <iostream>
+
+#include "clebsch.h"
+
+#include "su2clebsch.hpp"
+
+using namespace std;
+
+namespace su2clebsch {
+
+  double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
+             int irrep_2j, int irrep_2m) 
+  {
+    const int factor1_n = (factor1_2j + factor1_2m) / 2;
+    const int factor2_n = (factor2_2j + factor2_2m) / 2;
+    const int irrep_n =   (irrep_2j + irrep_2m) / 2;
+
+    clebsch::weight factor1_S(2);
+    factor1_S(1) = factor1_2j;
+    if ( ! (0 <= factor1_n && factor1_n < factor1_S.dimension()) ) {
+      //std::cerr << "clebsch warning: 0 <= factor1_state && factor1_state < factor1_dimension" << std::endl;
+      return 0.0;
+    }
+
+    clebsch::weight factor2_S(2);
+    factor2_S(1) = factor2_2j;
+    if ( ! (0 <= factor2_n && factor2_n < factor2_S.dimension()) ) {
+      //std::cerr << "clebsch warning: 0 <= factor2_state && factor2_state < factor2_dimension" << std::endl;
+      return 0.0;
+    }
+
+    clebsch::weight irrep_S(2);
+    irrep_S(1) = irrep_2j;
+    if ( ! (0 <= irrep_n && irrep_n < irrep_S.dimension()) ) {
+      //std::cerr << "clebsch warning: 0 <= irrep_state && irrep_state < irrep_dimension" << std::endl;
+      return 0.0;
+    }
+
+    clebsch::decomposition decomp(factor1_S, factor2_S);
+
+    bool b = true;
+    for (int i = 0; i < decomp.size(); ++i) {
+      if (decomp(i) == irrep_S) {
+       b = false;
+       break;
+      }
+    }
+    if (b) {
+      cerr << "WARNING: 2j=" << irrep_2j 
+          << " does not occour in decomposition!" << endl;
+      return 0.0;
+    }
+
+    const clebsch::coefficients C(irrep_S, factor1_S, factor2_S);
+
+    return C( factor1_n, factor2_n, 0, irrep_n );
+  }
+
+};
+

diff --git a/su2clebsch.hpp b/su2clebsch.hpp
new file mode 100644 (file)
index 0000000..689c0a1
--- /dev/null
+++ b/su2clebsch.hpp
@@ -0,0 +1,11 @@
+#ifndef SU2CLEBSCH_HPP
+#define SU2CLEBSCH_HPP
+
+namespace su2clebsch {
+
+  double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
+             int irrep_2j, int irrep_2m);
+
+};
+
+#endif