Utilities/interface/Primitive.h

0001 #ifndef PhysicsTools_Utilities_Primitive_h
0002 #define PhysicsTools_Utilities_Primitive_h
0003 #include "PhysicsTools/Utilities/interface/Functions.h"
0004 #include "PhysicsTools/Utilities/interface/Operations.h"
0005 #include "PhysicsTools/Utilities/interface/Fraction.h"
0006 #include "PhysicsTools/Utilities/interface/Derivative.h"
0007 #include "PhysicsTools/Utilities/interface/Parameter.h"
0008 #include "PhysicsTools/Utilities/interface/Identity.h"
0009 #include <type_traits>
0010
0011 #include "PhysicsTools/Utilities/interface/Simplify_begin.h"
0012
0013 namespace funct {
0014
0015   struct no_var;
0016
0017   struct UndefinedIntegral {};
0018
0019   template <typename X, typename F, bool independent = Independent<X, F>::value>
0020   struct ConstPrimitive {
0021     typedef UndefinedIntegral type;
0022     inline static type get(const F& f) { return type(); }
0023   };
0024
0025   //  /
0026   //  | c dx = c * x
0027   //  /
0028   template <typename X, typename F>
0029   struct ConstPrimitive<X, F, true> {
0030     typedef PROD(F, X) type;
0031     inline static type get(const F& f) { return f * X(); }
0032   };
0033
0034   template <typename F, typename X = no_var>
0035   struct Primitive : public ConstPrimitive<X, F> {};
0036
0037   template <typename F>
0038   struct Primitive<F> {};
0039
0040   template <typename X, typename F>
0041   typename Primitive<F, X>::type primitive(const F& f) {
0042     return Primitive<F, X>::get(f);
0043   }
0044
0045   template <typename F>
0046   typename Primitive<F>::type primitive(const F& f) {
0047     return Primitive<F>::get(f);
0048   }
0049
0050   //  /
0051   //  | f ^ g dx : UNDEFINED
0052   //  /
0053   PRIMIT_RULE(TYPXT2, POWER_S(A, B), UndefinedIntegral, type());
0054
0055   //  /
0056   //  | x dx = x ^ 2 / 2
0057   //  /
0058   PRIMIT_RULE(TYPX, X, RATIO(POWER(X, NUM(2)), NUM(2)), pow(_, num<2>()) / num<2>());
0059
0060   //  /
0061   //  | x ^ n dx = x ^ (n + 1) / (n + 1)
0062   //  /
0063   PRIMIT_RULE(TYPXN1,
0064               POWER_S(X, NUM(n)),
0065               RATIO(POWER(X, NUM(n + 1)), NUM(n + 1)),
0066               pow(_._1, num<n + 1>()) / num<n + 1>());
0067
0068   //  /
0069   //  | 1 / x ^ n dx = (- 1) / ((n - 1)  x ^ (n - 1))
0070   //  /
0071   PRIMIT_RULE(TYPXN1,
0072               RATIO_S(NUM(1), POWER_S(X, NUM(n))),
0073               RATIO(NUM(-1), PROD(NUM(n - 1), POWER(X, NUM(n - 1)))),
0074               num<-1>() / (num<n - 1>() * pow(_._2._1, num<n - 1>())));
0075
0076   PRIMIT_RULE(TYPXN1,
0077               POWER_S(RATIO_S(NUM(1), X), NUM(n)),
0078               RATIO(NUM(-1), PROD(NUM(n - 1), POWER(X, NUM(n - 1)))),
0079               num<-1>() / (num<n - 1>() * pow(_._1._2, num<n - 1>())));
0080
0081   //  /
0082   //  | x ^ n/m dx = m / (n + m) (x)^ (n + m / m)
0083   //  /
0084   PRIMIT_RULE(TYPXN2,
0085               POWER_S(X, FRACT_S(n, m)),
0086               PROD(FRACT(m, n + m), POWER(X, FRACT(n + m, m))),
0087               (fract<m, n + m>() * pow(_._1, fract<n + m, m>())));
0088
0089   //  /
0090   //  | sqrt(x) dx = 2/3 (x)^ 3/2
0091   //  /
0092   PRIMIT_RULE(TYPX, SQRT_S(X), PRIMIT(X, POWER_S(X, FRACT_S(1, 2))), (fract<2, 3>() * pow(_._, fract<3, 2>())));
0093
0094   //  /
0095   //  | exp(x) dx = exp(x)
0096   //  /
0097   PRIMIT_RULE(TYPX, EXP_S(X), EXP(X), _);
0098
0099   //  /
0100   //  | log(x) dx = x(log(x) - 1)
0101   //  /
0102   PRIMIT_RULE(TYPX, LOG_S(X), PROD(X, DIFF(LOG(X), NUM(1))), _._*(_ - num<1>()));
0103
0104   //  /
0105   //  | sgn(x) dx = abs(x)
0106   //  /
0107   PRIMIT_RULE(TYPX, SGN_S(X), ABS(X), abs(_._));
0108
0109   //  /
0110   //  | sin(x) dx = - cos(x)
0111   //  /
0112   PRIMIT_RULE(TYPX, SIN_S(X), MINUS(COS(X)), -cos(_._));
0113
0114   //  /
0115   //  | cos(x) dx = sin(x)
0116   //  /
0117   PRIMIT_RULE(TYPX, COS_S(X), SIN(X), sin(_._));
0118
0119   //  /
0120   //  | tan(x) dx = - log(abs(cos(x)))
0121   //  /
0122   PRIMIT_RULE(TYPX, TAN_S(X), MINUS(LOG(ABS(COS(X)))), -log(abs(cos(_._))));
0123
0124   //  /
0125   //  | 1 / x dx = log(abs(x))
0126   //  /
0127   PRIMIT_RULE(TYPX, RATIO_S(NUM(1), X), LOG(ABS(X)), log(abs(_._2)));
0128
0129   PRIMIT_RULE(TYPX, POWER_S(X, NUM(-1)), LOG(ABS(X)), log(abs(_._1)));
0130
0131   //  /
0132   //  | 1 / cos(x)^2 dx = tan(x)
0133   //  /
0134   PRIMIT_RULE(TYPX, RATIO_S(NUM(1), POWER_S(COS_S(X), NUM(2))), TAN(X), tan(_._2._1._));
0135
0136   //  /
0137   //  | 1 / sin(x)^2 dx = - 1 / tan(x)
0138   //  /
0139   PRIMIT_RULE(TYPX, RATIO_S(NUM(1), POWER_S(SIN_S(X), NUM(2))), RATIO(NUM(-1), TAN(X)), num<-1>() / tan(_._2._1._));
0140
0141   // composite primitives
0142
0143   //  /                    /           /
0144   //  | (f(x) + g(x)) dx = | f(x) dx + | g(x) dx
0145   //  /                    /           /
0146   PRIMIT_RULE(TYPXT2, SUM_S(A, B), SUM(PRIMIT(X, A), PRIMIT(X, B)), primitive<X>(_._1) + primitive<X>(_._2));
0147
0148   //  /                 /
0149   //  | (- f(x)) dx = - | f(x) dx
0150   //  /                 /
0151   PRIMIT_RULE(TYPXT1, MINUS_S(A), MINUS(PRIMIT(X, A)), -primitive<X>(_._));
0152
0153   //  /
0154   //  | f * g dx : defined only for f or g indep. of x or part. int.
0155   //  /
0156
0157   template <TYPXT2,
0158             bool bint = not ::std::is_same<PRIMIT(X, B), UndefinedIntegral>::value,
0159             bool aint = not ::std::is_same<PRIMIT(X, A), UndefinedIntegral>::value>
0160   struct PartIntegral {
0161     typedef UndefinedIntegral type;
0162     GET(PROD_S(A, B), type());
0163   };
0164
0165   TEMPL(XT2) struct PartIntegral<X, A, B, true, false> {
0166     typedef PRIMIT(X, B) B1;
0167     typedef DERIV(X, A) A1;
0168     typedef PRIMIT(X, PROD(A1, B1)) AB1;
0169     typedef DIFF(PROD(A, B1), PRIMIT(X, PROD(A1, B1))) type;
0170     inline static type get(const PROD_S(A, B) & _) {
0171       const A& a = _._1;
0172       B1 b = primitive<X>(_._2);
0173       return a * b - primitive<X>(derivative<X>(a) * b);
0174     }
0175   };
0176
0177   TEMPL(XT2) struct PartIntegral<X, B, A, false, true> {
0178     typedef PRIMIT(X, B) B1;
0179     typedef DERIV(X, A) A1;
0180     typedef PRIMIT(X, PROD(A1, B1)) AB1;
0181     typedef DIFF(PROD(A, B1), PRIMIT(X, PROD(A1, B1))) type;
0182     inline static type get(const PROD_S(B, A) & _) {
0183       const A& a = _._2;
0184       B1 b = primitive<X>(_._1);
0185       return a * b - primitive<X>(derivative<X>(a) * b);
0186     }
0187   };
0188
0189   TEMPL(XT2) struct PartIntegral<X, A, B, true, true> : public PartIntegral<X, A, B, true, false> {};
0190
0191   template <TYPXT2, bool indepf = Independent<X, A>::value, bool indepg = Independent<X, B>::value>
0192   struct ProductPrimitive : public PartIntegral<X, A, B> {};
0193
0194   TEMPL(XT2) struct ProductPrimitive<X, A, B, true, false> {
0195     typedef PROD(A, PRIMIT(X, B)) type;
0196     GET(PROD_S(A, B), _._1* primitive<X>(_._2));
0197   };
0198
0199   TEMPL(XT2) struct ProductPrimitive<X, A, B, false, true> {
0200     typedef PROD(B, PRIMIT(X, A)) type;
0201     GET(PROD_S(A, B), _._2* primitive<X>(_._1));
0202   };
0203
0204   TEMPL(XT2) struct ProductPrimitive<X, A, B, true, true> {
0205     typedef PROD(PROD(A, B), X) type;
0206     GET(PROD_S(A, B), _* X());
0207   };
0208
0209   TEMPL(XT2) struct Primitive<PROD_S(A, B), X> : public ProductPrimitive<X, A, B> {};
0210
0211   //  /
0212   //  | f / g dx : defined only for f or g indep. of x; try part. int.
0213   //  /
0214
0215   template <TYPXT2,
0216             bool bint = not ::std::is_same<PRIMIT(X, RATIO(NUM(1), B)), UndefinedIntegral>::value,
0217             bool aint = not ::std::is_same<PRIMIT(X, A), UndefinedIntegral>::value>
0218   struct PartIntegral2 {
0219     typedef UndefinedIntegral type;
0220     GET(RATIO_S(A, B), type());
0221   };
0222
0223   TEMPL(XT2) struct PartIntegral2<X, A, B, true, false> {
0224     typedef PRIMIT(X, RATIO(NUM(1), B)) B1;
0225     typedef DERIV(X, A) A1;
0226     typedef PRIMIT(X, PROD(A1, B1)) AB1;
0227     typedef DIFF(PROD(A, B1), AB1) type;
0228     inline static type get(const RATIO_S(A, B) & _) {
0229       const A& a = _._1;
0230       B1 b = primitive<X>(num<1>() / _._2);
0231       return a * b - primitive<X>(derivative<X>(a) * b);
0232     }
0233   };
0234
0235   TEMPL(XT2) struct PartIntegral2<X, B, A, false, true> {
0236     typedef PRIMIT(X, RATIO(NUM(1), B)) B1;
0237     typedef DERIV(X, A) A1;
0238     typedef PRIMIT(X, PROD(A1, B1)) AB1;
0239     typedef DIFF(PROD(A, B1), AB1) type;
0240     inline static type get(const RATIO_S(B, A) & _) {
0241       const A& a = _._1;
0242       B1 b = primitive<X>(num<1>() / _._2);
0243       return a * b - primitive<X>(derivative<X>(a) * b);
0244     }
0245   };
0246
0247   // should be improved: try both...
0248   TEMPL(XT2) struct PartIntegral2<X, A, B, true, true> : public PartIntegral2<X, A, B, true, false> {};
0249
0250   template <TYPXT2, bool indepa = Independent<X, A>::value, bool indepb = Independent<X, B>::value>
0251   struct RatioPrimitive : public PartIntegral2<X, A, B> {};
0252
0253   TEMPL(XT2) struct RatioPrimitive<X, A, B, true, false> {
0254     typedef PROD(A, PRIMIT(X, RATIO(NUM(1), B))) type;
0255     GET(RATIO_S(A, B), _._1* primitive<X>(num<1> / _._2));
0256   };
0257
0258   TEMPL(XT2) struct RatioPrimitive<X, A, B, false, true> {
0259     typedef RATIO(PRIMIT(X, A), B) type;
0260     GET(RATIO_S(A, B), primitive<X>(_._1) / _._2);
0261   };
0262
0263   TEMPL(XT2) struct RatioPrimitive<X, A, B, true, true> {
0264     typedef RATIO(RATIO(A, B), X) type;
0265     GET(RATIO_S(A, B), _* X());
0266   };
0267
0268   TEMPL(XT2) struct Primitive<RATIO_S(A, B), X> : public RatioPrimitive<X, A, B> {};
0269
0270   // Function integrals
0271
0272   //  /
0273   //  | c dx = c * x
0274   //  /
0275   template <>
0276   struct Primitive<Parameter> {
0277     typedef Product<Parameter, Identity>::type type;
0278     inline static type get(const Parameter& p) { return p * Identity(); }
0279   };
0280
0281 }  // namespace funct
0282
0283 #define DECLARE_PRIMITIVE(X, F, P)                             \
0284   namespace funct {                                            \
0285     template <typename X>                                      \
0286     struct Primitive<F, X> {                                   \
0287       typedef P type;                                          \
0288       inline static type get(const F& _) { return type(_._); } \
0289     };                                                         \
0290   }                                                            \
0291   struct __useless_ignoreme
0292
0293 #define DECLARE_FUNCT_PRIMITIVE(F, P)                       \
0294   namespace funct {                                         \
0295     template <>                                             \
0296     struct Primitive<F> {                                   \
0297       typedef P type;                                       \
0298       inline static type get(const F& _) { return type(); } \
0299     };                                                      \
0300   }                                                         \
0301   struct __useless_ignoreme
0302
0303 #include "PhysicsTools/Utilities/interface/Simplify_end.h"
0304
0305 #endif