HepLib
Pair.cpp
Go to the documentation of this file.
1 
6 #include "HEP.h"
7 
8 namespace HepLib {
9 
10  //-----------------------------------------------------------
11  // Pair Class
12  //-----------------------------------------------------------
13  //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Pair, basic,print_func<print_dflt>(&Pair::print).print_func<FormFormat>(&Pair::form_print).print_func<FCFormat>(&Pair::fc_print))
14  //GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eps, basic,print_func<print_dflt>(&Eps::print).print_func<FormFormat>(&Eps::form_print).print_func<FCFormat>(&Eps::fc_print))
15  GiNaC::registered_class_info & Pair::get_class_info_static() { return reg_info; }
17  Pair * Pair::duplicate() const { Pair * bp = new Pair(*this); bp->setflag(GiNaC::status_flags::dynallocated); return bp; }
18  void Pair::accept(GiNaC::visitor & v) const { if (visitor *p = dynamic_cast<visitor *>(&v)) p->visit(*this); else inherited::accept(v); }
19  const GiNaC::registered_class_info &Pair::get_class_info() const { return get_class_info_static(); }
20  GiNaC::registered_class_info &Pair::get_class_info() { return get_class_info_static(); }
21  const char *Pair::class_name() const { return get_class_info_static().options.get_name(); }
22  //GINAC_IMPLEMENT_REGISTERED_CLASS END
23 
27 
28 
34  Pair::Pair(const Vector &p1, const Vector &p2) {
35  if(p1.compare(p2)>0) {
36  lr[0]=p1;lr[1]=p2;
37  } else {
38  lr[0]=p2;lr[1]=p1;
39  }
40  }
41 
48  Pair::Pair(const Index &i1, const Index &i2) {
49  if(i1.compare(i2)>0) {
50  lr[0]=i1;lr[1]=i2;
51  } else {
52  lr[0]=i2;lr[1]=i1;
53  }
54  }
55 
62  Pair::Pair(const Vector &p, const Index &i) : lr{p, i} { }
63 
70  Pair::Pair(const Index &i, const Vector &p) : lr{p, i} { }
71 
72  int Pair::compare_same_type(const basic &other) const {
73  if(!is_a<Pair>(other)) throw Error("Pair::compare_same_type");
74  const Pair &o = static_cast<const Pair &>(other);
75  int c1 = lr[0].compare(o.lr[0]);
76  if(c1!=0) return c1;
77  int c2 = lr[1].compare(o.lr[1]);
78  return c2;
79  }
80 
81  bool Pair::is_equal_same_type(const basic & other) const {
82  if(!is_a<Pair>(other)) throw Error("Pair::is_equal_same_type");
83  const Pair &o = static_cast<const Pair &>(other);
84  if(!lr[0].is_equal(o.lr[0])) return false;
85  return lr[1].is_equal(o.lr[1]);
86  }
87 
93  void Pair::print(const print_dflt &c, unsigned level) const {
94  c.s << lr[0] << "." << lr[1];
95  }
96 
102  void Pair::form_print(const FormFormat &c, unsigned level) const {
103  if(is_a<Vector>(lr[0]) && is_a<Vector>(lr[1])) c << lr[0] << "." << lr[1];
104  else if(is_a<Vector>(lr[0]) && is_a<Index>(lr[1])) c << lr[0] << "(" << lr[1] << ")";
105  else if(is_a<Index>(lr[0]) && is_a<Index>(lr[1])) c << "d_(" << lr[0] << "," << lr[1] << ")";
106  }
107 
113  void Pair::fc_print(const FCFormat &c, unsigned level) const {
114  if(is_a<Vector>(lr[0]) && is_a<Vector>(lr[1])) c << "SPD[" << lr[0] << "," << lr[1] << "]";
115  else if(is_a<Vector>(lr[0]) && is_a<Index>(lr[1])) c << "FVD[" << lr[0] << "," << lr[1] << "]";
116  else if(is_a<Index>(lr[0]) && is_a<Index>(lr[1])) {
117  auto ii = ex_to<Index>(lr[0]);
118  if(ii.type == Index::Type::VD) c << "MTD[" << lr[0] << "," << lr[1] << "]";
119  else if(ii.type == Index::Type::CF) c << "SUNFDelta[" << lr[0] << "," << lr[1] << "]";
120  else if(ii.type == Index::Type::CA) c << "SUNDelta[" << lr[0] << "," << lr[1] << "]";
121  else throw Error("Pair::fc_print unexpected.");
122  }
123  }
124 
125  size_t Pair::nops() const { return 2; }
126  ex Pair::op(size_t i) const {
127  return lr[i];
128  }
129  ex & Pair::let_op(size_t i) {
130  ensure_if_modifiable();
131  return lr[i];
132  }
133 
134  ex Pair::eval() const {
135  if(flags & status_flags::evaluated) return *this;
136  else if(!is_a<Vector>(lr[0]) && !is_a<Index>(lr[0])) return SP(lr[0],lr[1]);
137  else if(!is_a<Vector>(lr[1]) && !is_a<Index>(lr[1])) return SP(lr[0],lr[1]);
138  else if((is_a<Vector>(lr[0]) && is_a<Vector>(lr[1])))
139  return Pair(ex_to<Vector>(lr[0]),ex_to<Vector>(lr[1])).hold();
140  else if((is_a<Index>(lr[0]) && is_a<Index>(lr[1])))
141  return Pair(ex_to<Index>(lr[0]),ex_to<Index>(lr[1])).hold();
142  else if((is_a<Index>(lr[0]) && is_a<Vector>(lr[1])))
143  return Pair(ex_to<Vector>(lr[1]),ex_to<Index>(lr[0])).hold();
144  else return this->hold();
145  }
146 
147  void Pair::archive(archive_node & n) const {
148  inherited::archive(n);
149  for(int i=0; i<2; i++) n.add_ex("lr"+to_string(i), lr[i]);
150  }
151 
152  void Pair::read_archive(const archive_node& n) {
153  inherited::read_archive(n);
154  for(int i=0; i<2; i++) {
155  n.find_ex("lr"+to_string(i), lr[i]);
156  }
157  }
158 
159  ex Pair::derivative(const symbol & s) const {
160  return 0;
161  }
162 
163  //-----------------------------------------------------------
164  // SP function - ScalarProduct
165  //-----------------------------------------------------------
166  ex SP(const ex & a, bool use_map) { return SP(a,a,use_map); }
167 
175  ex SP(const ex & a, const ex & b, bool use_map) {
176  ex ret_sp;
177  if(is_a<Vector>(a) && is_a<Vector>(b)) ret_sp = Pair(ex_to<Vector>(a), ex_to<Vector>(b));
178  else if(is_a<Vector>(a) && is_a<Index>(b)) ret_sp = Pair(ex_to<Vector>(a), ex_to<Index>(b));
179  else if(is_a<Index>(a) && is_a<Vector>(b)) ret_sp = Pair(ex_to<Vector>(b), ex_to<Index>(a));
180  else if(is_a<Index>(a) && is_a<Index>(b)) ret_sp = Pair(ex_to<Index>(a), ex_to<Index>(b));
181  if(is_a<Pair>(ret_sp)) {
182  if(use_map) return ret_sp.subs(SP_map);
183  else return ret_sp;
184  }
185 
186  lst alst, blst;
187  auto aex = expand(a);
188  auto bex = expand(b);
189  if(is_zero(aex) || is_zero(bex)) return 0;
190  if(is_a<add>(aex)) {
191  for(auto item : aex) alst.append(item);
192  } else alst.append(aex);
193  for(int i=0; i<alst.nops(); i++) {
194  if(!is_a<mul>(alst.op(i))) alst.let_op(i) = lst{alst.op(i)};
195  ex c=1;
196  ex v=1;
197  for(auto ii : alst.op(i)) {
198  if(is_a<Vector>(ii) || is_a<Index>(ii)) {
199  if(!is_zero(v-1)) throw Error("Error Found in SP @1 : "+ex2str(alst));
200  v = ii;
201  } else c *= ii;
202  }
203  if(is_zero(v-1)) throw Error("Error Found in SP @2 : "+ex2str(alst));
204  alst.let_op(i) = lst{c,v};
205  }
206 
207  if(is_a<add>(bex)) {
208  for(auto item : bex) blst.append(item);
209  } else blst.append(bex);
210  for(int i=0; i<blst.nops(); i++) {
211  if(!is_a<mul>(blst.op(i))) blst.let_op(i) = lst{blst.op(i)};
212  ex c=1;
213  ex v=1;
214  for(auto ii : blst.op(i)) {
215  if(is_a<Vector>(ii) || is_a<Index>(ii)) {
216  if(!is_zero(v-1)) throw Error("Error Found in SP @3");
217  v = ii;
218  } else c *= ii;
219  }
220  if(is_zero(v-1)) throw Error("Error Found in SP @4");
221  blst.let_op(i) = lst{c,v};
222  }
223 
224  ex res = 0;
225  for(auto ai : alst) {
226  for(auto bi : blst) res += ai.op(0) * bi.op(0) * SP(ai.op(1), bi.op(1), use_map);
227  }
228  return res;
229  }
230 
237  ex sp(const ex & a, const ex & b) {
238  if(is_a<Vector>(a) && is_a<Vector>(b)) return ex_to<Vector>(a).name * ex_to<Vector>(b).name;
239  else throw Error("Error in sp(2).");
240  }
241  ex sp(const ex & a) {
242  if(is_a<Vector>(a)) return ex_to<Vector>(a).name * ex_to<Vector>(a).name;
243  else throw Error("Error in sp(1).");
244  }
245 
252  ex & letSP(const ex &p1, const ex &p2) {
253  if(!(is_a<Vector>(p1) || is_a<Index>(p1)) || !(is_a<Vector>(p2) || is_a<Index>(p2)))
254  throw Error("Invalide arguments for letSP (2).");
255  return SP_map[SP(p1,p2)];
256  }
257 
263  ex & letSP(const ex &p) {
264  if(!is_a<Vector>(p))
265  throw Error("Invalide arguments for letSP (1).");
266  return SP_map[SP(p)];
267  }
268 
274  void clearSP(const ex &p1, const ex &p2) {
275  if(!(is_a<Vector>(p1) || is_a<Index>(p1)) || !(is_a<Vector>(p2) || is_a<Index>(p2)))
276  throw Error("Invalide arguments for resetSP (2).");
277  SP_map.erase(SP(p1,p2));
278  }
279 
284  void clearSP(const ex &p) {
285  if(!is_a<Vector>(p))
286  throw Error("Invalide arguments for resetSP (1).");
287  SP_map.erase(SP(p));
288  }
289 
293  void clearSP() {
294  SP_map.clear();
295  }
296 
302  ex SP2sp(const ex & exin) {
303  auto ret = exin.subs(SP_map);
304  lst vecs = Vector::all(ret);
305  exmap spmap;
306  for(auto vp1 : vecs) {
307  for(auto vp2 : vecs) {
308  spmap[SP(vp1,vp2)]=sp(vp1,vp2);
309  }
310  }
311  return ret.subs(spmap);
312  }
313 
318  exmap sp_map() {
319  exmap ret;
320  for(auto kv : SP_map) {
321  auto key = kv.first;
322  if(key.nops()!=2 || !is_a<Vector>(key.op(0)) || !is_a<Vector>(key.op(1))) continue;
323  auto p1 = ex_to<Vector>(key.op(0));
324  auto p2 = ex_to<Vector>(key.op(1));
325  ret[p1.name * p2.name] = kv.second.subs(SP_map);
326  }
327  return ret;
328  }
329 
330 }
331 
int * a
Definition: Functions.cpp:234
#define IMPLEMENT_HAS(classname)
Definition: BASIC.h:24
#define DEFAULT_CTOR(classname)
Definition: BASIC.h:21
#define IMPLEMENT_ALL(classname)
Definition: BASIC.h:30
HEP header file.
class used to wrap error message
Definition: BASIC.h:242
class for FCFormat Output
Definition: HEP.h:71
class for FormFormat Output
Definition: HEP.h:44
class for index object
Definition: HEP.h:104
virtual ~visitor()
Definition: Pair.cpp:16
class for Pair object
Definition: HEP.h:322
void read_archive(const archive_node &n) override
Definition: Pair.cpp:152
int compare_same_type(const GiNaC::basic &other) const override
Definition: Pair.cpp:72
const GiNaC::registered_class_info & get_class_info() const override
Definition: Pair.cpp:19
ex & let_op(size_t i) override
Definition: Pair.cpp:129
static GiNaC::registered_class_info & get_class_info_static()
Definition: Pair.cpp:15
ex eval() const override
Definition: Pair.cpp:134
ex op(size_t i) const override
Definition: Pair.cpp:126
void archive(archive_node &n) const override
Definition: Pair.cpp:147
size_t nops() const override
Definition: Pair.cpp:125
void accept(GiNaC::visitor &v) const override
Definition: Pair.cpp:18
ex derivative(const symbol &s) const override
Definition: Pair.cpp:159
bool is_equal_same_type(const basic &other) const override
Definition: Pair.cpp:81
void form_print(const FormFormat &c, unsigned level=0) const
FormFormat print function.
Definition: Pair.cpp:102
void print(const print_dflt &c, unsigned level=0) const
default print function
Definition: Pair.cpp:93
const char * class_name() const override
Definition: Pair.cpp:21
void fc_print(const FCFormat &c, unsigned level=0) const
FCFormat print function.
Definition: Pair.cpp:113
Pair * duplicate() const override
Definition: Pair.cpp:17
class for vector object
Definition: HEP.h:149
static lst all(const ex &e)
HepLib namespace.
Definition: BASIC.cpp:17
ex SP2sp(const ex &exin)
convert SP(a,b) to sp(a,b)
Definition: Pair.cpp:302
ex sp(const ex &a, const ex &b)
translated the vector dot a.b to a*b, useful in SecDec
Definition: Pair.cpp:237
exmap sp_map()
the SP_map with SP(a,b) replaced to sp(a,b)
Definition: Pair.cpp:318
exmap SP_map
Definition: Init.cpp:179
const Symbol CA
ex & letSP(const ex &p1, const ex &p2)
return the reference of p1.p2
Definition: Pair.cpp:252
const Symbol CF
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition: BASIC.cpp:715
void clearSP(const ex &p1, const ex &p2)
delete the assignment of p1.p2
Definition: Pair.cpp:274
ex SP(const ex &a, bool use_map=false)
Definition: Pair.cpp:166