METOOLS/Loops/Divergence_Array.H

0001 #ifndef METOOLS_Loops_Divergence_Array_H
0002 #define METOOLS_Loops_Divergence_Array_H
0003
0004 #include "ATOOLS/Math/MyComplex.H"
0005 #include "ATOOLS/Math/MathTools.H"
0006 #include "ATOOLS/Org/Exception.H"
0007 #include <vector>
0008
0009 using namespace ATOOLS;
0010
0011 namespace METOOLS {
0012
0013   template<typename T>
0014   class Divergence_Array {
0015   private:
0016     std::vector<T> m_result;
0017
0018   public:
0019     Divergence_Array() : m_result(6) {
0020     }
0021     Divergence_Array(const T& t) : m_result(6, t) {
0022     }
0023     Divergence_Array(const T& uv, const T& ir, const T& ir2,
0024                      const T& fin, const T& eps, const T& eps2) {
0025       m_result.reserve(6);
0026       m_result.push_back(uv);
0027       m_result.push_back(ir);
0028       m_result.push_back(ir2);
0029       m_result.push_back(fin);
0030       m_result.push_back(eps);
0031       m_result.push_back(eps2);
0032     }
0033     Divergence_Array(const std::vector<T>& res) : m_result(res) {
0034       if (m_result.size()!=6) THROW(fatal_error, "Internal error.");
0035     }
0036     ~Divergence_Array() {}
0037
0038     // extraction of entries
0039     inline const T&   GetUV()       const { return m_result[0]; }
0040     inline const T&   GetIR()       const { return m_result[1]; }
0041     inline const T&   GetIR2()      const { return m_result[2]; }
0042     inline const T&   GetFinite()   const { return m_result[3]; }
0043     inline const T&   GetEpsilon()  const { return m_result[4]; }
0044     inline const T&   GetEpsilon2() const { return m_result[5]; }
0045     // extraction of entries
0046     inline T&   UV()       { return m_result[0]; }
0047     inline T&   IR()       { return m_result[1]; }
0048     inline T&   IR2()      { return m_result[2]; }
0049     inline T&   Finite()   { return m_result[3]; }
0050     inline T&   Epsilon()  { return m_result[4]; }
0051     inline T&   Epsilon2() { return m_result[5]; }
0052
0053     inline const std::vector<T>& GetResult() const { return m_result; }
0054
0055     inline T&       operator[] (unsigned short int i)        {
0056       return m_result[i];
0057     }
0058
0059     inline const T& operator[] (unsigned short int i) const  {
0060       return m_result[i];
0061     }
0062
0063     // sign flip
0064     inline Divergence_Array<T> operator- () const {
0065       return Divergence_Array<T>(-m_result[0],-m_result[1],-m_result[2],
0066                                  -m_result[3],-m_result[4],-m_result[5]);
0067     }
0068
0069     // addition/subtraction
0070     template<typename T1>
0071     inline Divergence_Array<PROMOTE(T,T1)>
0072     operator+ (const Divergence_Array<T1>& da) const {
0073       return Divergence_Array<PROMOTE(T,T1)>(m_result[0]+da[0],
0074                                              m_result[1]+da[1],
0075                                              m_result[2]+da[2],
0076                                              m_result[3]+da[3],
0077                                              m_result[4]+da[4],
0078                                              m_result[5]+da[5]);
0079     }
0080
0081     inline Divergence_Array<T>&
0082     operator+=(const Divergence_Array<T>& da) {
0083       m_result[0]+=da[0];
0084       m_result[1]+=da[1];
0085       m_result[2]+=da[2];
0086       m_result[3]+=da[3];
0087       m_result[4]+=da[4];
0088       m_result[5]+=da[5];
0089       return *this;
0090     }
0091
0092
0093     template<typename T1>
0094     inline Divergence_Array<PROMOTE(T,T1)>
0095     operator- (const Divergence_Array<T1>& da) const {
0096       return Divergence_Array<PROMOTE(T,T1)>(m_result[0]-da[0],
0097                                              m_result[1]-da[1],
0098                                              m_result[2]-da[2],
0099                                              m_result[3]-da[3],
0100                                              m_result[4]-da[4],
0101                                              m_result[5]-da[5]);
0102     }
0103
0104     template<typename T1>
0105     inline Divergence_Array<PROMOTE(T,T1)>
0106     operator+ (const T1& c) const {
0107       return Divergence_Array<PROMOTE(T,T1)>(m_result[0],
0108                                              m_result[1],
0109                                              m_result[2],
0110                                              m_result[3]+c,
0111                                              m_result[4],
0112                                              m_result[5]);
0113     }
0114
0115     inline Divergence_Array<T>&
0116     operator+=(const T& c) {
0117       m_result[3]+=c;
0118       return *this;
0119     }
0120
0121     template<typename T1>
0122     inline Divergence_Array<PROMOTE(T,T1)>
0123     operator- (const T1& c) const {
0124       return Divergence_Array<PROMOTE(T,T1)>(m_result[0],
0125                                              m_result[1],
0126                                              m_result[2],
0127                                              m_result[3]-c,
0128                                              m_result[4],
0129                                              m_result[5]);
0130     }
0131
0132     // multiplication
0133     // keep in mind:
0134     //   - terms of order 1/epsUV^2, 1/epsIR^3, 1/epsIR^4, eps^2 are not calculated
0135     //   - also overlapping IR/UV divergences are not calculated yet,
0136     //     -> their treatment needs to be thought about
0137     //   - terms 1/epsIR^2 * epsUV, 1/epsIR^2 * epsIR, etc. are not differentiated properly
0138     template<typename T1>
0139     inline Divergence_Array<PROMOTE(T,T1)>
0140     operator* (const Divergence_Array<T1>& da) const {
0141       return Divergence_Array<PROMOTE(T,T1)>(m_result[0]*da[3]+da[0]*m_result[3],
0142                                              m_result[1]*da[3]+da[1]*m_result[3]
0143                                             +m_result[2]*da[4]+da[2]*m_result[4],
0144                                              m_result[2]*da[3]+da[2]*m_result[3]
0145                                             +m_result[1]*da[1]+da[1]*m_result[1],
0146                                              m_result[3]*da[3]+da[3]*m_result[3]
0147                                             +m_result[0]*da[4]+da[0]*m_result[4]
0148                                             +m_result[1]*da[4]+da[1]*m_result[4]
0149                                             +m_result[2]*da[5]+da[2]*m_result[5],
0150                                              m_result[4]*da[3]+da[4]*m_result[3]
0151                                             +m_result[0]*da[5]+da[0]*m_result[5]
0152                                             +m_result[1]*da[5]+da[1]*m_result[5],
0153                                              m_result[5]*da[3]+da[5]*m_result[3]);
0154     }
0155
0156     // multiplication with scalar
0157     template<typename T1>
0158     inline Divergence_Array<PROMOTE(T,T1)>
0159     operator* (const T1& s) const {
0160       return Divergence_Array<PROMOTE(T,T1)>(s*m_result[0],
0161                                              s*m_result[1],
0162                                              s*m_result[2],
0163                                              s*m_result[3],
0164                                              s*m_result[4],
0165                                              s*m_result[5]);
0166     }
0167
0168     inline Divergence_Array<T>&
0169     operator*=(const T& s) {
0170       m_result[0]*=s;
0171       m_result[1]*=s;
0172       m_result[2]*=s;
0173       m_result[3]*=s;
0174       m_result[4]*=s;
0175       m_result[5]*=s;
0176       return *this;
0177     }
0178
0179     template<typename T1>
0180     inline Divergence_Array<PROMOTE(T,T1)>
0181     operator/ (const T1& s) const {
0182       return Divergence_Array<PROMOTE(T,T1)>(m_result[0]/s,
0183                                              m_result[1]/s,
0184                                              m_result[2]/s,
0185                                              m_result[3]/s,
0186                                              m_result[4]/s,
0187                                              m_result[5]/s);
0188     }
0189
0190   }; // end of class Divergence_Array
0191
0192   // typedefs
0193   typedef Divergence_Array<double>  DivArrD;
0194   typedef Divergence_Array<Complex> DivArrC;
0195
0196   template<typename T1, typename T2>
0197   inline Divergence_Array<PROMOTE(T1,T2)>
0198   operator+ (const T1& s, const Divergence_Array<T2>& a) {
0199     return a+s;
0200   }
0201
0202   template<typename T1, typename T2>
0203   inline Divergence_Array<PROMOTE(T1,T2)>
0204   operator- (const T1& s, const Divergence_Array<T2>& a) {
0205     return -a+s;
0206   }
0207
0208   template<typename T1, typename T2>
0209   inline Divergence_Array<PROMOTE(T1,T2)>
0210   operator* (const T1& s, const Divergence_Array<T2>& a) {
0211     return a*s;
0212   }
0213
0214   template<typename T1, typename T2>
0215   inline Divergence_Array<PROMOTE(T1,T2)>
0216   operator/ (const T1& s, const Divergence_Array<T2>& a) {
0217     // only implemented for DivArrs with zero 1/eps parts
0218     // 1/(a+b*eps+c*eps^2)
0219     // = 1/a - b/a^2 * eps + (2b^2-ca)/(2a^3) * eps^2
0220     return s*Divergence_Array<PROMOTE(T1,T2)>(0.,
0221                                               0.,
0222                                               0.,
0223                                               1./a[3],
0224                                               -a[4]/sqr(a[3]),
0225                                               (2.*sqr(a[4])-a[5]*a[3])
0226                                                 /(2.*a[3]*a[3]*a[3]));
0227   }
0228
0229   template<typename T1, typename T2>
0230   inline Divergence_Array<PROMOTE(T1,T2)>
0231   operator/ (const Divergence_Array<T1>& a1, const Divergence_Array<T2>& a2) {
0232     return a1*(1./a2);
0233   }
0234
0235   inline DivArrC conj(DivArrC divarr) {
0236     return DivArrC(std::conj(divarr[0]),
0237                    std::conj(divarr[1]),
0238                    std::conj(divarr[2]),
0239                    std::conj(divarr[3]),
0240                    std::conj(divarr[4]),
0241                    std::conj(divarr[5]));
0242   }
0243
0244   inline DivArrD real(DivArrC divarr) {
0245     return DivArrD(std::real(divarr[0]),
0246                    std::real(divarr[1]),
0247                    std::real(divarr[2]),
0248                    std::real(divarr[3]),
0249                    std::real(divarr[4]),
0250                    std::real(divarr[5]));
0251   }
0252
0253   inline DivArrD imag(DivArrC divarr) {
0254     return DivArrD(std::imag(divarr[0]),
0255                    std::imag(divarr[1]),
0256                    std::imag(divarr[2]),
0257                    std::imag(divarr[3]),
0258                    std::imag(divarr[4]),
0259                    std::imag(divarr[5]));
0260   }
0261
0262 } // end of namespace METOOLS
0263
0264 #define RED(ARG) om::red<<ARG<<om::reset
0265 #define GREEN(ARG) om::green<<ARG<<om::reset
0266 #define BLUE(ARG) om::blue<<ARG<<om::reset
0267 #define YELLOW(ARG) om::brown<<ARG<<om::reset
0268 #define BLACK(ARG) ARG
0269 #define BOLD(ARG) om::bold<<ARG<<om::reset
0270
0271 namespace ATOOLS {
0272
0273   template<typename T>
0274   std::ostream& operator<<(std::ostream& s, const METOOLS::Divergence_Array<T>& divarr) {
0275     return s<<BOLD('(')<<BLUE(divarr[0])<<BOLD(',')<<RED(divarr[1])<<BOLD(',')
0276                <<YELLOW(divarr[2])<<BOLD(',')<<BOLD(divarr[3])<<BOLD(',')
0277                <<GREEN(divarr[4])<<BOLD(',')<<GREEN(divarr[5])<<BOLD(')');
0278   }
0279
0280   template<>
0281   inline bool IsNan<METOOLS::DivArrD>(const METOOLS::DivArrD& x) {
0282     for (size_t i=0; i<x.GetResult().size(); ++i)
0283       if (ATOOLS::IsNan(x[i])) return true;
0284     return false;
0285   }
0286
0287   template<>
0288   inline bool IsNan<METOOLS::DivArrC>(const METOOLS::DivArrC& x) {
0289     for (size_t i=0; i<x.GetResult().size(); ++i)
0290       if (ATOOLS::IsNan(x[i])) return true;
0291     return false;
0292   }
0293
0294   template<>
0295   inline bool IsZero<METOOLS::DivArrD>(const METOOLS::DivArrD& x) {
0296     for (size_t i=0; i<x.GetResult().size(); ++i)
0297       if (!ATOOLS::IsZero(x[i])) return false;
0298     return true;
0299   }
0300
0301   template<>
0302   inline bool IsZero<METOOLS::DivArrC>(const METOOLS::DivArrC& x) {
0303     for (size_t i=0; i<x.GetResult().size(); ++i)
0304       if (!ATOOLS::IsZero(x[i])) return false;
0305     return true;
0306   }
0307 }
0308
0309
0310 #endif