CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10/src/PhysicsTools/Utilities/interface/Factorize.h

Go to the documentation of this file.
00001 #ifndef PhysicsTools_Utilities_Factorize_h
00002 #define PhysicsTools_Utilities_Factorize_h
00003 
00004 #include "PhysicsTools/Utilities/interface/Numerical.h"
00005 #include "PhysicsTools/Utilities/interface/ParametricTrait.h"
00006 #include "PhysicsTools/Utilities/interface/DecomposePower.h"
00007 #include "PhysicsTools/Utilities/interface/DecomposeProduct.h"
00008 #include <boost/type_traits.hpp>
00009 #include <boost/mpl/and.hpp>
00010 #include <boost/mpl/not.hpp>
00011 #include <boost/mpl/int.hpp>
00012  
00013 #include "PhysicsTools/Utilities/interface/Simplify_begin.h"
00014 
00015 namespace funct {
00016   // find a common divider
00017 
00018   TEMPL(T1) struct Divides0 {
00019     static const bool value = false;
00020     typedef A arg;
00021     typedef NUM(1) type;
00022     GET(arg, num<1>());
00023   };
00024 
00025   TEMPL(T2) struct Divides : 
00026     public Divides0<A> { };
00027 
00028   TEMPL(T2) struct Divides<PROD_S(A, B), void> : 
00029     public Divides0<PROD_S(A, B)> { };
00030   
00031   template <TYPT1, bool par = Parametric<A>::value> struct ParametricDiv1 :
00032     public Divides<A, void> { };
00033 
00034   TEMPL(T1) struct ParametricDiv1<A, false> {
00035     static const bool value = true;
00036     typedef A arg;
00037     typedef arg type;
00038     GET(arg,  _);
00039   };
00040 
00041   TEMPL(T1) struct Divides<A, A> : 
00042     public ParametricDiv1<A> { };
00043   
00044   TEMPL(T2) struct Divides<PROD_S(A, B), PROD_S(A, B)> : 
00045     public ParametricDiv1<PROD_S(A, B)> { };
00046   
00047   TEMPL(N1T1) struct Divides<POWER_S(A, NUM(n)), 
00048     POWER_S(A, NUM(n))> : 
00049     public ParametricDiv1<POWER_S(A, NUM(n))> { };
00050   
00051   template <TYPN2T1, bool par = Parametric<A>::value> struct ParametricDivN :
00052     public Divides<POWER(A, NUM(n)), void> { };
00053   
00054   TEMPL(N2T1) struct ParametricDivN<n, m, A, false> {
00055     static const bool value = true;
00056     typedef POWER(A, NUM(n)) arg;
00057     static const int p = ::boost::mpl::if_c<(n <m), 
00058       ::boost::mpl::int_<n>, ::boost::mpl::int_<m> >::type::value;
00059     typedef POWER(A, NUM(p)) type;
00060     typedef DecomposePower<A, NUM(n)> Dec;
00061     GET(arg, pow(Dec::getBase(_), num<p>()));
00062   };
00063 
00064   TEMPL(N2T1) struct Divides<POWER_S(A, NUM(n)), 
00065     POWER_S(A, NUM(m))> : 
00066     public ParametricDivN<n, m, A> { };
00067   
00068   TEMPL(N1T1) struct Divides<A, POWER_S(A, NUM(n))> : 
00069     public ParametricDivN<1, n, A> { };
00070   
00071   TEMPL(N1T1) struct Divides<POWER_S(A, NUM(n)), A> : 
00072     public ParametricDivN<n, 1, A> { };
00073   
00074   namespace tmpl {
00075     
00076     template<int n, bool positive = (n>= 0)>
00077       struct abs { enum { value = n }; };
00078     
00079     template<int n>
00080       struct abs<n, false> { enum { value = -n }; };
00081     
00082   }
00083   
00084   TEMPL(N2) struct Divides<NUM(n), NUM(m)> { 
00085     enum { gcd = boost::math::static_gcd<tmpl::abs<n>::value, tmpl::abs<m>::value>::value };
00086     static const bool value = (gcd != 1);
00087     typedef NUM(n) arg;
00088     typedef NUM(gcd) type;
00089     GET(arg, num<gcd>());
00090   };
00091 
00092   TEMPL(N1) struct Divides<NUM(n), NUM(n)> { 
00093     static const bool value = true;
00094     typedef NUM(n) arg;
00095     typedef arg type;
00096     GET(arg, _);
00097   };
00098 
00099   TEMPL(T2) struct Divides<A, MINUS_S(B)> : 
00100     public Divides<A, B> { };
00101   
00102   TEMPL(T3) struct Divides<PROD_S(A, B), MINUS_S(C)> : 
00103     public Divides<PROD_S(A, B), C> { };
00104   
00105   TEMPL(T2) struct Divides<MINUS_S(A), B> { 
00106     typedef Divides<A, B> Div;
00107     static const bool value = Div::value;
00108     typedef MINUS_S(A) arg;
00109     typedef typename Div::type type;
00110     GET(arg, Div::get(_._));
00111   };
00112 
00113   TEMPL(T2) struct Divides<MINUS_S(A), MINUS_S(B)> {
00114     typedef Divides<A, B> Div;
00115     static const bool value = Div::value;
00116     typedef MINUS_S(A) arg;
00117     typedef typename Div::type type;
00118     GET(arg, Div::get(_._));
00119   };
00120 
00121   TEMPL(T3) struct Divides<MINUS_S(A), PROD_S(B, C)> {
00122     typedef Divides<A, PROD_S(B, C)> Div;
00123     static const bool value = Div::value;
00124     typedef MINUS_S(A) arg;
00125     typedef typename Div::type type;
00126     GET(arg, Div::get(_._));
00127   };
00128 
00129   TEMPL(T3) struct Divides<A, PROD_S(B, C)> {
00130     typedef A arg;
00131     typedef Divides<arg, void> D0;
00132     typedef Divides<arg, B> D1;
00133     typedef Divides<arg, C> D2;
00134     typedef typename ::boost::mpl::if_<D1, D1,
00135       typename ::boost::mpl::if_ <D2, D2, D0>::type>::type Div;
00136     static const bool value = Div::value;
00137     typedef typename Div::type type;
00138     GET(arg, Div::get(_));
00139   };
00140 
00141   TEMPL(T3) struct Divides<PROD_S(A, B), C> {
00142     typedef PROD_S(A, B) arg;
00143     typedef Divides<arg, void> D0;
00144     typedef Divides<A, C> D1;
00145     typedef Divides<B, C> D2;
00146     typedef typename ::boost::mpl::if_<D1, D1,
00147       typename ::boost::mpl::if_ <D2, D2, D0>::type>::type Div;
00148     typedef typename Div::type type;
00149     static const bool value = Div::value;
00150     typedef DecomposeProduct<arg, typename Div::arg> D;
00151     GET(arg, Div::get(D::get(_)));
00152   };
00153   
00154   TEMPL(T4) struct Divides<PROD_S(A, B), PROD_S(C, D)> {
00155     typedef PROD_S(A, B) arg;
00156     typedef Divides<arg, void> D0;
00157     typedef Divides<arg, C> D1;
00158     typedef Divides<arg, D> D2;
00159     typedef typename ::boost::mpl::if_<D1, D1,
00160       typename ::boost::mpl::if_ <D2, D2, D0>::type>::type Div;
00161     static const bool value = Div::value;
00162     typedef typename Div::type type;
00163     GET(arg, Div::get(_));
00164   };
00165 
00166   /*
00167     TEMPL(T4) struct Divides<RATIO_S(A, B), RATIO_S(C, D)> {
00168     typedef RATIO_S(A, B) arg;
00169     typedef Divides<B, D> Div;
00170     static const bool value = Div::value;
00171     DEF_TYPE(RATIO(NUM(1), typename Div::type))
00172     GET(arg, num(1) / Div::get(_))
00173     };
00174   */
00175   
00176   // general factorization algorithm
00177   template <TYPT2, bool factor = Divides<A, B>::value>
00178     struct FactorizeSum {
00179       typedef SUM_S(A, B) type;
00180       COMBINE(A, B, type(_1, _2));
00181     };
00182 
00183   TEMPL(T2) struct FactorizeSum<A, B, true> {
00184     typedef typename Divides<A, B>::type F;
00185     typedef PROD(F, SUM(RATIO(A, F), RATIO(B, F))) type;
00186     inline static type combine(const A& _1, const B& _2) { 
00187       const F& f = Divides<A, B>::get(_1);       
00188       return f * ((_1 / f) + (_2 / f));
00189     }
00190   };
00191 
00192   TEMPL(T3) struct Sum<PROD_S(A, B), C> : 
00193     public FactorizeSum<PROD_S(A, B), C> { };
00194   
00195   TEMPL(T3) struct Sum<A, PROD_S(B, C)> : 
00196     public FactorizeSum<A, PROD_S(B, C)> { };
00197   
00198   TEMPL(T3) struct Sum<MINUS_S(PROD_S(A, B)), C> : 
00199     public FactorizeSum<MINUS_S(PROD_S(A, B)), C> { };
00200   
00201   TEMPL(T3) struct Sum<A, MINUS_S(PROD_S(B, C))> : 
00202     public FactorizeSum<A, MINUS_S(PROD_S(B, C))> { };
00203   
00204   TEMPL(T4) struct Sum<PROD_S(A, B), PROD_S(C, D)> : 
00205     public FactorizeSum<PROD_S(A, B), PROD_S(C, D)> { };
00206   
00207   TEMPL(T4) struct Sum<MINUS_S(PROD_S(A, B)), 
00208     MINUS_S(PROD_S(C, D))> : 
00209     public FactorizeSum<MINUS_S(PROD_S(A, B)), 
00210     MINUS_S(PROD_S(C, D))> { };
00211   
00212   TEMPL(T4) struct Sum<PROD_S(A, B), 
00213     MINUS_S(PROD_S(C, D))> : 
00214     public FactorizeSum<PROD_S(A, B), 
00215     MINUS_S(PROD_S(C, D))> { };
00216   
00217   TEMPL(T4) struct Sum<MINUS_S(PROD_S(A, B)), 
00218     PROD_S(C, D)> : 
00219     public FactorizeSum<MINUS_S(PROD_S(A, B)), 
00220     PROD_S(C, D)> { };
00221 
00222   TEMPL(N1T2) struct Sum<PROD_S(A, B), NUM(n)> : 
00223     public FactorizeSum<PROD_S(A, B), NUM(n)> { };
00224   
00225   TEMPL(N1T2) struct Sum<NUM(n), PROD_S(A, B)> : 
00226     public FactorizeSum<NUM(n), PROD_S(A, B)> { };
00227   
00228   /*
00229     TEMPL(T4) struct Sum<SUM_S(A, B), PROD_S(C, D)> : 
00230     public FactorizeSum<SUM_S(A, B), PROD_S(C, D)> { };
00231     
00232     TEMPL(T4) struct Sum<SUM_S(A, B), MINUS_S(PROD_S(C, D))> : 
00233     public FactorizeSum<SUM_S(A, B), MINUS_S(PROD_S(C, D))> { };
00234     
00235     TEMPL(T4) struct Sum<PROD_S(A, B), SUM_S(C, D)> : 
00236     public FactorizeSum<PROD_S(A, B), SUM_S(C, D)> { };
00237   */
00238   
00239   /*
00240     template <TYPT4, bool factor = Divides<B, D>::value>
00241     struct FactorizeSumRatio {
00242     DEF_TYPE(SUM_S(RATIO_S(A, B), RATIO_S(C, D)))
00243     COMBINE(A, B, type(_1, _2))
00244     };
00245     
00246     TEMPL(T4) struct FactorizeSumRatio<A, B, C, D, true> {
00247     typedef typename Divides<B, D>::type F;
00248     DEF_TYPE(PROD(RATIO(NUM(1), F), 
00249     SUM(RATIO(PROD(A, F), B), 
00250     RATIO(PROD(C, F), D))))
00251     inline static type combine(const RATIO_S(A, B)& _1, 
00252     const RATIO_S(C, D)& _2) { 
00253     const F& f = Divides<B, D>::get(_1);       
00254     return (num(1) / f) * ((_1._1 * f) / _1._2 + (_2._1 * f) / _2._2);
00255     }
00256     };
00257     
00258     TEMPL(T4) struct Sum<RATIO_S(A, B), RATIO_S(C, D)> : 
00259     public FactorizeSum<RATIO_S(A, B), RATIO_S(C, D)> { };
00260   */
00261 }
00262 
00263 #include "PhysicsTools/Utilities/interface/Simplify_end.h"
00264 
00265 #endif