+exvector get_all_dummy_indices_safely(const ex & e)
+{
+ if (is_a<indexed>(e))
+ return ex_to<indexed>(e).get_dummy_indices();
+ else if (is_a<power>(e) && e.op(1)==2) {
+ return e.op(0).get_free_indices();
+ }
+ else if (is_a<mul>(e) || is_a<ncmul>(e)) {
+ exvector dummies;
+ exvector free_indices;
+ for (int i=0; i<e.nops(); ++i) {
+ exvector dummies_of_factor = get_all_dummy_indices_safely(e.op(i));
+ dummies.insert(dummies.end(), dummies_of_factor.begin(),
+ dummies_of_factor.end());
+ exvector free_of_factor = e.op(i).get_free_indices();
+ free_indices.insert(free_indices.begin(), free_of_factor.begin(),
+ free_of_factor.end());
+ }
+ exvector free_out, dummy_out;
+ find_free_and_dummy(free_indices.begin(), free_indices.end(), free_out,
+ dummy_out);
+ dummies.insert(dummies.end(), dummy_out.begin(), dummy_out.end());
+ return dummies;
+ }
+ else if(is_a<add>(e)) {
+ exvector result;
+ for(int i=0; i<e.nops(); ++i) {
+ exvector dummies_of_term = get_all_dummy_indices_safely(e.op(i));
+ sort(dummies_of_term.begin(), dummies_of_term.end());
+ exvector new_vec;
+ set_union(result.begin(), result.end(), dummies_of_term.begin(),
+ dummies_of_term.end(), std::back_inserter<exvector>(new_vec),
+ ex_is_less());
+ result.swap(new_vec);
+ }
+ return result;
+ }
+ return exvector();
+}
+