1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
#include "Cleanup.hh"
#include "algorithms/eliminate_kronecker.hh"
#include "properties/KroneckerDelta.hh"
#include "properties/Integer.hh"
using namespace cadabra;
eliminate_kronecker::eliminate_kronecker(const Kernel& k, Ex& tr)
: Algorithm(k, tr)
{
}
bool eliminate_kronecker::can_apply(iterator st)
{
if(*st->name!="\\prod")
if(!is_single_term(st))
return false;
return true;
}
Algorithm::result_t eliminate_kronecker::apply(iterator& st)
{
result_t ret=result_t::l_no_action;
prod_wrap_single_term(st);
const nset_t::iterator onept=name_set.insert("1").first;
int looping=0;
sibling_iterator it=tr.begin(st);
while(it!=tr.end(st)) { // Loop over all factors in a product, looking for Kroneckers
bool replaced=false;
// std::cerr << *it->name << std::endl;
const KroneckerDelta *kr=kernel.properties.get<KroneckerDelta>(it);
if(kr && number_of_indices(it)==2) {
// std::cerr << "KD " << Ex(it) << std::endl;
index_iterator ii1=begin_index(it);
index_iterator ii2=ii1;
++ii2;
if(subtree_compare(&kernel.properties, ii1, ii2, 1, false, -2, true)==0) { // a self-contracted Kronecker delta
// std::cerr << "self-contracted delta with " << Ex(ii1) << " = " << Ex(ii2) << std::endl;
const Integer *itg1=kernel.properties.get<Integer>(ii1, true);
const Integer *itg2=kernel.properties.get<Integer>(ii2, true);
if(itg1 && itg2) {
if(ii1->is_rational()==false && ii2->is_rational()==false) {
if(itg1->from.begin()!=itg1->from.end() && itg2->from.begin()!=itg2->from.end()) {
if(itg1->difference.begin()->name==onept) {
multiply(st->multiplier, *itg1->difference.begin()->multiplier);
it=tr.erase(it);
}
else {
it=tr.replace(it, itg1->difference.begin());
}
ret=result_t::l_applied;
}
else ++it;
}
else if(ii1->is_rational() && ii2->is_rational()) {
if(ii1->multiplier==ii2->multiplier) {
it=tr.erase(it);
}
else {
zero(it->multiplier);
}
}
}
else ++it;
}
else {
sibling_iterator oi=tr.begin(st);
++looping;
// iterate over all factors in the product
bool doing1=false;
bool doing2=false;
while(!replaced && oi!=tr.end(st)) {
if(oi!=it) { // this is not the delta node
// compare delta indices with all indices of this object
index_iterator ind=begin_index(oi);
while(ind!=end_index(oi)) {
index_iterator nxt=ind;
++nxt;
if(ii1->is_rational()==false && subtree_compare(&kernel.properties, ind, ii1, 1, false, -2, true)==0 ) {
if(! (replaced && doing2) ) {
multiplier_t mt=(*ind->multiplier) / (*ii1->multiplier);
iterator rep=tr.replace_index(ind, ii2);
rep->fl.parent_rel=ii2->fl.parent_rel;
multiply(rep->multiplier, mt);
replaced=true;
doing1=true;
}
// cannot 'break' here because that would miss cases when the
// delta multiplies a sum.
}
else if(ii2->is_rational()==false && subtree_compare(&kernel.properties, ind, ii2, 1, false, -2, true)==0) {
if(! (replaced && doing1) ) {
multiplier_t mt=(*ind->multiplier) / (*ii2->multiplier);
iterator rep=tr.replace_index(ind, ii1);
rep->fl.parent_rel=ii1->fl.parent_rel;
multiply(rep->multiplier, mt);
replaced=true;
doing2=true;
}
// no break here either.
}
ind=nxt;
}
}
if(!replaced)
++oi;
}
if(replaced) {
ret=result_t::l_applied;
iterator tmp=oi;
cleanup_dispatch(kernel, tr, tmp);
it=tr.erase(it);
}
else ++it;
}
}
else ++it;
}
// the product may have reduced to a single term or even just a constant
// txtout << "exiting eliminate" << std::endl;
// prod_unwrap_single_term(st);
// txtout << st.node << " " << tr.parent(st).node << std::endl;
// txtout << *st->name << " " << *(tr.parent(st)->name) << std::endl;
sibling_iterator ff=tr.begin(st);
if(ff==tr.end(st)) {
st->name=onept;
}
else {
++ff;
if(ff==tr.end(st)) {
tr.begin(st)->fl.bracket=st->fl.bracket;
tr.begin(st)->fl.parent_rel=st->fl.parent_rel;
tr.begin(st)->multiplier=st->multiplier;
tr.flatten(st);
st=tr.erase(st);
}
}
cleanup_dispatch(kernel, tr, st);
// cleanup_sums_products(tr, st);
// txtout << "looped " << looping << std::endl;
return ret;
}
|
