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
|
#include "Exceptions.hh"
#include "algorithms/reduce_delta.hh"
#include "properties/Integer.hh"
#include "properties/KroneckerDelta.hh"
using namespace cadabra;
reduce_delta::reduce_delta(const Kernel& k, Ex& e)
: Algorithm(k, e)
{
}
bool reduce_delta::can_apply(iterator st)
{
const KroneckerDelta *kr=kernel.properties.get<KroneckerDelta>(st);
if(kr) {
if(tr.number_of_children(st)>2)
return true;
}
return false;
}
Algorithm::result_t reduce_delta::apply(iterator& st)
{
result_t res=result_t::l_no_action;
int num=0;
while(one_step_(st)) {
++num;
res=result_t::l_applied;
if(tr.number_of_children(st)==0) {
st->name=name_set.insert("1").first;
break;
}
};
return res;
}
bool reduce_delta::one_step_(sibling_iterator dl)
{
sibling_iterator up=tr.begin(dl), dn;
int flip=999, masterflip=1;
while(up!=tr.end(dl)) {
flip=masterflip;
dn=tr.begin(dl);
++dn;
while(dn!=tr.end(dl)) {
if(up->name==dn->name) {
if(!up->is_rational()) {
goto found;
}
}
++dn;
++dn;
flip=-flip;
}
++up;
++up;
masterflip=-masterflip;
}
return false;
found:
// std::cerr << "reduce_delta: eliminating " << *up->name << " contraction." << std::endl;
// {unsigned int num=1;
// tr.print_recursive_treeform(debugout, dl, num) << std::endl;}
// FIXME: use properties for the dimension!
// txtout << *dl->multiplier << std::endl;
const Integer *itg=kernel.properties.get<Integer>(up, true);
int dim;
if(itg) {
const nset_t::iterator onept=name_set.insert("1").first;
if(itg->difference.begin()->name==onept)
dim=to_long(*itg->difference.begin()->multiplier);
else
throw ConsistencyException("Summation range for index is not an integer.");
}
else throw ConsistencyException("No dimension known for summation index "+*up->name+".");
int mult=flip*(dim-tr.number_of_children(dl)/2+1);
multiply(dl->multiplier, (multiplier_t)(mult));
multiply(dl->multiplier, multiplier_t(2)/((multiplier_t)(tr.number_of_children(dl))));
// txtout << "flip =" << flip << std::endl;
// remove the indices
sibling_iterator up2=up;
++up2;
++up2;
while(up2!=tr.end(dl)) {
up->name=up2->name;
up->multiplier=up2->multiplier;
++up;
++up;
++up2;
++up2;
}
sibling_iterator dn2=dn;
++dn2;
++dn2;
while(dn2!=tr.end(dl)) {
dn->name=dn2->name;
dn->multiplier=dn2->multiplier;
++dn;
++dn;
++dn2;
++dn2;
}
sibling_iterator lst=tr.end(dl);
--lst;
--lst;
lst=tr.erase(lst);
tr.erase(lst);
// {unsigned int num=1;
// tr.print_recursive_treeform(debugout, dl, num) << std::endl;}
return true;
}
|
