Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 10:29:31 

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 
0028 // For the digamma function, that is the derivative of lgamma. We get it via
0029 // mathmore from the GSL, so the pullbacks that use digamma are only available
0030 // with mathmore=ON.
0031 #ifdef R__HAS_MATHMORE
0032 #include "Math/SpecFuncMathMore.h"
0033 #endif
0034 
0035 #include <stdexcept>
0036 
0037 namespace clad {
0038 namespace custom_derivatives {
0039 namespace TMath {
0040 template <typename T>
0041 ValueAndPushforward<T, T> Abs_pushforward(T x, T d_x)
0042 {
0043    return {::TMath::Abs(x), ((x < 0) ? -1 : 1) * d_x};
0044 }
0045 
0046 template <typename T>
0047 ValueAndPushforward<T, T> ACos_pushforward(T x, T d_x)
0048 {
0049    return {::TMath::ACos(x), (-1. / ::TMath::Sqrt(1 - x * x)) * d_x};
0050 }
0051 
0052 template <typename T>
0053 ValueAndPushforward<T, T> ACosH_pushforward(T x, T d_x)
0054 {
0055    return {::TMath::ACosH(x), (1. / ::TMath::Sqrt(x * x - 1)) * d_x};
0056 }
0057 
0058 template <typename T>
0059 ValueAndPushforward<T, T> ASin_pushforward(T x, T d_x)
0060 {
0061    return {::TMath::ASin(x), (1. / ::TMath::Sqrt(1 - x * x)) * d_x};
0062 }
0063 
0064 template <typename T>
0065 ValueAndPushforward<T, T> ASinH_pushforward(T x, T d_x)
0066 {
0067    return {::TMath::ASinH(x), (1. / ::TMath::Sqrt(x * x + 1)) * d_x};
0068 }
0069 
0070 template <typename T>
0071 ValueAndPushforward<T, T> ATan_pushforward(T x, T d_x)
0072 {
0073    return {::TMath::ATan(x), (1. / (x * x + 1)) * d_x};
0074 }
0075 
0076 template <typename T>
0077 ValueAndPushforward<T, T> ATanH_pushforward(T x, T d_x)
0078 {
0079    return {::TMath::ATanH(x), (1. / (1 - x * x)) * d_x};
0080 }
0081 
0082 template <typename T>
0083 ValueAndPushforward<T, T> Cos_pushforward(T x, T d_x)
0084 {
0085    return {::TMath::Cos(x), -::TMath::Sin(x) * d_x};
0086 }
0087 
0088 template <typename T>
0089 ValueAndPushforward<T, T> CosH_pushforward(T x, T d_x)
0090 {
0091    return {::TMath::CosH(x), ::TMath::SinH(x) * d_x};
0092 }
0093 
0094 template <typename T>
0095 ValueAndPushforward<T, T> Erf_pushforward(T x, T d_x)
0096 {
0097    return {::TMath::Erf(x), (2 * ::TMath::Exp(-x * x) / ::TMath::Sqrt(::TMath::Pi())) * d_x};
0098 }
0099 
0100 template <typename T>
0101 ValueAndPushforward<T, T> Erfc_pushforward(T x, T d_x)
0102 {
0103    return {::TMath::Erfc(x), -Erf_pushforward(x, d_x).pushforward};
0104 }
0105 
0106 #ifdef R__HAS_MATHMORE
0107 
0108 template <typename T>
0109 ValueAndPushforward<T, T> LnGamma_pushforward(T z, T d_z)
0110 {
0111    return {::TMath::LnGamma(z), ::ROOT::Math::digamma(z) * d_z};
0112 }
0113 
0114 #endif
0115 
0116 template <typename T>
0117 ValueAndPushforward<T, T> Exp_pushforward(T x, T d_x)
0118 {
0119    return {::TMath::Exp(x), ::TMath::Exp(x) * d_x};
0120 }
0121 
0122 template <typename T>
0123 ValueAndPushforward<T, T> Hypot_pushforward(T x, T y, T d_x, T d_y)
0124 {
0125    return {::TMath::Hypot(x, y), x / ::TMath::Hypot(x, y) * d_x + y / ::TMath::Hypot(x, y) * d_y};
0126 }
0127 
0128 template <typename T, typename U>
0129 void Hypot_pullback(T x, T y, U p, clad::array_ref<T> d_x, clad::array_ref<T> d_y)
0130 {
0131    T h = ::TMath::Hypot(x, y);
0132    *d_x += x / h * p;
0133    *d_y += y / h * p;
0134 }
0135 
0136 template <typename T>
0137 ValueAndPushforward<T, T> Log_pushforward(T x, T d_x)
0138 {
0139    return {::TMath::Log(x), (1. / x) * d_x};
0140 }
0141 
0142 template <typename T>
0143 ValueAndPushforward<T, T> Log10_pushforward(T x, T d_x)
0144 {
0145    return {::TMath::Log10(x), (1.0 / (x * ::TMath::Ln10())) * d_x};
0146 }
0147 
0148 template <typename T>
0149 ValueAndPushforward<T, T> Log2_pushforward(T x, T d_x)
0150 {
0151    return {::TMath::Log2(x), (1.0 / (x * ::TMath::Log(2.0))) * d_x};
0152 }
0153 
0154 template <typename T, typename U>
0155 ValueAndPushforward<T, T> Power_pushforward(T x, U y, T d_x, U d_y)
0156 {
0157    T pushforward = y * ::TMath::Power(x, y - 1) * d_x;
0158    if (d_y) {
0159       pushforward += (::TMath::Power(x, y) * ::TMath::Log(x)) * d_y;
0160    }
0161    return {::TMath::Power(x, y), pushforward};
0162 }
0163 
0164 template <typename T, typename U, typename V>
0165 void Power_pullback(T x, U y, V p, clad::array_ref<T> d_x, clad::array_ref<U> d_y)
0166 {
0167    auto t = pow_pushforward(x, y, 1, 0);
0168    *d_x += t.pushforward * p;
0169    t = pow_pushforward(x, y, 0, 1);
0170    *d_y += t.pushforward * p;
0171 }
0172 
0173 template <typename T>
0174 ValueAndPushforward<T, T> Sin_pushforward(T x, T d_x)
0175 {
0176    return {::TMath::Sin(x), ::TMath::Cos(x) * d_x};
0177 }
0178 
0179 template <typename T>
0180 ValueAndPushforward<T, T> SinH_pushforward(T x, T d_x)
0181 {
0182    return {::TMath::SinH(x), ::TMath::CosH(x) * d_x};
0183 }
0184 
0185 template <typename T>
0186 ValueAndPushforward<T, T> Sq_pushforward(T x, T d_x)
0187 {
0188    return {::TMath::Sq(x), 2 * x * d_x};
0189 }
0190 
0191 template <typename T>
0192 ValueAndPushforward<T, T> Sqrt_pushforward(T x, T d_x)
0193 {
0194    return {::TMath::Sqrt(x), (0.5 / ::TMath::Sqrt(x)) * d_x};
0195 }
0196 
0197 template <typename T>
0198 ValueAndPushforward<T, T> Tan_pushforward(T x, T d_x)
0199 {
0200    return {::TMath::Tan(x), (1. / ::TMath::Sq(::TMath::Cos(x))) * d_x};
0201 }
0202 
0203 template <typename T>
0204 ValueAndPushforward<T, T> TanH_pushforward(T x, T d_x)
0205 {
0206    return {::TMath::TanH(x), (1. / ::TMath::Sq(::TMath::CosH(x))) * d_x};
0207 }
0208 
0209 #ifdef WIN32
0210 // Additional custom derivatives that can be removed
0211 // after Issue #12108 in ROOT is resolved
0212 // constexpr is removed
0213 ValueAndPushforward<Double_t, Double_t> Pi_pushforward()
0214 {
0215    return {3.1415926535897931, 0.};
0216 }
0217 // constexpr is removed
0218 ValueAndPushforward<Double_t, Double_t> Ln10_pushforward()
0219 {
0220    return {2.3025850929940459, 0.};
0221 }
0222 #endif
0223 } // namespace TMath
0224 
0225 namespace ROOT {
0226 namespace Math {
0227 
0228 inline void landau_pdf_pullback(double x, double xi, double x0, double d_out, double *d_x, double *d_xi, double *d_x0)
0229 {
0230    if (xi <= 0) {
0231       return;
0232    }
0233    // clang-format off
0234    static double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635,   0.001511162253};
0235    static double q1[5] = {1.0         ,-0.3388260629, 0.09594393323, -0.01608042283,    0.003778942063};
0236 
0237    static double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411,   0.0001283617211};
0238    static double q2[5] = {1.0         , 0.7428795082, 0.3153932961,   0.06694219548,    0.008790609714};
0239 
0240    static double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
0241    static double q3[5] = {1.0         , 0.6097809921, 0.2560616665,   0.04746722384,    0.006957301675};
0242 
0243    static double p4[5] = {0.9874054407, 118.6723273,  849.2794360,   -743.7792444,      427.0262186};
0244    static double q4[5] = {1.0         , 106.8615961,  337.6496214,    2016.712389,      1597.063511};
0245 
0246    static double p5[5] = {1.003675074,  167.5702434,  4789.711289,    21217.86767,     -22324.94910};
0247    static double q5[5] = {1.0         , 156.9424537,  3745.310488,    9834.698876,      66924.28357};
0248 
0249    static double p6[5] = {1.000827619,  664.9143136,  62972.92665,    475554.6998,     -5743609.109};
0250    static double q6[5] = {1.0         , 651.4101098,  56974.73333,    165917.4725,     -2815759.939};
0251 
0252    static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
0253 
0254    static double a2[2] = {-1.845568670,-4.284640743};
0255    // clang-format on
0256    const double _const0 = 0.3989422803;
0257    double v = (x - x0) / xi;
0258    double _d_v = 0;
0259    double _d_denlan = 0;
0260    if (v < -5.5) {
0261       double u = ::std::exp(v + 1.);
0262       double _d_u = 0;
0263       if (u >= 1.e-10) {
0264          const double ue = ::std::exp(-1 / u);
0265          const double us = ::std::sqrt(u);
0266          double _t3;
0267          double _d_ue = 0;
0268          double _d_us = 0;
0269          double denlan = _const0 * (ue / us) * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
0270          _d_denlan += d_out / xi;
0271          *d_xi += d_out * -(denlan / (xi * xi));
0272          denlan = _t3;
0273          double _r_d3 = _d_denlan;
0274          _d_denlan -= _r_d3;
0275          _d_ue += _const0 * _r_d3 * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u) / us;
0276          double _r5 = _const0 * _r_d3 * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u) * -(ue / (us * us));
0277          _d_us += _r5;
0278          _d_u += a1[2] * _const0 * (ue / us) * _r_d3 * u * u;
0279          _d_u += (a1[1] + a1[2] * u) * _const0 * (ue / us) * _r_d3 * u;
0280          _d_u += (a1[0] + (a1[1] + a1[2] * u) * u) * _const0 * (ue / us) * _r_d3;
0281          double _r_d2 = _d_us;
0282          _d_us -= _r_d2;
0283          double _r4 = 0;
0284          _r4 += _r_d2 * clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0285          _d_u += _r4;
0286          double _r_d1 = _d_ue;
0287          _d_ue -= _r_d1;
0288          double _r2 = 0;
0289          _r2 += _r_d1 * ::std::exp(-1 / u);
0290          double _r3 = _r2 * -(-1 / (u * u));
0291          _d_u += _r3;
0292       }
0293       double _r_d0 = _d_u;
0294       _d_u -= _r_d0;
0295       double _r1 = 0;
0296       _r1 += _r_d0 * ::std::exp(v + 1.);
0297       _d_v += _r1;
0298    } else if (v < -1) {
0299       double _t4;
0300       double u = ::std::exp(-v - 1);
0301       double _d_u = 0;
0302       double _t5;
0303       double _t8 = ::std::exp(-u);
0304       double _t7 = ::std::sqrt(u);
0305       double _t6 = (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
0306       double denlan = _t8 * _t7 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / _t6;
0307       _d_denlan += d_out / xi;
0308       *d_xi += d_out * -(denlan / (xi * xi));
0309       denlan = _t5;
0310       double _r_d5 = _d_denlan;
0311       _d_denlan -= _r_d5;
0312       double _r7 = 0;
0313       _r7 += _r_d5 / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) * _t7 * ::std::exp(-u);
0314       _d_u += -_r7;
0315       double _r8 = 0;
0316       _r8 += _t8 * _r_d5 / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) *
0317              clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0318       _d_u += _r8;
0319       _d_v += p1[4] * _t8 * _t7 * _r_d5 / _t6 * v * v * v;
0320       _d_v += (p1[3] + p1[4] * v) * _t8 * _t7 * _r_d5 / _t6 * v * v;
0321       _d_v += (p1[2] + (p1[3] + p1[4] * v) * v) * _t8 * _t7 * _r_d5 / _t6 * v;
0322       _d_v += (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * _t8 * _t7 * _r_d5 / _t6;
0323       double _r9 = _r_d5 * -(_t8 * _t7 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (_t6 * _t6));
0324       _d_v += q1[4] * _r9 * v * v * v;
0325       _d_v += (q1[3] + q1[4] * v) * _r9 * v * v;
0326       _d_v += (q1[2] + (q1[3] + q1[4] * v) * v) * _r9 * v;
0327       _d_v += (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * _r9;
0328       u = _t4;
0329       double _r_d4 = _d_u;
0330       _d_u -= _r_d4;
0331       double _r6 = 0;
0332       _r6 += _r_d4 * ::std::exp(-v - 1);
0333       _d_v += -_r6;
0334    } else if (v < 1) {
0335       double _t9;
0336       double _t10 = (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
0337       double denlan = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / _t10;
0338       _d_denlan += d_out / xi;
0339       *d_xi += d_out * -(denlan / (xi * xi));
0340       denlan = _t9;
0341       double _r_d6 = _d_denlan;
0342       _d_denlan -= _r_d6;
0343       _d_v += p2[4] * _r_d6 / _t10 * v * v * v;
0344       _d_v += (p2[3] + p2[4] * v) * _r_d6 / _t10 * v * v;
0345       _d_v += (p2[2] + (p2[3] + p2[4] * v) * v) * _r_d6 / _t10 * v;
0346       _d_v += (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * _r_d6 / _t10;
0347       double _r10 = _r_d6 * -((p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / (_t10 * _t10));
0348       _d_v += q2[4] * _r10 * v * v * v;
0349       _d_v += (q2[3] + q2[4] * v) * _r10 * v * v;
0350       _d_v += (q2[2] + (q2[3] + q2[4] * v) * v) * _r10 * v;
0351       _d_v += (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * _r10;
0352    } else if (v < 5) {
0353       double _t11;
0354       double _t12 = (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
0355       double denlan = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / _t12;
0356       _d_denlan += d_out / xi;
0357       *d_xi += d_out * -(denlan / (xi * xi));
0358       denlan = _t11;
0359       double _r_d7 = _d_denlan;
0360       _d_denlan -= _r_d7;
0361       _d_v += p3[4] * _r_d7 / _t12 * v * v * v;
0362       _d_v += (p3[3] + p3[4] * v) * _r_d7 / _t12 * v * v;
0363       _d_v += (p3[2] + (p3[3] + p3[4] * v) * v) * _r_d7 / _t12 * v;
0364       _d_v += (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * _r_d7 / _t12;
0365       double _r11 = _r_d7 * -((p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / (_t12 * _t12));
0366       _d_v += q3[4] * _r11 * v * v * v;
0367       _d_v += (q3[3] + q3[4] * v) * _r11 * v * v;
0368       _d_v += (q3[2] + (q3[3] + q3[4] * v) * v) * _r11 * v;
0369       _d_v += (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * _r11;
0370    } else if (v < 12) {
0371       double u = 1 / v;
0372       double _d_u = 0;
0373       double _t14;
0374       double _t15 = (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
0375       double denlan = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / _t15;
0376       _d_denlan += d_out / xi;
0377       *d_xi += d_out * -(denlan / (xi * xi));
0378       denlan = _t14;
0379       double _r_d9 = _d_denlan;
0380       _d_denlan -= _r_d9;
0381       _d_u += _r_d9 / _t15 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) * u;
0382       _d_u += u * _r_d9 / _t15 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u);
0383       _d_u += p4[4] * u * u * _r_d9 / _t15 * u * u * u;
0384       _d_u += (p4[3] + p4[4] * u) * u * u * _r_d9 / _t15 * u * u;
0385       _d_u += (p4[2] + (p4[3] + p4[4] * u) * u) * u * u * _r_d9 / _t15 * u;
0386       _d_u += (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u * u * _r_d9 / _t15;
0387       double _r13 = _r_d9 * -(u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / (_t15 * _t15));
0388       _d_u += q4[4] * _r13 * u * u * u;
0389       _d_u += (q4[3] + q4[4] * u) * _r13 * u * u;
0390       _d_u += (q4[2] + (q4[3] + q4[4] * u) * u) * _r13 * u;
0391       _d_u += (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * _r13;
0392       double _r_d8 = _d_u;
0393       _d_u -= _r_d8;
0394       double _r12 = _r_d8 * -(1 / (v * v));
0395       _d_v += _r12;
0396    } else if (v < 50) {
0397       double u = 1 / v;
0398       double _d_u = 0;
0399       double _t17;
0400       double _t18 = (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
0401       double denlan = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / _t18;
0402       _d_denlan += d_out / xi;
0403       *d_xi += d_out * -(denlan / (xi * xi));
0404       denlan = _t17;
0405       double _r_d11 = _d_denlan;
0406       _d_denlan -= _r_d11;
0407       _d_u += _r_d11 / _t18 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) * u;
0408       _d_u += u * _r_d11 / _t18 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u);
0409       _d_u += p5[4] * u * u * _r_d11 / _t18 * u * u * u;
0410       _d_u += (p5[3] + p5[4] * u) * u * u * _r_d11 / _t18 * u * u;
0411       _d_u += (p5[2] + (p5[3] + p5[4] * u) * u) * u * u * _r_d11 / _t18 * u;
0412       _d_u += (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u * u * _r_d11 / _t18;
0413       double _r15 = _r_d11 * -(u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / (_t18 * _t18));
0414       _d_u += q5[4] * _r15 * u * u * u;
0415       _d_u += (q5[3] + q5[4] * u) * _r15 * u * u;
0416       _d_u += (q5[2] + (q5[3] + q5[4] * u) * u) * _r15 * u;
0417       _d_u += (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * _r15;
0418       double _r_d10 = _d_u;
0419       _d_u -= _r_d10;
0420       double _r14 = _r_d10 * -(1 / (v * v));
0421       _d_v += _r14;
0422    } else if (v < 300) {
0423       double _t19;
0424       double u = 1 / v;
0425       double _d_u = 0;
0426       double _t20;
0427       double _t21 = (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
0428       double denlan = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / _t21;
0429       _d_denlan += d_out / xi;
0430       *d_xi += d_out * -(denlan / (xi * xi));
0431       denlan = _t20;
0432       double _r_d13 = _d_denlan;
0433       _d_denlan -= _r_d13;
0434       _d_u += _r_d13 / _t21 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) * u;
0435       _d_u += u * _r_d13 / _t21 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u);
0436       _d_u += p6[4] * u * u * _r_d13 / _t21 * u * u * u;
0437       _d_u += (p6[3] + p6[4] * u) * u * u * _r_d13 / _t21 * u * u;
0438       _d_u += (p6[2] + (p6[3] + p6[4] * u) * u) * u * u * _r_d13 / _t21 * u;
0439       _d_u += (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u * u * _r_d13 / _t21;
0440       double _r17 = _r_d13 * -(u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / (_t21 * _t21));
0441       _d_u += q6[4] * _r17 * u * u * u;
0442       _d_u += (q6[3] + q6[4] * u) * _r17 * u * u;
0443       _d_u += (q6[2] + (q6[3] + q6[4] * u) * u) * _r17 * u;
0444       _d_u += (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * _r17;
0445       u = _t19;
0446       double _r_d12 = _d_u;
0447       _d_u -= _r_d12;
0448       double _r16 = _r_d12 * -(1 / (v * v));
0449       _d_v += _r16;
0450    } else {
0451       double _t22;
0452       double _t25 = ::std::log(v);
0453       double _t24 = (v + 1);
0454       double _t23 = (v - v * _t25 / _t24);
0455       double u = 1 / _t23;
0456       double _d_u = 0;
0457       double _t26;
0458       double denlan = u * u * (1 + (a2[0] + a2[1] * u) * u);
0459       _d_denlan += d_out / xi;
0460       *d_xi += d_out * -(denlan / (xi * xi));
0461       denlan = _t26;
0462       double _r_d15 = _d_denlan;
0463       _d_denlan -= _r_d15;
0464       _d_u += _r_d15 * (1 + (a2[0] + a2[1] * u) * u) * u;
0465       _d_u += u * _r_d15 * (1 + (a2[0] + a2[1] * u) * u);
0466       _d_u += a2[1] * u * u * _r_d15 * u;
0467       _d_u += (a2[0] + a2[1] * u) * u * u * _r_d15;
0468       u = _t22;
0469       double _r_d14 = _d_u;
0470       _d_u -= _r_d14;
0471       double _r18 = _r_d14 * -(1 / (_t23 * _t23));
0472       _d_v += _r18;
0473       _d_v += -_r18 / _t24 * _t25;
0474       double _r19 = 0;
0475       _r19 += v * -_r18 / _t24 / v;
0476       _d_v += _r19;
0477       double _r20 = -_r18 * -(v * _t25 / (_t24 * _t24));
0478       _d_v += _r20;
0479    }
0480    *d_x += _d_v / xi;
0481    *d_x0 += -_d_v / xi;
0482    double _r0 = _d_v * -((x - x0) / (xi * xi));
0483    *d_xi += _r0;
0484 }
0485 
0486 inline void landau_cdf_pullback(double x, double xi, double x0, double d_out, double *d_x, double *d_xi, double *d_x0)
0487 {
0488    // clang-format off
0489    static double p1[5] = {0.2514091491e+0,-0.6250580444e-1, 0.1458381230e-1,-0.2108817737e-2, 0.7411247290e-3};
0490    static double q1[5] = {1.0            ,-0.5571175625e-2, 0.6225310236e-1,-0.3137378427e-2, 0.1931496439e-2};
0491 
0492    static double p2[4] = {0.2868328584e+0, 0.3564363231e+0, 0.1523518695e+0, 0.2251304883e-1};
0493    static double q2[4] = {1.0            , 0.6191136137e+0, 0.1720721448e+0, 0.2278594771e-1};
0494 
0495    static double p3[4] = {0.2868329066e+0, 0.3003828436e+0, 0.9950951941e-1, 0.8733827185e-2};
0496    static double q3[4] = {1.0            , 0.4237190502e+0, 0.1095631512e+0, 0.8693851567e-2};
0497 
0498    static double p4[4] = {0.1000351630e+1, 0.4503592498e+1, 0.1085883880e+2, 0.7536052269e+1};
0499    static double q4[4] = {1.0            , 0.5539969678e+1, 0.1933581111e+2, 0.2721321508e+2};
0500 
0501    static double p5[4] = {0.1000006517e+1, 0.4909414111e+2, 0.8505544753e+2, 0.1532153455e+3};
0502    static double q5[4] = {1.0            , 0.5009928881e+2, 0.1399819104e+3, 0.4200002909e+3};
0503 
0504    static double p6[4] = {0.1000000983e+1, 0.1329868456e+3, 0.9162149244e+3,-0.9605054274e+3};
0505    static double q6[4] = {1.0            , 0.1339887843e+3, 0.1055990413e+4, 0.5532224619e+3};
0506 
0507    static double a1[4] = {0              ,-0.4583333333e+0, 0.6675347222e+0,-0.1641741416e+1};
0508    static double a2[4] = {0              , 1.0            ,-0.4227843351e+0,-0.2043403138e+1};
0509    // clang-format on
0510 
0511    const double v = (x - x0) / xi;
0512    double _d_v = 0;
0513    if (v < -5.5) {
0514       double _d_u = 0;
0515       const double _const0 = 0.3989422803;
0516       double u = ::std::exp(v + 1);
0517       double _t3 = ::std::exp(-1. / u);
0518       double _t2 = ::std::sqrt(u);
0519       double _r2 = 0;
0520       _r2 += _const0 * d_out * (1 + (a1[1] + (a1[2] + a1[3] * u) * u) * u) * _t2 * ::std::exp(-1. / u);
0521       double _r3 = _r2 * -(-1. / (u * u));
0522       _d_u += _r3;
0523       double _r4 = 0;
0524       _r4 += _const0 * _t3 * d_out * (1 + (a1[1] + (a1[2] + a1[3] * u) * u) * u) *
0525              clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0526       _d_u += _r4;
0527       _d_u += a1[3] * _const0 * _t3 * _t2 * d_out * u * u;
0528       _d_u += (a1[2] + a1[3] * u) * _const0 * _t3 * _t2 * d_out * u;
0529       _d_u += (a1[1] + (a1[2] + a1[3] * u) * u) * _const0 * _t3 * _t2 * d_out;
0530       _d_v += _d_u * ::std::exp(v + 1);
0531    } else if (v < -1) {
0532       double _d_u = 0;
0533       double u = ::std::exp(-v - 1);
0534       double _t8 = ::std::exp(-u);
0535       double _t7 = ::std::sqrt(u);
0536       double _t6 = (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
0537       double _r6 = 0;
0538       _r6 += d_out / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / _t7 * ::std::exp(-u);
0539       _d_u += -_r6;
0540       double _r7 = d_out / _t6 * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) * -(_t8 / (_t7 * _t7));
0541       double _r8 = 0;
0542       _r8 += _r7 * clad::custom_derivatives::sqrt_pushforward(u, 1.).pushforward;
0543       _d_u += _r8;
0544       _d_v += p1[4] * (_t8 / _t7) * d_out / _t6 * v * v * v;
0545       _d_v += (p1[3] + p1[4] * v) * (_t8 / _t7) * d_out / _t6 * v * v;
0546       _d_v += (p1[2] + (p1[3] + p1[4] * v) * v) * (_t8 / _t7) * d_out / _t6 * v;
0547       _d_v += (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * (_t8 / _t7) * d_out / _t6;
0548       double _r9 = d_out * -((_t8 / _t7) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (_t6 * _t6));
0549       _d_v += q1[4] * _r9 * v * v * v;
0550       _d_v += (q1[3] + q1[4] * v) * _r9 * v * v;
0551       _d_v += (q1[2] + (q1[3] + q1[4] * v) * v) * _r9 * v;
0552       _d_v += (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * _r9;
0553       _d_v += -_d_u * ::std::exp(-v - 1);
0554    } else if (v < 1) {
0555       double _t10 = (q2[0] + (q2[1] + (q2[2] + q2[3] * v) * v) * v);
0556       _d_v += p2[3] * d_out / _t10 * v * v;
0557       _d_v += (p2[2] + p2[3] * v) * d_out / _t10 * v;
0558       _d_v += (p2[1] + (p2[2] + p2[3] * v) * v) * d_out / _t10;
0559       double _r10 = d_out * -((p2[0] + (p2[1] + (p2[2] + p2[3] * v) * v) * v) / (_t10 * _t10));
0560       _d_v += q2[3] * _r10 * v * v;
0561       _d_v += (q2[2] + q2[3] * v) * _r10 * v;
0562       _d_v += (q2[1] + (q2[2] + q2[3] * v) * v) * _r10;
0563    } else if (v < 4) {
0564       double _t12 = (q3[0] + (q3[1] + (q3[2] + q3[3] * v) * v) * v);
0565       _d_v += p3[3] * d_out / _t12 * v * v;
0566       _d_v += (p3[2] + p3[3] * v) * d_out / _t12 * v;
0567       _d_v += (p3[1] + (p3[2] + p3[3] * v) * v) * d_out / _t12;
0568       double _r11 = d_out * -((p3[0] + (p3[1] + (p3[2] + p3[3] * v) * v) * v) / (_t12 * _t12));
0569       _d_v += q3[3] * _r11 * v * v;
0570       _d_v += (q3[2] + q3[3] * v) * _r11 * v;
0571       _d_v += (q3[1] + (q3[2] + q3[3] * v) * v) * _r11;
0572    } else if (v < 12) {
0573       double _d_u = 0;
0574       double u = 1. / v;
0575       double _t15 = (q4[0] + (q4[1] + (q4[2] + q4[3] * u) * u) * u);
0576       _d_u += p4[3] * d_out / _t15 * u * u;
0577       _d_u += (p4[2] + p4[3] * u) * d_out / _t15 * u;
0578       _d_u += (p4[1] + (p4[2] + p4[3] * u) * u) * d_out / _t15;
0579       double _r13 = d_out * -((p4[0] + (p4[1] + (p4[2] + p4[3] * u) * u) * u) / (_t15 * _t15));
0580       _d_u += q4[3] * _r13 * u * u;
0581       _d_u += (q4[2] + q4[3] * u) * _r13 * u;
0582       _d_u += (q4[1] + (q4[2] + q4[3] * u) * u) * _r13;
0583       double _r12 = _d_u * -(1. / (v * v));
0584       _d_v += _r12;
0585    } else if (v < 50) {
0586       double _d_u = 0;
0587       double u = 1. / v;
0588       double _t18 = (q5[0] + (q5[1] + (q5[2] + q5[3] * u) * u) * u);
0589       _d_u += p5[3] * d_out / _t18 * u * u;
0590       _d_u += (p5[2] + p5[3] * u) * d_out / _t18 * u;
0591       _d_u += (p5[1] + (p5[2] + p5[3] * u) * u) * d_out / _t18;
0592       double _r15 = d_out * -((p5[0] + (p5[1] + (p5[2] + p5[3] * u) * u) * u) / (_t18 * _t18));
0593       _d_u += q5[3] * _r15 * u * u;
0594       _d_u += (q5[2] + q5[3] * u) * _r15 * u;
0595       _d_u += (q5[1] + (q5[2] + q5[3] * u) * u) * _r15;
0596       double _r14 = _d_u * -(1. / (v * v));
0597       _d_v += _r14;
0598    } else if (v < 300) {
0599       double _d_u = 0;
0600       double u = 1. / v;
0601       double _t21 = (q6[0] + (q6[1] + (q6[2] + q6[3] * u) * u) * u);
0602       _d_u += p6[3] * d_out / _t21 * u * u;
0603       _d_u += (p6[2] + p6[3] * u) * d_out / _t21 * u;
0604       _d_u += (p6[1] + (p6[2] + p6[3] * u) * u) * d_out / _t21;
0605       double _r17 = d_out * -((p6[0] + (p6[1] + (p6[2] + p6[3] * u) * u) * u) / (_t21 * _t21));
0606       _d_u += q6[3] * _r17 * u * u;
0607       _d_u += (q6[2] + q6[3] * u) * _r17 * u;
0608       _d_u += (q6[1] + (q6[2] + q6[3] * u) * u) * _r17;
0609       double _r16 = _d_u * -(1. / (v * v));
0610       _d_v += _r16;
0611    } else {
0612       double _d_u = 0;
0613       double _t25 = ::std::log(v);
0614       double _t24 = (v + 1);
0615       double _t23 = (v - v * _t25 / _t24);
0616       double u = 1. / _t23;
0617       double _t26;
0618       _d_u += a2[3] * -d_out * u * u;
0619       _d_u += (a2[2] + a2[3] * u) * -d_out * u;
0620       _d_u += (a2[1] + (a2[2] + a2[3] * u) * u) * -d_out;
0621       double _r18 = _d_u * -(1. / (_t23 * _t23));
0622       _d_v += _r18;
0623       _d_v += -_r18 / _t24 * _t25;
0624       double _r19 = 0;
0625       _r19 += v * -_r18 / _t24 / v;
0626       _d_v += _r19;
0627       double _r20 = -_r18 * -(v * _t25 / (_t24 * _t24));
0628       _d_v += _r20;
0629    }
0630 
0631    *d_x += _d_v / xi;
0632    *d_x0 += -_d_v / xi;
0633    *d_xi += _d_v * -((x - x0) / (xi * xi));
0634 }
0635 
0636 #ifdef R__HAS_MATHMORE
0637 
0638 inline void inc_gamma_c_pullback(double a, double x, double _d_y, double *_d_a, double *_d_x);
0639 
0640 inline void inc_gamma_pullback(double a, double x, double _d_y, double *_d_a, double *_d_x)
0641 {
0642    // Synced with SpecFuncCephes.h
0643    constexpr double kMACHEP = 1.11022302462515654042363166809e-16;
0644    constexpr double kMAXLOG = 709.782712893383973096206318587;
0645    constexpr double kMINLOG = -708.396418532264078748994506896;
0646    constexpr double kMAXSTIR = 108.116855767857671821730036754;
0647    constexpr double kMAXLGM = 2.556348e305;
0648    constexpr double kBig = 4.503599627370496e15;
0649    constexpr double kBiginv = 2.22044604925031308085e-16;
0650 
0651    double _d_ans = 0, _d_ax = 0, _d_c = 0, _d_r = 0;
0652    double _t1;
0653    double _t2;
0654    double _t3;
0655    double _t4;
0656    double _t5;
0657    clad::tape<double> _t7 = {};
0658    clad::tape<double> _t8 = {};
0659    clad::tape<double> _t9 = {};
0660    double ans, ax, c, r;
0661    if (a <= 0)
0662       return;
0663    if (x <= 0)
0664       return;
0665    if ((x > 1.) && (x > a)) {
0666       double _r0 = 0;
0667       double _r1 = 0;
0668       inc_gamma_c_pullback(a, x, -_d_y, &_r0, &_r1);
0669       *_d_a += _r0;
0670       *_d_x += _r1;
0671       return;
0672    }
0673    _t1 = ::std::log(x);
0674    ax = a * _t1 - x - ::std::lgamma(a);
0675    if (ax < -kMAXLOG) {
0676       *_d_x += (a * _d_ax / x) - _d_ax;
0677       *_d_a +=
0678          _d_ax *
0679          (_t1 - ::ROOT::Math::digamma(a)); // numerical_diff::forward_central_difference(::std::lgamma, a, 0, 0, a);
0680       _d_ax = 0.;
0681       return;
0682    }
0683    _t2 = ax;
0684    ax = ::std::exp(ax);
0685    _t3 = r;
0686    r = a;
0687    _t4 = c;
0688    c = 1.;
0689    _t5 = ans;
0690    ans = 1.;
0691    unsigned long _t6 = 0;
0692    do {
0693       _t6++;
0694       clad::push(_t7, r);
0695       r += 1.;
0696       clad::push(_t8, c);
0697       c *= x / r;
0698       clad::push(_t9, ans);
0699       ans += c;
0700    } while (c / ans > kMACHEP);
0701    {
0702       _d_ans += _d_y / a * ax;
0703       _d_ax += ans * _d_y / a;
0704       double _r6 = _d_y * -(ans * ax / (a * a));
0705       *_d_a += _r6;
0706    }
0707    do {
0708       {
0709          {
0710             ans = clad::pop(_t9);
0711             double _r_d7 = _d_ans;
0712             _d_c += _r_d7;
0713          }
0714          {
0715             c = clad::pop(_t8);
0716             double _r_d6 = _d_c;
0717             _d_c -= _r_d6;
0718             _d_c += _r_d6 * x / r;
0719             *_d_x += c * _r_d6 / r;
0720             double _r5 = c * _r_d6 * -(x / (r * r));
0721             _d_r += _r5;
0722          }
0723          {
0724             r = clad::pop(_t7);
0725             double _r_d5 = _d_r;
0726          }
0727       }
0728       _t6--;
0729    } while (_t6);
0730    {
0731       ans = _t5;
0732       double _r_d4 = _d_ans;
0733       _d_ans -= _r_d4;
0734    }
0735    {
0736       c = _t4;
0737       double _r_d3 = _d_c;
0738       _d_c -= _r_d3;
0739    }
0740    {
0741       r = _t3;
0742       double _r_d2 = _d_r;
0743       _d_r -= _r_d2;
0744       *_d_a += _r_d2;
0745    }
0746    {
0747       ax = _t2;
0748       double _r_d1 = _d_ax;
0749       _d_ax -= _r_d1;
0750       double _r4 = 0;
0751       _r4 += _r_d1 * ::std::exp(ax);
0752       _d_ax += _r4;
0753    }
0754    {
0755       *_d_x += (a * _d_ax / x) - _d_ax;
0756       *_d_a +=
0757          _d_ax *
0758          (_t1 - ::ROOT::Math::digamma(a)); // numerical_diff::forward_central_difference(::std::lgamma, a, 0, 0, a);
0759       _d_ax = 0.;
0760    }
0761 }
0762 
0763 inline void inc_gamma_c_pullback(double a, double x, double _d_y, double *_d_a, double *_d_x)
0764 {
0765    // Synced with SpecFuncCephes.h
0766    constexpr double kMACHEP = 1.11022302462515654042363166809e-16;
0767    constexpr double kMAXLOG = 709.782712893383973096206318587;
0768    constexpr double kMINLOG = -708.396418532264078748994506896;
0769    constexpr double kMAXSTIR = 108.116855767857671821730036754;
0770    constexpr double kMAXLGM = 2.556348e305;
0771    constexpr double kBig = 4.503599627370496e15;
0772    constexpr double kBiginv = 2.22044604925031308085e-16;
0773 
0774    double _d_ans = 0, _d_ax = 0, _d_c = 0, _d_yc = 0, _d_r = 0, _d_t = 0, _d_y0 = 0, _d_z = 0;
0775    double _d_pk = 0, _d_pkm1 = 0, _d_pkm2 = 0, _d_qk = 0, _d_qkm1 = 0, _d_qkm2 = 0;
0776    double _t1;
0777    double _t2;
0778    double _t3;
0779    double _t4;
0780    double _t5;
0781    double _t6;
0782    double _t7;
0783    double _t8;
0784    double _t9;
0785    double _t10;
0786    unsigned long _t11;
0787    clad::tape<double> _t12 = {};
0788    clad::tape<double> _t13 = {};
0789    clad::tape<double> _t14 = {};
0790    clad::tape<double> _t15 = {};
0791    clad::tape<double> _t16 = {};
0792    clad::tape<double> _t17 = {};
0793    clad::tape<double> _t19 = {};
0794    clad::tape<double> _t20 = {};
0795    clad::tape<double> _t21 = {};
0796    clad::tape<double> _t22 = {};
0797    clad::tape<double> _t23 = {};
0798    clad::tape<double> _t24 = {};
0799    clad::tape<double> _t25 = {};
0800    clad::tape<double> _t26 = {};
0801    clad::tape<double> _t27 = {};
0802    clad::tape<bool> _t29 = {};
0803    clad::tape<double> _t30 = {};
0804    clad::tape<double> _t31 = {};
0805    clad::tape<double> _t32 = {};
0806    clad::tape<double> _t33 = {};
0807    double ans, ax, c, yc, r, t, y, z;
0808    double pk, pkm1, pkm2, qk, qkm1, qkm2;
0809    if (a <= 0)
0810       return;
0811    if (x <= 0)
0812       return;
0813    if ((x < 1.) || (x < a)) {
0814       double _r0 = 0;
0815       double _r1 = 0;
0816       inc_gamma_pullback(a, x, -_d_y, &_r0, &_r1);
0817       *_d_a += _r0;
0818       *_d_x += _r1;
0819       return;
0820    }
0821    _t1 = ::std::log(x);
0822    ax = a * _t1 - x - ::std::lgamma(a);
0823    if (ax < -kMAXLOG) {
0824       *_d_x += a * _d_ax / x - _d_ax;
0825       *_d_a +=
0826          _d_ax *
0827          (_t1 - ::ROOT::Math::digamma(a)); // numerical_diff::forward_central_difference(::std::lgamma, a, 0, 0, a);
0828       _d_ax = 0.;
0829       return;
0830    }
0831    _t2 = ax;
0832    ax = ::std::exp(ax);
0833    _t3 = y;
0834    y = 1. - a;
0835    _t4 = z;
0836    z = x + y + 1.;
0837    _t5 = c;
0838    c = 0.;
0839    _t6 = pkm2;
0840    pkm2 = 1.;
0841    _t7 = qkm2;
0842    qkm2 = x;
0843    _t8 = pkm1;
0844    pkm1 = x + 1.;
0845    _t9 = qkm1;
0846    qkm1 = z * x;
0847    _t10 = ans;
0848    ans = pkm1 / qkm1;
0849    _t11 = 0;
0850    do {
0851       _t11++;
0852       clad::push(_t12, c);
0853       c += 1.;
0854       clad::push(_t13, y);
0855       y += 1.;
0856       clad::push(_t14, z);
0857       z += 2.;
0858       clad::push(_t15, yc);
0859       yc = y * c;
0860       clad::push(_t16, pk);
0861       pk = pkm1 * z - pkm2 * yc;
0862       clad::push(_t17, qk);
0863       qk = qkm1 * z - qkm2 * yc;
0864       double _t18 = qk;
0865       {
0866          if (_t18) {
0867             clad::push(_t20, r);
0868             r = pk / qk;
0869             clad::push(_t21, t);
0870             t = ::std::abs((ans - r) / r);
0871             clad::push(_t22, ans);
0872             ans = r;
0873          } else {
0874             clad::push(_t23, t);
0875             t = 1.;
0876          }
0877          clad::push(_t19, _t18);
0878       }
0879       clad::push(_t24, pkm2);
0880       pkm2 = pkm1;
0881       clad::push(_t25, pkm1);
0882       pkm1 = pk;
0883       clad::push(_t26, qkm2);
0884       qkm2 = qkm1;
0885       clad::push(_t27, qkm1);
0886       qkm1 = qk;
0887       bool _t28 = ::std::abs(pk) > kBig;
0888       {
0889          if (_t28) {
0890             clad::push(_t30, pkm2);
0891             pkm2 *= kBiginv;
0892             clad::push(_t31, pkm1);
0893             pkm1 *= kBiginv;
0894             clad::push(_t32, qkm2);
0895             qkm2 *= kBiginv;
0896             clad::push(_t33, qkm1);
0897             qkm1 *= kBiginv;
0898          }
0899          clad::push(_t29, _t28);
0900       }
0901    } while (t > kMACHEP);
0902    {
0903       _d_ans += _d_y * ax;
0904       _d_ax += ans * _d_y;
0905    }
0906    do {
0907       {
0908          if (clad::pop(_t29)) {
0909             {
0910                qkm1 = clad::pop(_t33);
0911                double _r_d27 = _d_qkm1;
0912                _d_qkm1 -= _r_d27;
0913                _d_qkm1 += _r_d27 * kBiginv;
0914             }
0915             {
0916                qkm2 = clad::pop(_t32);
0917                double _r_d26 = _d_qkm2;
0918                _d_qkm2 -= _r_d26;
0919                _d_qkm2 += _r_d26 * kBiginv;
0920             }
0921             {
0922                pkm1 = clad::pop(_t31);
0923                double _r_d25 = _d_pkm1;
0924                _d_pkm1 -= _r_d25;
0925                _d_pkm1 += _r_d25 * kBiginv;
0926             }
0927             {
0928                pkm2 = clad::pop(_t30);
0929                double _r_d24 = _d_pkm2;
0930                _d_pkm2 -= _r_d24;
0931                _d_pkm2 += _r_d24 * kBiginv;
0932             }
0933          }
0934          {
0935             qkm1 = clad::pop(_t27);
0936             double _r_d23 = _d_qkm1;
0937             _d_qkm1 -= _r_d23;
0938             _d_qk += _r_d23;
0939          }
0940          {
0941             qkm2 = clad::pop(_t26);
0942             double _r_d22 = _d_qkm2;
0943             _d_qkm2 -= _r_d22;
0944             _d_qkm1 += _r_d22;
0945          }
0946          {
0947             pkm1 = clad::pop(_t25);
0948             double _r_d21 = _d_pkm1;
0949             _d_pkm1 -= _r_d21;
0950             _d_pk += _r_d21;
0951          }
0952          {
0953             pkm2 = clad::pop(_t24);
0954             double _r_d20 = _d_pkm2;
0955             _d_pkm2 -= _r_d20;
0956             _d_pkm1 += _r_d20;
0957          }
0958          if (clad::pop(_t19)) {
0959             {
0960                ans = clad::pop(_t22);
0961                double _r_d18 = _d_ans;
0962                _d_ans -= _r_d18;
0963                _d_r += _r_d18;
0964             }
0965             {
0966                t = clad::pop(_t21);
0967                double _r_d17 = _d_t;
0968                _d_t -= _r_d17;
0969                double _r7 = 0;
0970                _r7 += _r_d17 * clad::custom_derivatives::std::abs_pushforward((ans - r) / r, 1.).pushforward;
0971                _d_ans += _r7 / r;
0972                _d_r += -_r7 / r;
0973                double _r8 = _r7 * -((ans - r) / (r * r));
0974                _d_r += _r8;
0975             }
0976             {
0977                r = clad::pop(_t20);
0978                double _r_d16 = _d_r;
0979                _d_r -= _r_d16;
0980                _d_pk += _r_d16 / qk;
0981                double _r6 = _r_d16 * -(pk / (qk * qk));
0982                _d_qk += _r6;
0983             }
0984          } else {
0985             t = clad::pop(_t23);
0986             double _r_d19 = _d_t;
0987             _d_t -= _r_d19;
0988          }
0989          {
0990             qk = clad::pop(_t17);
0991             double _r_d15 = _d_qk;
0992             _d_qk -= _r_d15;
0993             _d_qkm1 += _r_d15 * z;
0994             _d_z += qkm1 * _r_d15;
0995             _d_qkm2 += -_r_d15 * yc;
0996             _d_yc += qkm2 * -_r_d15;
0997          }
0998          {
0999             pk = clad::pop(_t16);
1000             double _r_d14 = _d_pk;
1001             _d_pk -= _r_d14;
1002             _d_pkm1 += _r_d14 * z;
1003             _d_z += pkm1 * _r_d14;
1004             _d_pkm2 += -_r_d14 * yc;
1005             _d_yc += pkm2 * -_r_d14;
1006          }
1007          {
1008             yc = clad::pop(_t15);
1009             double _r_d13 = _d_yc;
1010             _d_yc -= _r_d13;
1011             _d_y0 += _r_d13 * c;
1012             _d_c += y * _r_d13;
1013          }
1014          {
1015             z = clad::pop(_t14);
1016             double _r_d12 = _d_z;
1017          }
1018          {
1019             y = clad::pop(_t13);
1020             double _r_d11 = _d_y0;
1021          }
1022          {
1023             c = clad::pop(_t12);
1024             double _r_d10 = _d_c;
1025          }
1026       }
1027       _t11--;
1028    } while (_t11);
1029    {
1030       ans = _t10;
1031       double _r_d9 = _d_ans;
1032       _d_ans -= _r_d9;
1033       _d_pkm1 += _r_d9 / qkm1;
1034       double _r5 = _r_d9 * -(pkm1 / (qkm1 * qkm1));
1035       _d_qkm1 += _r5;
1036    }
1037    {
1038       qkm1 = _t9;
1039       double _r_d8 = _d_qkm1;
1040       _d_qkm1 -= _r_d8;
1041       _d_z += _r_d8 * x;
1042       *_d_x += z * _r_d8;
1043    }
1044    {
1045       pkm1 = _t8;
1046       double _r_d7 = _d_pkm1;
1047       _d_pkm1 -= _r_d7;
1048       *_d_x += _r_d7;
1049    }
1050    {
1051       qkm2 = _t7;
1052       double _r_d6 = _d_qkm2;
1053       _d_qkm2 -= _r_d6;
1054       *_d_x += _r_d6;
1055    }
1056    {
1057       pkm2 = _t6;
1058       double _r_d5 = _d_pkm2;
1059       _d_pkm2 -= _r_d5;
1060    }
1061    {
1062       c = _t5;
1063       double _r_d4 = _d_c;
1064       _d_c -= _r_d4;
1065    }
1066    {
1067       z = _t4;
1068       double _r_d3 = _d_z;
1069       _d_z -= _r_d3;
1070       *_d_x += _r_d3;
1071       _d_y0 += _r_d3;
1072    }
1073    {
1074       y = _t3;
1075       double _r_d2 = _d_y0;
1076       _d_y0 -= _r_d2;
1077       *_d_a += -_r_d2;
1078    }
1079    {
1080       ax = _t2;
1081       double _r_d1 = _d_ax;
1082       _d_ax -= _r_d1;
1083       double _r4 = _r_d1 * ::std::exp(ax);
1084       _d_ax += _r4;
1085    }
1086    {
1087       *_d_x += a * _d_ax / x - _d_ax;
1088       *_d_a +=
1089          _d_ax *
1090          (_t1 - ::ROOT::Math::digamma(a)); // numerical_diff::forward_central_difference(::std::lgamma, a, 0, 0, a);
1091       _d_ax = 0.;
1092    }
1093 }
1094 
1095 #endif // R__HAS_MATHMORE
1096 
1097 } // namespace Math
1098 } // namespace ROOT
1099 
1100 } // namespace custom_derivatives
1101 } // namespace clad
1102 
1103 // Forward declare BLAS functions.
1104 extern "C" void sgemm_(const char *transa, const char *transb, const int *m, const int *n, const int *k,
1105                        const float *alpha, const float *A, const int *lda, const float *B, const int *ldb,
1106                        const float *beta, float *C, const int *ldc);
1107 
1108 namespace clad::custom_derivatives {
1109 
1110 namespace TMVA::Experimental::SOFIE {
1111 
1112 inline void Gemm_Call_pullback(float *output, bool transa, bool transb, int m, int n, int k, float alpha,
1113                                const float *A, const float *B, float beta, const float *C, float *_d_output, bool *,
1114                                bool *, int *, int *, int *, float *_d_alpha, float *_d_A, float *_d_B, float *_d_beta,
1115                                float *_d_C)
1116 {
1117    // TODO:
1118    //    - handle transa and transb cases correctly
1119    if (transa || transb) {
1120       return;
1121    }
1122 
1123    char ct = 't';
1124    char cn = 'n';
1125 
1126    // beta needs to be one because we want to add to _d_A and _d_B instead of
1127    // overwriting it.
1128    float one = 1.;
1129 
1130    // _d_A, _d_B
1131    // note: beta needs to be one because we want to add to _d_A and _d_B instead of overwriting it.
1132    ::sgemm_(&cn, &ct, &m, &k, &n, &alpha, _d_output, &m, B, &k, &one, _d_A, &m);
1133    ::sgemm_(&ct, &cn, &k, &n, &m, &alpha, A, &m, _d_output, &m, &one, _d_B, &k);
1134 
1135    // _d_alpha, _d_beta, _d_C
1136    int sizeC = n * m;
1137    for (int i = 0; i < sizeC; ++i) {
1138       *_d_alpha += _d_output[i] * (output[i] - beta * C[i]);
1139       *_d_beta += _d_output[i] * C[i];
1140       _d_C[i] += _d_output[i] * beta;
1141    }
1142 }
1143 
1144 } // namespace TMVA::Experimental::SOFIE
1145 
1146 } // namespace clad::custom_derivatives
1147 
1148 #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!