/** Clear some status_flags. */
const basic & clearflag(unsigned f) const {flags &= ~f; return *this;}
- /** Get relative precedence level (useful for implementing pretty-printed
- * output). */
- unsigned get_precedence(void) const {return precedence;}
-
protected:
void ensure_if_modifiable(void) const;
{
exvector s;
s.reserve(v.size());
- unsigned rl = ex_to_clifford(v[0]).get_representation_label();
// Remove superfluous ONEs
exvector::const_iterator cit = v.begin(), citend = v.end();
const ex & ib = b.op(1);
if (ia.is_equal(ib)) {
a = lorentz_g(ia, ib);
- b = dirac_ONE(rl);
+ b = dirac_ONE(representation_label);
something_changed = true;
}
}
}
if (s.size() == 0)
- return clifford(diracone(), rl) * sign;
+ return clifford(diracone(), representation_label) * sign;
if (something_changed)
return nonsimplified_ncmul(s) * sign;
else
v.push_back(e.op(i));
// Stupid bubble sort because we only want to swap adjacent gammas
- exvector::iterator itstart = v.begin(), itend = v.end(), next_to_last = itend - 1;
- if (is_ex_of_type(itstart->op(0), diracgamma5))
- itstart++;
- while (next_to_last != itstart) {
- exvector::iterator it = itstart;
- while (it != next_to_last) {
- if (it[0].op(1).compare(it[1].op(1)) > 0) {
- ex save0 = it[0], save1 = it[1];
- it[0] = lorentz_g(it[0].op(1), it[1].op(1));
- it[1] = _ex2();
- ex sum = ncmul(v);
- it[0] = save1;
- it[1] = save0;
- sum -= ncmul(v);
- return canonicalize_clifford(sum);
- }
- it++;
+ exvector::iterator it = v.begin(), next_to_last = v.end() - 1;
+ if (is_ex_of_type(it->op(0), diracgamma5))
+ it++;
+ while (it != next_to_last) {
+ if (it[0].op(1).compare(it[1].op(1)) > 0) {
+ ex save0 = it[0], save1 = it[1];
+ it[0] = lorentz_g(it[0].op(1), it[1].op(1));
+ it[1] = _ex2();
+ ex sum = ncmul(v);
+ it[0] = save1;
+ it[1] = save0;
+ sum -= ncmul(v);
+ return canonicalize_clifford(sum);
}
- next_to_last--;
+ it++;
}
return ncmul(v);
}
{
exvector s;
s.reserve(v.size());
- unsigned rl = ex_to_color(v[0]).get_representation_label();
// Remove superfluous ONEs
exvector::const_iterator it = v.begin(), itend = v.end();
}
if (s.size() == 0)
- return color(su3one(), rl);
+ return color(su3one(), representation_label);
else
return simplified_ncmul(s);
}
if (is_ex_exactly_of_type(other->op(0), su3d)) {
// Find the dummy indices of the contraction
- exvector dummy_indices;
- dummy_indices = ex_to_indexed(*self).get_dummy_indices(ex_to_indexed(*other));
+ exvector self_indices = ex_to_indexed(*self).get_indices();
+ exvector other_indices = ex_to_indexed(*other).get_indices();
+ exvector all_indices = self_indices;
+ all_indices.insert(all_indices.end(), other_indices.begin(), other_indices.end());
+ exvector free_indices, dummy_indices;
+ find_free_and_dummy(all_indices, free_indices, dummy_indices);
// d.abc*d.abc=40/3
if (dummy_indices.size() == 3) {
// d.akl*d.bkl=5/3*delta.ab
} else if (dummy_indices.size() == 2) {
- exvector a = index_set_difference(ex_to_indexed(*self).get_indices(), dummy_indices);
- exvector b = index_set_difference(ex_to_indexed(*other).get_indices(), dummy_indices);
- GINAC_ASSERT(a.size() > 0);
- GINAC_ASSERT(b.size() > 0);
- *self = numeric(5, 3) * delta_tensor(a[0], b[0]);
+ exvector a;
+ back_insert_iterator<exvector> ita(a);
+ ita = set_difference(self_indices.begin(), self_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less());
+ ita = set_difference(other_indices.begin(), other_indices.end(), dummy_indices.begin(), dummy_indices.end(), ita, ex_is_less());
+ GINAC_ASSERT(a.size() == 2);
+ *self = numeric(5, 3) * delta_tensor(a[0], a[1]);
*other = _ex1();
return true;
}
*/
#include <stdexcept>
+#include <algorithm>
#include "idx.h"
#include "symbol.h"
return is_dummy_pair(ex_to_idx(e1), ex_to_idx(e2));
}
-/** Bring a vector of indices into a canonic order. Dummy indices will lie
- * next to each other after the sorting. */
-static void sort_index_vector(exvector &v)
+// Shaker sort is sufficient for the expected small number of indices
+template <class It, class Cmp>
+inline void shaker_sort(It first, It last, Cmp comp)
{
- // Nothing to sort if less than 2 elements
- if (v.size() < 2)
+ if (first == last)
return;
-
- // Simple bubble sort algorithm should be sufficient for the small
- // number of indices expected
- exvector::iterator it1 = v.begin(), itend = v.end(), next_to_last_idx = itend - 1;
- while (it1 != next_to_last_idx) {
- exvector::iterator it2 = it1 + 1;
- while (it2 != itend) {
- if (it1->compare(*it2) > 0)
- it1->swap(*it2);
- it2++;
+ --last;
+ if (first == last)
+ return;
+ It flag = first;
+ do {
+ It i;
+ for (i=last; i>first; --i) {
+ if (comp(*i, i[-1])) {
+ iter_swap(i-1, i);
+ flag = i - 1;
+ }
}
- it1++;
- }
+ ++flag;
+ first = flag;
+ for (i=first; i<last; ++i) {
+ if (comp(i[1], *i)) {
+ iter_swap(i, i+1);
+ flag = i + 1;
+ }
+ }
+ last = flag - 1;
+ } while (first <= last);
}
-
void find_free_and_dummy(exvector::const_iterator it, exvector::const_iterator itend, exvector & out_free, exvector & out_dummy)
{
out_free.clear();
// Sort index vector. This will cause dummy indices come to lie next
// to each other (because the sort order is defined to guarantee this).
exvector v(it, itend);
- sort_index_vector(v);
+ shaker_sort(v.begin(), v.end(), ex_is_less());
// Find dummy pairs and free indices
it = v.begin(); itend = v.end();
out_free.push_back(*last);
}
-exvector index_set_difference(const exvector & set1, const exvector & set2)
-{
- exvector ret;
-
- exvector::const_iterator ait = set1.begin(), aitend = set1.end();
- while (ait != aitend) {
- exvector::const_iterator bit = set2.begin(), bitend = set2.end();
- bool found = false;
- while (bit != bitend) {
- if (ait->is_equal(*bit)) {
- found = true;
- break;
- }
- bit++;
- }
- if (!found)
- ret.push_back(*ait);
- ait++;
- }
-
- return ret;
-}
-
} // namespace GiNaC
return free_indices.size();
}
-/** Given two index vectors, find those indices that appear in the first
- * vector but not in the second one (asymmetric set difference). */
-exvector index_set_difference(const exvector & set1, const exvector & set2);
-
} // namespace GiNaC
#endif // ndef __GINAC_IDX_H__
return inherited::info(inf);
}
+struct idx_is_not : public binary_function<ex, unsigned, bool> {
+ bool operator() (const ex & e, unsigned inf) const {
+ return !(ex_to_idx(e).get_value().info(inf));
+ }
+};
+
bool indexed::all_index_values_are(unsigned inf) const
{
// No indices? Then no property can be fulfilled
return false;
// Check all indices
- exvector::const_iterator it = seq.begin() + 1, itend = seq.end();
- while (it != itend) {
- GINAC_ASSERT(is_ex_of_type(*it, idx));
- if (!ex_to_idx(*it).get_value().info(inf))
- return false;
- it++;
- }
- return true;
+ return find_if(seq.begin() + 1, seq.end(), bind2nd(idx_is_not(), inf)) == seq.end();
}
int indexed::compare_same_type(const basic & other) const
if (v1.size() != v2.size())
return false;
- // And also the indices themselves
- exvector::const_iterator ait = v1.begin(), aitend = v1.end(),
- bit = v2.begin(), bitend = v2.end();
- while (ait != aitend) {
- if (!ait->is_equal(*bit))
- return false;
- ait++; bit++;
- }
- return true;
+ return equal(v1.begin(), v1.end(), v2.begin(), ex_is_equal());
}
exvector indexed::get_indices(void) const
return basis.get_free_indices();
}
-/* Function object for STL sort() */
-struct ex_is_less {
- bool operator() (const ex &lh, const ex &rh) const
- {
- return lh.compare(rh) < 0;
- }
-};
-
/** Rename dummy indices in an expression.
*
* @param e Expression to be worked on
if (local_size == 0)
return e;
- sort(local_dummy_indices.begin(), local_dummy_indices.end(), ex_is_less());
-
if (global_size < local_size) {
// More local indices than we encountered before, add the new ones
int remaining = local_size - global_size;
exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
while (it != itend && remaining > 0) {
- exvector::const_iterator git = global_dummy_indices.begin(), gitend = global_dummy_indices.end();
- while (git != gitend) {
- if (it->is_equal(*git))
- goto found;
- git++;
+ if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(ex_is_equal(), *it)) == global_dummy_indices.end()) {
+ global_dummy_indices.push_back(*it);
+ global_size++;
+ remaining--;
}
- global_dummy_indices.push_back(*it);
- global_size++;
- remaining--;
-found: it++;
+ it++;
}
- sort(global_dummy_indices.begin(), global_dummy_indices.end(), ex_is_less());
}
// Replace index symbols in expression
for (unsigned i=0; i<local_size; i++) {
ex loc_sym = local_dummy_indices[i].op(0);
ex glob_sym = global_dummy_indices[i].op(0);
- if (!loc_sym.is_equal(glob_sym))
+ if (!loc_sym.is_equal(glob_sym)) {
all_equal = false;
- local_syms.append(loc_sym);
- global_syms.append(glob_sym);
+ local_syms.append(loc_sym);
+ global_syms.append(glob_sym);
+ }
}
if (all_equal)
return e;
}
// Find free indices (concatenate them all and call find_free_and_dummy())
- exvector un, local_dummy_indices;
+ // and all dummy indices that appear
+ exvector un, individual_dummy_indices;
it1 = v.begin(); itend = v.end();
while (it1 != itend) {
- exvector free_indices_of_factor = it1->get_free_indices();
+ exvector free_indices_of_factor;
+ if (is_ex_of_type(*it1, indexed)) {
+ exvector dummy_indices_of_factor;
+ find_free_and_dummy(ex_to_indexed(*it1).seq.begin() + 1, ex_to_indexed(*it1).seq.end(), free_indices_of_factor, dummy_indices_of_factor);
+ individual_dummy_indices.insert(individual_dummy_indices.end(), dummy_indices_of_factor.begin(), dummy_indices_of_factor.end());
+ } else
+ free_indices_of_factor = it1->get_free_indices();
un.insert(un.end(), free_indices_of_factor.begin(), free_indices_of_factor.end());
it1++;
}
+ exvector local_dummy_indices;
find_free_and_dummy(un, free_indices, local_dummy_indices);
+ local_dummy_indices.insert(local_dummy_indices.end(), individual_dummy_indices.begin(), individual_dummy_indices.end());
ex r;
if (something_changed)
ex e_expanded = e.expand();
// Simplification of single indexed object: just find the free indices
- // (and perform dummy index renaming if
+ // and perform dummy index renaming
if (is_ex_of_type(e_expanded, indexed)) {
const indexed &i = ex_to_indexed(e_expanded);
exvector local_dummy_indices;
unsigned row; ///< number of rows
unsigned col; ///< number of columns
exvector m; ///< representation (cols indexed first)
- static unsigned precedence;
};
cout << "k_cum[" << i << "]=" << k_cum[i] << endl;
cout << "upper_limit[" << i << "]=" << upper_limit[i] << endl;
}
- for (exvector::const_iterator cit=term.begin(); cit!=term.end(); ++cit) {
- cout << *cit << endl;
- }
+ for_each(term.begin(), term.end(), ostream_iterator<ex>(cout, "\n"));
cout << "end term" << endl;
*/
classname(const classname & other); \
const classname & operator=(const classname & other); \
basic * duplicate() const; \
+ unsigned get_precedence(void) const {return precedence;} \
protected: \
void copy(const classname & other); \
void destroy(bool call_parent); \
}
#endif
-/** Append one exvector to another */
-void append_exvector_to_exvector(exvector & dest, const exvector & source)
-{
- dest.reserve(dest.size() + source.size());
- dest.insert(dest.end(), source.begin(), source.end());
-}
-
//////////
// `construct on first use' chest of numbers
//////////
#include <string>
#include <stdexcept>
+#include <functional>
#if defined(HAVE_SSTREAM)
#include <sstream>
#elif defined(HAVE_STRSTREAM)
return sigma;
}
-void append_exvector_to_exvector(exvector & dest, const exvector & source);
+/* Function objects for STL sort() etc. */
+struct ex_is_less : public binary_function<ex, ex, bool> {
+ bool operator() (const ex &lh, const ex &rh) const { return lh.compare(rh) < 0; }
+};
+
+struct ex_is_equal : public binary_function<ex, ex, bool> {
+ bool operator() (const ex &lh, const ex &rh) const { return lh.is_equal(rh); }
+};
// Collection of `construct on first use' wrappers for safely avoiding
// internal object replication without running into the `static
git://www.ginac.de / ginac.git / commitdiff