Add analytical derivative support for std::beta

Author: leetcodez

include/clad/Differentiator/BuiltinDerivatives.h

@@ -1062,6 +1062,48 @@ CUDA_HOST_DEVICE void hypot_pullback(T x, T y, U d_z, T* d_x, T* d_y) {
   *d_y += (y / h) * d_z;
 }
 
+// 7. Special Functions
+#if __cplusplus >= 201703L
+template <typename T>
+CUDA_HOST_DEVICE inline T clad_digamma(T x) {
+  if (x <= 0.0) {
+    if (x == ::std::floor(x)) return (T)NAN; 
+    return clad_digamma(1.0 - x) - ::std::acos((T)-1.0) / ::std::tan(::std::acos((T)-1.0) * x);
+  }
+  T result = 0.0;
+  while (x < 8.0) {
+    result -= 1.0 / x;
+    x += 1.0;
+  }
+  T inv_x = 1.0 / x;
+  T inv_x2 = inv_x * inv_x;
+  result += ::std::log(x) - 0.5 * inv_x 
+            - inv_x2 * (1.0/12.0 
+            - inv_x2 * (1.0/120.0 
+            - inv_x2 * (1.0/252.0 
+            - inv_x2 * (1.0/240.0))));
+  return result;
+}
+
+template <typename T, typename dT>
+CUDA_HOST_DEVICE ValueAndPushforward<T, dT> beta_pushforward(T x, T y, dT d_x, dT d_y) {
+  T b = ::std::beta(x, y);
+  T psi_xy = clad_digamma(x + y);
+  dT pushforward = 0;
+  if (d_x) pushforward += b * (clad_digamma(x) - psi_xy) * d_x;
+  if (d_y) pushforward += b * (clad_digamma(y) - psi_xy) * d_y;
+  return {b, pushforward};
+}
+
+template <typename T, typename U>
+CUDA_HOST_DEVICE void beta_pullback(T x, T y, U d_z, T* d_x, T* d_y) {
+  T b = ::std::beta(x, y);
+  T psi_xy = clad_digamma(x + y);
+  if (d_x) *d_x += b * (clad_digamma(x) - psi_xy) * d_z;
+  if (d_y) *d_y += b * (clad_digamma(y) - psi_xy) * d_z;
+}
+#endif
+
 } // namespace std
 
 CUDA_HOST_DEVICE inline ValueAndPushforward<float, float>
@@ -1327,6 +1369,12 @@ using std::pow_pullback;
 using std::pow_pushforward;
 using std::sqrt_pushforward;
 
+// 7. Special Functions
+#if __cplusplus >= 201703L
+using std::beta_pullback;
+using std::beta_pushforward;
+#endif
+
 namespace class_functions {
 template <typename T, typename U>
 void constructor_pullback(ValueAndPushforward<T, U> rhs,

test/FirstDerivative/BuiltinDerivatives.C

@@ -515,6 +515,14 @@ double f_custom_min(double x, double y) { return std::min(x, y, std::greater<dou
 // CHECK-NEXT:      return _t0.pushforward;
 // CHECK-NEXT:  }
 
+double f_beta(double x, double y) { return std::beta(x, y); }
+// CHECK: double f_beta_darg0(double x, double y) {
+// CHECK-NEXT:     double _d_x = 1;
+// CHECK-NEXT:     double _d_y = 0;
+// CHECK-NEXT:     {{.*}}ValueAndPushforward<double, double> _t0 = {{.*}}beta_pushforward(x, y, _d_x, _d_y);
+// CHECK-NEXT:     return _t0.pushforward;
+// CHECK-NEXT: }
+
 int main () { //expected-no-diagnostics
   float f_result[2];
   double d_result[2];
@@ -695,5 +703,8 @@ int main () { //expected-no-diagnostics
   auto d_custom_min =  clad::differentiate(f_custom_min, 0);
   printf("Result is = %.6f\n", d_custom_min.execute(2, 3)); // CHECK-EXEC: Result is = 0.000000
 
+  auto d_beta = clad::differentiate(f_beta, 0);
+  printf("Result is = %.6f\n", d_beta.execute(2.0, 3.0)); // CHECK-EXEC: Result is = -0.090278
+
   return 0;
 }