Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 10:22:02 

0001 /// \file CladDerivator.h
0002 ///
0003 /// \brief The file is a bridge between ROOT and clad automatic differentiation
0004 /// plugin.
0005 ///
0006 /// \author Vassil Vassilev <vvasilev@cern.ch>
0007 ///
0008 /// \date July, 2018
0009 
0010 /*************************************************************************
0011  * Copyright (C) 1995-2018, Rene Brun and Fons Rademakers.               *
0012  * All rights reserved.                                                  *
0013  *                                                                       *
0014  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0015  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0016  *************************************************************************/
0017 
0018 #ifndef CLAD_DERIVATOR
0019 #define CLAD_DERIVATOR
0020 
0021 #ifndef __CLING__
0022 #error "This file must not be included by compiled programs."
0023 #endif //__CLING__
0024 
0025 #include <plugins/include/clad/Differentiator/Differentiator.h>
0026 #include "TMath.h"
0027 namespace clad {
0028 namespace custom_derivatives {
0029 namespace TMath {
0030 template <typename T>
0031 ValueAndPushforward<T, T> Abs_pushforward(T x, T d_x)
0032 {
0033    return {::TMath::Abs(x), ((x < 0) ? -1 : 1) * d_x};
0034 }
0035 
0036 template <typename T>
0037 ValueAndPushforward<T, T> ACos_pushforward(T x, T d_x)
0038 {
0039    return {::TMath::ACos(x), (-1. / ::TMath::Sqrt(1 - x * x)) * d_x};
0040 }
0041 
0042 template <typename T>
0043 ValueAndPushforward<T, T> ACosH_pushforward(T x, T d_x)
0044 {
0045    return {::TMath::ACosH(x), (1. / ::TMath::Sqrt(x * x - 1)) * d_x};
0046 }
0047 
0048 template <typename T>
0049 ValueAndPushforward<T, T> ASin_pushforward(T x, T d_x)
0050 {
0051    return {::TMath::ASin(x), (1. / ::TMath::Sqrt(1 - x * x)) * d_x};
0052 }
0053 
0054 template <typename T>
0055 ValueAndPushforward<T, T> ASinH_pushforward(T x, T d_x)
0056 {
0057    return {::TMath::ASinH(x), (1. / ::TMath::Sqrt(x * x + 1)) * d_x};
0058 }
0059 
0060 template <typename T>
0061 ValueAndPushforward<T, T> ATan_pushforward(T x, T d_x)
0062 {
0063    return {::TMath::ATan(x), (1. / (x * x + 1)) * d_x};
0064 }
0065 
0066 template <typename T>
0067 ValueAndPushforward<T, T> ATanH_pushforward(T x, T d_x)
0068 {
0069    return {::TMath::ATanH(x), (1. / (1 - x * x)) * d_x};
0070 }
0071 
0072 template <typename T>
0073 ValueAndPushforward<T, T> Cos_pushforward(T x, T d_x)
0074 {
0075    return {::TMath::Cos(x), -::TMath::Sin(x) * d_x};
0076 }
0077 
0078 template <typename T>
0079 ValueAndPushforward<T, T> CosH_pushforward(T x, T d_x)
0080 {
0081    return {::TMath::CosH(x), ::TMath::SinH(x) * d_x};
0082 }
0083 
0084 template <typename T>
0085 ValueAndPushforward<T, T> Erf_pushforward(T x, T d_x)
0086 {
0087    return {::TMath::Erf(x), (2 * ::TMath::Exp(-x * x) / ::TMath::Sqrt(::TMath::Pi())) * d_x};
0088 }
0089 
0090 template <typename T>
0091 ValueAndPushforward<T, T> Erfc_pushforward(T x, T d_x)
0092 {
0093    return {::TMath::Erfc(x), -Erf_pushforward(x, d_x).pushforward};
0094 }
0095 
0096 template <typename T>
0097 ValueAndPushforward<T, T> Exp_pushforward(T x, T d_x)
0098 {
0099    return {::TMath::Exp(x), ::TMath::Exp(x) * d_x};
0100 }
0101 
0102 template <typename T>
0103 ValueAndPushforward<T, T> Hypot_pushforward(T x, T y, T d_x, T d_y)
0104 {
0105    return {::TMath::Hypot(x, y), x / ::TMath::Hypot(x, y) * d_x + y / ::TMath::Hypot(x, y) * d_y};
0106 }
0107 
0108 template <typename T, typename U>
0109 void Hypot_pullback(T x, T y, U p, clad::array_ref<T> d_x, clad::array_ref<T> d_y)
0110 {
0111    T h = ::TMath::Hypot(x, y);
0112    *d_x += x / h * p;
0113    *d_y += y / h * p;
0114 }
0115 
0116 template <typename T>
0117 ValueAndPushforward<T, T> Log_pushforward(T x, T d_x)
0118 {
0119    return {::TMath::Log(x), (1. / x) * d_x};
0120 }
0121 
0122 template <typename T>
0123 ValueAndPushforward<T, T> Log10_pushforward(T x, T d_x)
0124 {
0125    return {::TMath::Log10(x), (1.0 / (x * ::TMath::Ln10())) * d_x};
0126 }
0127 
0128 template <typename T>
0129 ValueAndPushforward<T, T> Log2_pushforward(T x, T d_x)
0130 {
0131    return {::TMath::Log2(x), (1.0 / (x * ::TMath::Log(2.0))) * d_x};
0132 }
0133 
0134 template <typename T>
0135 ValueAndPushforward<T, T> Max_pushforward(T x, T y, T d_x, T d_y)
0136 {
0137    T derivative = 0;
0138    if (x >= y)
0139       derivative = d_x;
0140    else
0141       derivative = d_y;
0142    return {::TMath::Max(x, y), derivative};
0143 }
0144 
0145 template <typename T, typename U>
0146 void Max_pullback(T a, T b, U p, clad::array_ref<T> d_a, clad::array_ref<T> d_b)
0147 {
0148    if (a >= b)
0149       *d_a += p;
0150    else
0151       *d_b += p;
0152 }
0153 
0154 template <typename T>
0155 ValueAndPushforward<T, T> Min_pushforward(T x, T y, T d_x, T d_y)
0156 {
0157    T derivative = 0;
0158    if (x <= y)
0159       derivative = d_x;
0160    else
0161       derivative = d_y;
0162    return {::TMath::Min(x, y), derivative};
0163 }
0164 
0165 template <typename T, typename U>
0166 void Min_pullback(T a, T b, U p, clad::array_ref<T> d_a, clad::array_ref<T> d_b)
0167 {
0168    if (a <= b)
0169       *d_a += p;
0170    else
0171       *d_b += p;
0172 }
0173 
0174 template <typename T>
0175 ValueAndPushforward<T, T> Power_pushforward(T x, T y, T d_x, T d_y)
0176 {
0177    T pushforward = y * ::TMath::Power(x, y - 1) * d_x;
0178    if (d_y) {
0179       pushforward += (::TMath::Power(x, y) * ::TMath::Log(x)) * d_y;
0180    }
0181    return {::TMath::Power(x, y), pushforward};
0182 }
0183 
0184 template <typename T, typename U>
0185 void Power_pullback(T x, T y, U p, clad::array_ref<T> d_x, clad::array_ref<T> d_y)
0186 {
0187    auto t = pow_pushforward(x, y, 1, 0);
0188    *d_x += t.pushforward * p;
0189    t = pow_pushforward(x, y, 0, 1);
0190    *d_y += t.pushforward * p;
0191 }
0192 
0193 template <typename T>
0194 ValueAndPushforward<T, T> Sin_pushforward(T x, T d_x)
0195 {
0196    return {::TMath::Sin(x), ::TMath::Cos(x) * d_x};
0197 }
0198 
0199 template <typename T>
0200 ValueAndPushforward<T, T> SinH_pushforward(T x, T d_x)
0201 {
0202    return {::TMath::SinH(x), ::TMath::CosH(x) * d_x};
0203 }
0204 
0205 template <typename T>
0206 ValueAndPushforward<T, T> Sq_pushforward(T x, T d_x)
0207 {
0208    return {::TMath::Sq(x), 2 * x * d_x};
0209 }
0210 
0211 template <typename T>
0212 ValueAndPushforward<T, T> Sqrt_pushforward(T x, T d_x)
0213 {
0214    return {::TMath::Sqrt(x), (0.5 / ::TMath::Sqrt(x)) * d_x};
0215 }
0216 
0217 template <typename T>
0218 ValueAndPushforward<T, T> Tan_pushforward(T x, T d_x)
0219 {
0220    return {::TMath::Tan(x), (1. / ::TMath::Sq(::TMath::Cos(x))) * d_x};
0221 }
0222 
0223 template <typename T>
0224 ValueAndPushforward<T, T> TanH_pushforward(T x, T d_x)
0225 {
0226    return {::TMath::TanH(x), (1. / ::TMath::Sq(::TMath::CosH(x))) * d_x};
0227 }
0228 
0229 #ifdef WIN32
0230 // Additional custom derivatives that can be removed
0231 // after Issue #12108 in ROOT is resolved
0232 // constexpr is removed
0233 ValueAndPushforward<Double_t, Double_t> Pi_pushforward()
0234 {
0235    return {3.1415926535897931, 0.};
0236 }
0237 // constexpr is removed
0238 ValueAndPushforward<Double_t, Double_t> Ln10_pushforward()
0239 {
0240    return {2.3025850929940459, 0.};
0241 }
0242 #endif
0243 } // namespace TMath
0244 
0245 
0246 namespace ROOT {
0247 namespace Math {
0248 
0249 inline void landau_pdf_pullback(double x, double xi, double x0, double d_out, double *d_x, double *d_xi, double *d_x0)
0250 {
0251    if (xi <= 0) {
0252       return;
0253    }
0254    // clang-format off
0255    static double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635,   0.001511162253};
0256    static double q1[5] = {1.0         ,-0.3388260629, 0.09594393323, -0.01608042283,    0.003778942063};
0257 
0258    static double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411,   0.0001283617211};
0259    static double q2[5] = {1.0         , 0.7428795082, 0.3153932961,   0.06694219548,    0.008790609714};
0260 
0261    static double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
0262    static double q3[5] = {1.0         , 0.6097809921, 0.2560616665,   0.04746722384,    0.006957301675};
0263 
0264    static double p4[5] = {0.9874054407, 118.6723273,  849.2794360,   -743.7792444,      427.0262186};
0265    static double q4[5] = {1.0         , 106.8615961,  337.6496214,    2016.712389,      1597.063511};
0266 
0267    static double p5[5] = {1.003675074,  167.5702434,  4789.711289,    21217.86767,     -22324.94910};
0268    static double q5[5] = {1.0         , 156.9424537,  3745.310488,    9834.698876,      66924.28357};
0269 
0270    static double p6[5] = {1.000827619,  664.9143136,  62972.92665,    475554.6998,     -5743609.109};
0271    static double q6[5] = {1.0         , 651.4101098,  56974.73333,    165917.4725,     -2815759.939};
0272 
0273    static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
0274 
0275    static double a2[2] = {-1.845568670,-4.284640743};
0276    // clang-format on
0277    const double _const0 = 0.3989422803;
0278    double v = (x - x0) / xi;
0279    double _d_v = 0;
0280    double _d_denlan = 0;
0281    if (v < -5.5) {
0282       double _t0;
0283       double u = ::std::exp(v + 1.);
0284       double _d_u = 0;
0285       if (u >= 1.e-10) {
0286          const double ue = ::std::exp(-1 / u);
0287          const double us = ::std::sqrt(u);
0288          double _t3;
0289          double _d_ue = 0;
0290          double _d_us = 0;
0291          double denlan = _const0 * (ue / us) * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
0292          _d_denlan += d_out / xi;
0293          *d_xi += d_out * -(denlan / (xi * xi));
0294          denlan = _t3;
0295          double _r_d3 = _d_denlan;
0296          _d_denlan -= _r_d3;
0297          _d_ue += _const0 * _r_d3 * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u) / us;
0298          double _r5 = _const0 * _r_d3 * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u) * -(ue / (us * us));
0299          _d_us += _r5;
0300          _d_u += a1[2] * _const0 * (ue / us) * _r_d3 * u * u;
0301          _d_u += (a1[1] + a1[2] * u) * _const0 * (ue / us) * _r_d3 * u;
0302          _d_u += (a1[0] + (a1[1] + a1[2] * u) * u) * _const0 * (ue / us) * _r_d3;
0303          double _r_d2 = _d_us;
0304          _d_us -= _r_d2;
0305          double _r4 = 0;
0306          _r4 += _r_d2 * clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0307          _d_u += _r4;
0308          double _r_d1 = _d_ue;
0309          _d_ue -= _r_d1;
0310          double _r2 = 0;
0311          _r2 += _r_d1 * clad::custom_derivatives::exp_pushforward(-1 / u, 1.).pushforward;
0312          double _r3 = _r2 * -(-1 / (u * u));
0313          _d_u += _r3;
0314       }
0315       u = _t0;
0316       double _r_d0 = _d_u;
0317       _d_u -= _r_d0;
0318       double _r1 = 0;
0319       _r1 += _r_d0 * clad::custom_derivatives::exp_pushforward(v + 1., 1.).pushforward;
0320       _d_v += _r1;
0321    } else if (v < -1) {
0322       double _t4;
0323       double u = ::std::exp(-v - 1);
0324       double _d_u = 0;
0325       double _t5;
0326       double _t8 = ::std::exp(-u);
0327       double _t7 = ::std::sqrt(u);
0328       double _t6 = (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
0329       double denlan = _t8 * _t7 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / _t6;
0330       _d_denlan += d_out / xi;
0331       *d_xi += d_out * -(denlan / (xi * xi));
0332       denlan = _t5;
0333       double _r_d5 = _d_denlan;
0334       _d_denlan -= _r_d5;
0335       double _r7 = 0;
0336       _r7 += _r_d5 / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) * _t7 *
0337              clad::custom_derivatives::exp_pushforward(-u, 1.).pushforward;
0338       _d_u += -_r7;
0339       double _r8 = 0;
0340       _r8 += _t8 * _r_d5 / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) *
0341              clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0342       _d_u += _r8;
0343       _d_v += p1[4] * _t8 * _t7 * _r_d5 / _t6 * v * v * v;
0344       _d_v += (p1[3] + p1[4] * v) * _t8 * _t7 * _r_d5 / _t6 * v * v;
0345       _d_v += (p1[2] + (p1[3] + p1[4] * v) * v) * _t8 * _t7 * _r_d5 / _t6 * v;
0346       _d_v += (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * _t8 * _t7 * _r_d5 / _t6;
0347       double _r9 = _r_d5 * -(_t8 * _t7 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (_t6 * _t6));
0348       _d_v += q1[4] * _r9 * v * v * v;
0349       _d_v += (q1[3] + q1[4] * v) * _r9 * v * v;
0350       _d_v += (q1[2] + (q1[3] + q1[4] * v) * v) * _r9 * v;
0351       _d_v += (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * _r9;
0352       u = _t4;
0353       double _r_d4 = _d_u;
0354       _d_u -= _r_d4;
0355       double _r6 = 0;
0356       _r6 += _r_d4 * clad::custom_derivatives::exp_pushforward(-v - 1, 1.).pushforward;
0357       _d_v += -_r6;
0358    } else if (v < 1) {
0359       double _t9;
0360       double _t10 = (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
0361       double denlan = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / _t10;
0362       _d_denlan += d_out / xi;
0363       *d_xi += d_out * -(denlan / (xi * xi));
0364       denlan = _t9;
0365       double _r_d6 = _d_denlan;
0366       _d_denlan -= _r_d6;
0367       _d_v += p2[4] * _r_d6 / _t10 * v * v * v;
0368       _d_v += (p2[3] + p2[4] * v) * _r_d6 / _t10 * v * v;
0369       _d_v += (p2[2] + (p2[3] + p2[4] * v) * v) * _r_d6 / _t10 * v;
0370       _d_v += (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * _r_d6 / _t10;
0371       double _r10 = _r_d6 * -((p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / (_t10 * _t10));
0372       _d_v += q2[4] * _r10 * v * v * v;
0373       _d_v += (q2[3] + q2[4] * v) * _r10 * v * v;
0374       _d_v += (q2[2] + (q2[3] + q2[4] * v) * v) * _r10 * v;
0375       _d_v += (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * _r10;
0376    } else if (v < 5) {
0377       double _t11;
0378       double _t12 = (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
0379       double denlan = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / _t12;
0380       _d_denlan += d_out / xi;
0381       *d_xi += d_out * -(denlan / (xi * xi));
0382       denlan = _t11;
0383       double _r_d7 = _d_denlan;
0384       _d_denlan -= _r_d7;
0385       _d_v += p3[4] * _r_d7 / _t12 * v * v * v;
0386       _d_v += (p3[3] + p3[4] * v) * _r_d7 / _t12 * v * v;
0387       _d_v += (p3[2] + (p3[3] + p3[4] * v) * v) * _r_d7 / _t12 * v;
0388       _d_v += (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * _r_d7 / _t12;
0389       double _r11 = _r_d7 * -((p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / (_t12 * _t12));
0390       _d_v += q3[4] * _r11 * v * v * v;
0391       _d_v += (q3[3] + q3[4] * v) * _r11 * v * v;
0392       _d_v += (q3[2] + (q3[3] + q3[4] * v) * v) * _r11 * v;
0393       _d_v += (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * _r11;
0394    } else if (v < 12) {
0395       double u = 1 / v;
0396       double _d_u = 0;
0397       double _t14;
0398       double _t15 = (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
0399       double denlan = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / _t15;
0400       _d_denlan += d_out / xi;
0401       *d_xi += d_out * -(denlan / (xi * xi));
0402       denlan = _t14;
0403       double _r_d9 = _d_denlan;
0404       _d_denlan -= _r_d9;
0405       _d_u += _r_d9 / _t15 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) * u;
0406       _d_u += u * _r_d9 / _t15 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u);
0407       _d_u += p4[4] * u * u * _r_d9 / _t15 * u * u * u;
0408       _d_u += (p4[3] + p4[4] * u) * u * u * _r_d9 / _t15 * u * u;
0409       _d_u += (p4[2] + (p4[3] + p4[4] * u) * u) * u * u * _r_d9 / _t15 * u;
0410       _d_u += (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u * u * _r_d9 / _t15;
0411       double _r13 = _r_d9 * -(u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / (_t15 * _t15));
0412       _d_u += q4[4] * _r13 * u * u * u;
0413       _d_u += (q4[3] + q4[4] * u) * _r13 * u * u;
0414       _d_u += (q4[2] + (q4[3] + q4[4] * u) * u) * _r13 * u;
0415       _d_u += (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * _r13;
0416       double _r_d8 = _d_u;
0417       _d_u -= _r_d8;
0418       double _r12 = _r_d8 * -(1 / (v * v));
0419       _d_v += _r12;
0420    } else if (v < 50) {
0421       double u = 1 / v;
0422       double _d_u = 0;
0423       double _t17;
0424       double _t18 = (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
0425       double denlan = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / _t18;
0426       _d_denlan += d_out / xi;
0427       *d_xi += d_out * -(denlan / (xi * xi));
0428       denlan = _t17;
0429       double _r_d11 = _d_denlan;
0430       _d_denlan -= _r_d11;
0431       _d_u += _r_d11 / _t18 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) * u;
0432       _d_u += u * _r_d11 / _t18 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u);
0433       _d_u += p5[4] * u * u * _r_d11 / _t18 * u * u * u;
0434       _d_u += (p5[3] + p5[4] * u) * u * u * _r_d11 / _t18 * u * u;
0435       _d_u += (p5[2] + (p5[3] + p5[4] * u) * u) * u * u * _r_d11 / _t18 * u;
0436       _d_u += (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u * u * _r_d11 / _t18;
0437       double _r15 = _r_d11 * -(u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / (_t18 * _t18));
0438       _d_u += q5[4] * _r15 * u * u * u;
0439       _d_u += (q5[3] + q5[4] * u) * _r15 * u * u;
0440       _d_u += (q5[2] + (q5[3] + q5[4] * u) * u) * _r15 * u;
0441       _d_u += (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * _r15;
0442       double _r_d10 = _d_u;
0443       _d_u -= _r_d10;
0444       double _r14 = _r_d10 * -(1 / (v * v));
0445       _d_v += _r14;
0446    } else if (v < 300) {
0447       double _t19;
0448       double u = 1 / v;
0449       double _d_u = 0;
0450       double _t20;
0451       double _t21 = (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
0452       double denlan = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / _t21;
0453       _d_denlan += d_out / xi;
0454       *d_xi += d_out * -(denlan / (xi * xi));
0455       denlan = _t20;
0456       double _r_d13 = _d_denlan;
0457       _d_denlan -= _r_d13;
0458       _d_u += _r_d13 / _t21 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) * u;
0459       _d_u += u * _r_d13 / _t21 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u);
0460       _d_u += p6[4] * u * u * _r_d13 / _t21 * u * u * u;
0461       _d_u += (p6[3] + p6[4] * u) * u * u * _r_d13 / _t21 * u * u;
0462       _d_u += (p6[2] + (p6[3] + p6[4] * u) * u) * u * u * _r_d13 / _t21 * u;
0463       _d_u += (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u * u * _r_d13 / _t21;
0464       double _r17 = _r_d13 * -(u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / (_t21 * _t21));
0465       _d_u += q6[4] * _r17 * u * u * u;
0466       _d_u += (q6[3] + q6[4] * u) * _r17 * u * u;
0467       _d_u += (q6[2] + (q6[3] + q6[4] * u) * u) * _r17 * u;
0468       _d_u += (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * _r17;
0469       u = _t19;
0470       double _r_d12 = _d_u;
0471       _d_u -= _r_d12;
0472       double _r16 = _r_d12 * -(1 / (v * v));
0473       _d_v += _r16;
0474    } else {
0475       double _t22;
0476       double _t25 = ::std::log(v);
0477       double _t24 = (v + 1);
0478       double _t23 = (v - v * _t25 / _t24);
0479       double u = 1 / _t23;
0480       double _d_u = 0;
0481       double _t26;
0482       double denlan = u * u * (1 + (a2[0] + a2[1] * u) * u);
0483       _d_denlan += d_out / xi;
0484       *d_xi += d_out * -(denlan / (xi * xi));
0485       denlan = _t26;
0486       double _r_d15 = _d_denlan;
0487       _d_denlan -= _r_d15;
0488       _d_u += _r_d15 * (1 + (a2[0] + a2[1] * u) * u) * u;
0489       _d_u += u * _r_d15 * (1 + (a2[0] + a2[1] * u) * u);
0490       _d_u += a2[1] * u * u * _r_d15 * u;
0491       _d_u += (a2[0] + a2[1] * u) * u * u * _r_d15;
0492       u = _t22;
0493       double _r_d14 = _d_u;
0494       _d_u -= _r_d14;
0495       double _r18 = _r_d14 * -(1 / (_t23 * _t23));
0496       _d_v += _r18;
0497       _d_v += -_r18 / _t24 * _t25;
0498       double _r19 = 0;
0499       _r19 += v * -_r18 / _t24 * clad::custom_derivatives::log_pushforward(v, 1.).pushforward;
0500       _d_v += _r19;
0501       double _r20 = -_r18 * -(v * _t25 / (_t24 * _t24));
0502       _d_v += _r20;
0503    }
0504    *d_x += _d_v / xi;
0505    *d_x0 += -_d_v / xi;
0506    double _r0 = _d_v * -((x - x0) / (xi * xi));
0507    *d_xi += _r0;
0508 }
0509 
0510 inline void landau_cdf_pullback(double x, double xi, double x0, double d_out, double *d_x, double *d_xi, double *d_x0)
0511 {
0512    // clang-format off
0513    static double p1[5] = {0.2514091491e+0,-0.6250580444e-1, 0.1458381230e-1,-0.2108817737e-2, 0.7411247290e-3};
0514    static double q1[5] = {1.0            ,-0.5571175625e-2, 0.6225310236e-1,-0.3137378427e-2, 0.1931496439e-2};
0515 
0516    static double p2[4] = {0.2868328584e+0, 0.3564363231e+0, 0.1523518695e+0, 0.2251304883e-1};
0517    static double q2[4] = {1.0            , 0.6191136137e+0, 0.1720721448e+0, 0.2278594771e-1};
0518 
0519    static double p3[4] = {0.2868329066e+0, 0.3003828436e+0, 0.9950951941e-1, 0.8733827185e-2};
0520    static double q3[4] = {1.0            , 0.4237190502e+0, 0.1095631512e+0, 0.8693851567e-2};
0521 
0522    static double p4[4] = {0.1000351630e+1, 0.4503592498e+1, 0.1085883880e+2, 0.7536052269e+1};
0523    static double q4[4] = {1.0            , 0.5539969678e+1, 0.1933581111e+2, 0.2721321508e+2};
0524 
0525    static double p5[4] = {0.1000006517e+1, 0.4909414111e+2, 0.8505544753e+2, 0.1532153455e+3};
0526    static double q5[4] = {1.0            , 0.5009928881e+2, 0.1399819104e+3, 0.4200002909e+3};
0527 
0528    static double p6[4] = {0.1000000983e+1, 0.1329868456e+3, 0.9162149244e+3,-0.9605054274e+3};
0529    static double q6[4] = {1.0            , 0.1339887843e+3, 0.1055990413e+4, 0.5532224619e+3};
0530 
0531    static double a1[4] = {0              ,-0.4583333333e+0, 0.6675347222e+0,-0.1641741416e+1};
0532    static double a2[4] = {0              , 1.0            ,-0.4227843351e+0,-0.2043403138e+1};
0533    // clang-format on
0534 
0535    const double v = (x - x0) / xi;
0536    double _d_v = 0;
0537    if (v < -5.5) {
0538       double _d_u = 0;
0539       const double _const0 = 0.3989422803;
0540       double u = ::std::exp(v + 1);
0541       double _t3 = ::std::exp(-1. / u);
0542       double _t2 = ::std::sqrt(u);
0543       double _r2 = 0;
0544       _r2 += _const0 * d_out * (1 + (a1[1] + (a1[2] + a1[3] * u) * u) * u) * _t2 *
0545              clad::custom_derivatives::exp_pushforward(-1. / u, 1.).pushforward;
0546       double _r3 = _r2 * -(-1. / (u * u));
0547       _d_u += _r3;
0548       double _r4 = 0;
0549       _r4 += _const0 * _t3 * d_out * (1 + (a1[1] + (a1[2] + a1[3] * u) * u) * u) *
0550              clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0551       _d_u += _r4;
0552       _d_u += a1[3] * _const0 * _t3 * _t2 * d_out * u * u;
0553       _d_u += (a1[2] + a1[3] * u) * _const0 * _t3 * _t2 * d_out * u;
0554       _d_u += (a1[1] + (a1[2] + a1[3] * u) * u) * _const0 * _t3 * _t2 * d_out;
0555       _d_v += _d_u * clad::custom_derivatives::exp_pushforward(v + 1, 1.).pushforward;
0556    } else if (v < -1) {
0557       double _d_u = 0;
0558       double u = ::std::exp(-v - 1);
0559       double _t8 = ::std::exp(-u);
0560       double _t7 = ::std::sqrt(u);
0561       double _t6 = (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
0562       double _r6 = 0;
0563       _r6 += d_out / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / _t7 *
0564              clad::custom_derivatives::exp_pushforward(-u, 1.).pushforward;
0565       _d_u += -_r6;
0566       double _r7 = d_out / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) * -(_t8 / (_t7 * _t7));
0567       double _r8 = 0;
0568       _r8 += _r7 * clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0569       _d_u += _r8;
0570       _d_v += p1[4] * (_t8 / _t7) * d_out / _t6 * v * v * v;
0571       _d_v += (p1[3] + p1[4] * v) * (_t8 / _t7) * d_out / _t6 * v * v;
0572       _d_v += (p1[2] + (p1[3] + p1[4] * v) * v) * (_t8 / _t7) * d_out / _t6 * v;
0573       _d_v += (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * (_t8 / _t7) * d_out / _t6;
0574       double _r9 = d_out * -((_t8 / _t7) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (_t6 * _t6));
0575       _d_v += q1[4] * _r9 * v * v * v;
0576       _d_v += (q1[3] + q1[4] * v) * _r9 * v * v;
0577       _d_v += (q1[2] + (q1[3] + q1[4] * v) * v) * _r9 * v;
0578       _d_v += (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * _r9;
0579       _d_v += -_d_u * clad::custom_derivatives::exp_pushforward(-v - 1, 1.).pushforward;
0580    } else if (v < 1) {
0581       double _t10 = (q2[0] + (q2[1] + (q2[2] + q2[3] * v) * v) * v);
0582       _d_v += p2[3] * d_out / _t10 * v * v;
0583       _d_v += (p2[2] + p2[3] * v) * d_out / _t10 * v;
0584       _d_v += (p2[1] + (p2[2] + p2[3] * v) * v) * d_out / _t10;
0585       double _r10 = d_out * -((p2[0] + (p2[1] + (p2[2] + p2[3] * v) * v) * v) / (_t10 * _t10));
0586       _d_v += q2[3] * _r10 * v * v;
0587       _d_v += (q2[2] + q2[3] * v) * _r10 * v;
0588       _d_v += (q2[1] + (q2[2] + q2[3] * v) * v) * _r10;
0589    } else if (v < 4) {
0590       double _t12 = (q3[0] + (q3[1] + (q3[2] + q3[3] * v) * v) * v);
0591       _d_v += p3[3] * d_out / _t12 * v * v;
0592       _d_v += (p3[2] + p3[3] * v) * d_out / _t12 * v;
0593       _d_v += (p3[1] + (p3[2] + p3[3] * v) * v) * d_out / _t12;
0594       double _r11 = d_out * -((p3[0] + (p3[1] + (p3[2] + p3[3] * v) * v) * v) / (_t12 * _t12));
0595       _d_v += q3[3] * _r11 * v * v;
0596       _d_v += (q3[2] + q3[3] * v) * _r11 * v;
0597       _d_v += (q3[1] + (q3[2] + q3[3] * v) * v) * _r11;
0598    } else if (v < 12) {
0599       double _d_u = 0;
0600       double u = 1. / v;
0601       double _t15 = (q4[0] + (q4[1] + (q4[2] + q4[3] * u) * u) * u);
0602       _d_u += p4[3] * d_out / _t15 * u * u;
0603       _d_u += (p4[2] + p4[3] * u) * d_out / _t15 * u;
0604       _d_u += (p4[1] + (p4[2] + p4[3] * u) * u) * d_out / _t15;
0605       double _r13 = d_out * -((p4[0] + (p4[1] + (p4[2] + p4[3] * u) * u) * u) / (_t15 * _t15));
0606       _d_u += q4[3] * _r13 * u * u;
0607       _d_u += (q4[2] + q4[3] * u) * _r13 * u;
0608       _d_u += (q4[1] + (q4[2] + q4[3] * u) * u) * _r13;
0609       double _r12 = _d_u * -(1. / (v * v));
0610       _d_v += _r12;
0611    } else if (v < 50) {
0612       double _d_u = 0;
0613       double u = 1. / v;
0614       double _t18 = (q5[0] + (q5[1] + (q5[2] + q5[3] * u) * u) * u);
0615       _d_u += p5[3] * d_out / _t18 * u * u;
0616       _d_u += (p5[2] + p5[3] * u) * d_out / _t18 * u;
0617       _d_u += (p5[1] + (p5[2] + p5[3] * u) * u) * d_out / _t18;
0618       double _r15 = d_out * -((p5[0] + (p5[1] + (p5[2] + p5[3] * u) * u) * u) / (_t18 * _t18));
0619       _d_u += q5[3] * _r15 * u * u;
0620       _d_u += (q5[2] + q5[3] * u) * _r15 * u;
0621       _d_u += (q5[1] + (q5[2] + q5[3] * u) * u) * _r15;
0622       double _r14 = _d_u * -(1. / (v * v));
0623       _d_v += _r14;
0624    } else if (v < 300) {
0625       double _d_u = 0;
0626       double u = 1. / v;
0627       double _t21 = (q6[0] + (q6[1] + (q6[2] + q6[3] * u) * u) * u);
0628       _d_u += p6[3] * d_out / _t21 * u * u;
0629       _d_u += (p6[2] + p6[3] * u) * d_out / _t21 * u;
0630       _d_u += (p6[1] + (p6[2] + p6[3] * u) * u) * d_out / _t21;
0631       double _r17 = d_out * -((p6[0] + (p6[1] + (p6[2] + p6[3] * u) * u) * u) / (_t21 * _t21));
0632       _d_u += q6[3] * _r17 * u * u;
0633       _d_u += (q6[2] + q6[3] * u) * _r17 * u;
0634       _d_u += (q6[1] + (q6[2] + q6[3] * u) * u) * _r17;
0635       double _r16 = _d_u * -(1. / (v * v));
0636       _d_v += _r16;
0637    } else {
0638       double _d_u = 0;
0639       double _t25 = ::std::log(v);
0640       double _t24 = (v + 1);
0641       double _t23 = (v - v * _t25 / _t24);
0642       double u = 1. / _t23;
0643       double _t26;
0644       _d_u += a2[3] * -d_out * u * u;
0645       _d_u += (a2[2] + a2[3] * u) * -d_out * u;
0646       _d_u += (a2[1] + (a2[2] + a2[3] * u) * u) * -d_out;
0647       double _r18 = _d_u * -(1. / (_t23 * _t23));
0648       _d_v += _r18;
0649       _d_v += -_r18 / _t24 * _t25;
0650       double _r19 = 0;
0651       _r19 += v * -_r18 / _t24 * clad::custom_derivatives::log_pushforward(v, 1.).pushforward;
0652       _d_v += _r19;
0653       double _r20 = -_r18 * -(v * _t25 / (_t24 * _t24));
0654       _d_v += _r20;
0655    }
0656 
0657    *d_x += _d_v / xi;
0658    *d_x0 += -_d_v / xi;
0659    *d_xi += _d_v * -((x - x0) / (xi * xi));
0660 }
0661 
0662 } // namespace Math
0663 } // namespace ROOT
0664 
0665 } // namespace custom_derivatives
0666 } // namespace clad
0667 
0668 #endif // CLAD_DERIVATOR
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!