--- /dev/null
+#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;
+}
+
--- /dev/null
+#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
--- /dev/null
+#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;
+}
projects / phys / su2clebsch.git / commitdiff