Utilities/interface/Factorize.h

0001 #ifndef PhysicsTools_Utilities_Factorize_h
0002 #define PhysicsTools_Utilities_Factorize_h
0003
0004 #include "PhysicsTools/Utilities/interface/Numerical.h"
0005 #include "PhysicsTools/Utilities/interface/ParametricTrait.h"
0006 #include "PhysicsTools/Utilities/interface/DecomposePower.h"
0007 #include "PhysicsTools/Utilities/interface/DecomposeProduct.h"
0008 #include "PhysicsTools/Utilities/interface/Simplify_begin.h"
0009
0010 #include <type_traits>
0011
0012 namespace funct {
0013   // find a common divider
0014
0015   TEMPL(T1) struct Divides0 {
0016     static const bool value = false;
0017     typedef A arg;
0018     typedef NUM(1) type;
0019     GET(arg, num<1>());
0020   };
0021
0022   TEMPL(T2) struct Divides : public Divides0<A> {};
0023
0024   TEMPL(T2) struct Divides<PROD_S(A, B), void> : public Divides0<PROD_S(A, B)> {};
0025
0026   template <TYPT1, bool par = Parametric<A>::value>
0027   struct ParametricDiv1 : public Divides<A, void> {};
0028
0029   TEMPL(T1) struct ParametricDiv1<A, false> {
0030     static const bool value = true;
0031     typedef A arg;
0032     typedef arg type;
0033     GET(arg, _);
0034   };
0035
0036   TEMPL(T1) struct Divides<A, A> : public ParametricDiv1<A> {};
0037
0038   TEMPL(T2) struct Divides<PROD_S(A, B), PROD_S(A, B)> : public ParametricDiv1<PROD_S(A, B)> {};
0039
0040   TEMPL(N1T1) struct Divides<POWER_S(A, NUM(n)), POWER_S(A, NUM(n))> : public ParametricDiv1<POWER_S(A, NUM(n))> {};
0041
0042   template <TYPN2T1, bool par = Parametric<A>::value>
0043   struct ParametricDivN : public Divides<POWER(A, NUM(n)), void> {};
0044
0045   TEMPL(N2T1) struct ParametricDivN<n, m, A, false> {
0046     static const bool value = true;
0047     typedef POWER(A, NUM(n)) arg;
0048     static const int p = []() {
0049       if constexpr (n < m)
0050         return n;
0051       else
0052         return m;
0053     }();
0054     typedef POWER(A, NUM(p)) type;
0055     typedef DecomposePower<A, NUM(n)> Dec;
0056     GET(arg, pow(Dec::getBase(_), num<p>()));
0057   };
0058
0059   TEMPL(N2T1) struct Divides<POWER_S(A, NUM(n)), POWER_S(A, NUM(m))> : public ParametricDivN<n, m, A> {};
0060
0061   TEMPL(N1T1) struct Divides<A, POWER_S(A, NUM(n))> : public ParametricDivN<1, n, A> {};
0062
0063   TEMPL(N1T1) struct Divides<POWER_S(A, NUM(n)), A> : public ParametricDivN<n, 1, A> {};
0064
0065   namespace tmpl {
0066
0067     template <int n, bool positive = (n >= 0)>
0068     struct abs {
0069       enum { value = n };
0070     };
0071
0072     template <int n>
0073     struct abs<n, false> {
0074       enum { value = -n };
0075     };
0076
0077   }  // namespace tmpl
0078
0079   TEMPL(N2) struct Divides<NUM(n), NUM(m)> {
0080     enum { gcd = std::gcd(tmpl::abs<n>::value, tmpl::abs<m>::value) };
0081     static const bool value = (gcd != 1);
0082     typedef NUM(n) arg;
0083     typedef NUM(gcd) type;
0084     GET(arg, num<gcd>());
0085   };
0086
0087   TEMPL(N1) struct Divides<NUM(n), NUM(n)> {
0088     static const bool value = true;
0089     typedef NUM(n) arg;
0090     typedef arg type;
0091     GET(arg, _);
0092   };
0093
0094   TEMPL(T2) struct Divides<A, MINUS_S(B)> : public Divides<A, B> {};
0095
0096   TEMPL(T3) struct Divides<PROD_S(A, B), MINUS_S(C)> : public Divides<PROD_S(A, B), C> {};
0097
0098   TEMPL(T2) struct Divides<MINUS_S(A), B> {
0099     typedef Divides<A, B> Div;
0100     static const bool value = Div::value;
0101     typedef MINUS_S(A) arg;
0102     typedef typename Div::type type;
0103     GET(arg, Div::get(_._));
0104   };
0105
0106   TEMPL(T2) struct Divides<MINUS_S(A), MINUS_S(B)> {
0107     typedef Divides<A, B> Div;
0108     static const bool value = Div::value;
0109     typedef MINUS_S(A) arg;
0110     typedef typename Div::type type;
0111     GET(arg, Div::get(_._));
0112   };
0113
0114   TEMPL(T3) struct Divides<MINUS_S(A), PROD_S(B, C)> {
0115     typedef Divides<A, PROD_S(B, C)> Div;
0116     static const bool value = Div::value;
0117     typedef MINUS_S(A) arg;
0118     typedef typename Div::type type;
0119     GET(arg, Div::get(_._));
0120   };
0121
0122   TEMPL(T3) struct Divides<A, PROD_S(B, C)> {
0123     typedef A arg;
0124     typedef Divides<arg, void> D0;
0125     typedef Divides<arg, B> D1;
0126     typedef Divides<arg, C> D2;
0127     typedef typename std::conditional<D1::value, D1, typename std::conditional<D2::value, D2, D0>::type>::type Div;
0128     static const bool value = Div::value;
0129     typedef typename Div::type type;
0130     GET(arg, Div::get(_));
0131   };
0132
0133   TEMPL(T3) struct Divides<PROD_S(A, B), C> {
0134     typedef PROD_S(A, B) arg;
0135     typedef Divides<arg, void> D0;
0136     typedef Divides<A, C> D1;
0137     typedef Divides<B, C> D2;
0138     typedef typename std::conditional<D1::value, D1, typename std::conditional<D2::value, D2, D0>::type>::type Div;
0139     typedef typename Div::type type;
0140     static const bool value = Div::value;
0141     typedef DecomposeProduct<arg, typename Div::arg> D;
0142     GET(arg, Div::get(D::get(_)));
0143   };
0144
0145   TEMPL(T4) struct Divides<PROD_S(A, B), PROD_S(C, D)> {
0146     typedef PROD_S(A, B) arg;
0147     typedef Divides<arg, void> D0;
0148     typedef Divides<arg, C> D1;
0149     typedef Divides<arg, D> D2;
0150     typedef typename std::conditional<D1::value, D1, typename std::conditional<D2::value, D2, D0>::type>::type Div;
0151     static const bool value = Div::value;
0152     typedef typename Div::type type;
0153     GET(arg, Div::get(_));
0154   };
0155
0156   /*
0157     TEMPL(T4) struct Divides<RATIO_S(A, B), RATIO_S(C, D)> {
0158     typedef RATIO_S(A, B) arg;
0159     typedef Divides<B, D> Div;
0160     static const bool value = Div::value;
0161     DEF_TYPE(RATIO(NUM(1), typename Div::type))
0162     GET(arg, num(1) / Div::get(_))
0163     };
0164   */
0165
0166   // general factorization algorithm
0167   template <TYPT2, bool factor = Divides<A, B>::value>
0168   struct FactorizeSum {
0169     typedef SUM_S(A, B) type;
0170     COMBINE(A, B, type(_1, _2));
0171   };
0172
0173   TEMPL(T2) struct FactorizeSum<A, B, true> {
0174     typedef typename Divides<A, B>::type F;
0175     typedef PROD(F, SUM(RATIO(A, F), RATIO(B, F))) type;
0176     inline static type combine(const A& _1, const B& _2) {
0177       const F& f = Divides<A, B>::get(_1);
0178       return f * ((_1 / f) + (_2 / f));
0179     }
0180   };
0181
0182   TEMPL(T3) struct Sum<PROD_S(A, B), C> : public FactorizeSum<PROD_S(A, B), C> {};
0183
0184   TEMPL(T3) struct Sum<A, PROD_S(B, C)> : public FactorizeSum<A, PROD_S(B, C)> {};
0185
0186   TEMPL(T3) struct Sum<MINUS_S(PROD_S(A, B)), C> : public FactorizeSum<MINUS_S(PROD_S(A, B)), C> {};
0187
0188   TEMPL(T3) struct Sum<A, MINUS_S(PROD_S(B, C))> : public FactorizeSum<A, MINUS_S(PROD_S(B, C))> {};
0189
0190   TEMPL(T4) struct Sum<PROD_S(A, B), PROD_S(C, D)> : public FactorizeSum<PROD_S(A, B), PROD_S(C, D)> {};
0191
0192   TEMPL(T4)
0193   struct Sum<MINUS_S(PROD_S(A, B)), MINUS_S(PROD_S(C, D))>
0194       : public FactorizeSum<MINUS_S(PROD_S(A, B)), MINUS_S(PROD_S(C, D))> {};
0195
0196   TEMPL(T4)
0197   struct Sum<PROD_S(A, B), MINUS_S(PROD_S(C, D))> : public FactorizeSum<PROD_S(A, B), MINUS_S(PROD_S(C, D))> {};
0198
0199   TEMPL(T4)
0200   struct Sum<MINUS_S(PROD_S(A, B)), PROD_S(C, D)> : public FactorizeSum<MINUS_S(PROD_S(A, B)), PROD_S(C, D)> {};
0201
0202   TEMPL(N1T2) struct Sum<PROD_S(A, B), NUM(n)> : public FactorizeSum<PROD_S(A, B), NUM(n)> {};
0203
0204   TEMPL(N1T2) struct Sum<NUM(n), PROD_S(A, B)> : public FactorizeSum<NUM(n), PROD_S(A, B)> {};
0205
0206   /*
0207     TEMPL(T4) struct Sum<SUM_S(A, B), PROD_S(C, D)> :
0208     public FactorizeSum<SUM_S(A, B), PROD_S(C, D)> { };
0209
0210     TEMPL(T4) struct Sum<SUM_S(A, B), MINUS_S(PROD_S(C, D))> :
0211     public FactorizeSum<SUM_S(A, B), MINUS_S(PROD_S(C, D))> { };
0212
0213     TEMPL(T4) struct Sum<PROD_S(A, B), SUM_S(C, D)> :
0214     public FactorizeSum<PROD_S(A, B), SUM_S(C, D)> { };
0215   */
0216
0217   /*
0218     template <TYPT4, bool factor = Divides<B, D>::value>
0219     struct FactorizeSumRatio {
0220     DEF_TYPE(SUM_S(RATIO_S(A, B), RATIO_S(C, D)))
0221     COMBINE(A, B, type(_1, _2))
0222     };
0223
0224     TEMPL(T4) struct FactorizeSumRatio<A, B, C, D, true> {
0225     typedef typename Divides<B, D>::type F;
0226     DEF_TYPE(PROD(RATIO(NUM(1), F),
0227     SUM(RATIO(PROD(A, F), B),
0228     RATIO(PROD(C, F), D))))
0229     inline static type combine(const RATIO_S(A, B)& _1,
0230     const RATIO_S(C, D)& _2) {
0231     const F& f = Divides<B, D>::get(_1);
0232     return (num(1) / f) * ((_1._1 * f) / _1._2 + (_2._1 * f) / _2._2);
0233     }
0234     };
0235
0236     TEMPL(T4) struct Sum<RATIO_S(A, B), RATIO_S(C, D)> :
0237     public FactorizeSum<RATIO_S(A, B), RATIO_S(C, D)> { };
0238   */
0239 }  // namespace funct
0240
0241 #include "PhysicsTools/Utilities/interface/Simplify_end.h"
0242
0243 #endif