Expanded numerical stability and performance for std_normal_lcdf, reused across distributions & functions

stan/math/opencl/kernel_generator/elt_function_cl.hpp

@@ -7,6 +7,7 @@
 #include <stan/math/opencl/kernels/device_functions/beta.hpp>
 #include <stan/math/opencl/kernels/device_functions/digamma.hpp>
 #include <stan/math/opencl/kernels/device_functions/erfcx.hpp>
+#include <stan/math/opencl/kernels/device_functions/exp_mod_normal_lcdf.hpp>
 #include <stan/math/opencl/kernels/device_functions/inv_logit.hpp>
 #include <stan/math/opencl/kernels/device_functions/inv_Phi.hpp>
 #include <stan/math/opencl/kernels/device_functions/inv_square.hpp>
@@ -313,17 +314,19 @@ ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_logit,
                                  opencl_kernels::inv_logit_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::log1m_device_function,
                                  opencl_kernels::logit_device_function)
-ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::phi_device_function)
+ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::erfcx_device_function,
+                                 opencl_kernels::phi_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi_approx,
                                  opencl_kernels::inv_logit_device_function,
                                  opencl_kernels::phi_approx_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(
-    std_normal_lcdf_scaled_impl,
+    std_normal_lcdf_impl, opencl_kernels::erfcx_device_function,
     opencl_kernels::std_normal_lcdf_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(
-    std_normal_lcdf_dscaled_impl,
+    std_normal_lcdf_derivative, opencl_kernels::erfcx_device_function,
     opencl_kernels::std_normal_lcdf_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_Phi, opencl_kernels::log1m_device_function,
+                                 opencl_kernels::erfcx_device_function,
                                  opencl_kernels::phi_device_function,
                                  opencl_kernels::inv_phi_device_function)
 ADD_UNARY_FUNCTION_WITH_INCLUDES(
@@ -355,6 +358,32 @@ ADD_BINARY_FUNCTION_WITH_INCLUDES(copysign)
 
 ADD_BINARY_FUNCTION_WITH_INCLUDES(
     erfcx_derivative, stan::math::opencl_kernels::erfcx_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lcdf_value,
+    stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lcdf_d_z, stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lcdf_d_v, stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lcdf_d_s, stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lccdf_value,
+    stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lccdf_d_z, stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lccdf_d_v, stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
+ADD_BINARY_FUNCTION_WITH_INCLUDES(
+    exp_mod_normal_lccdf_d_s, stan::math::opencl_kernels::erfcx_device_function,
+    stan::math::opencl_kernels::exp_mod_normal_lcdf_device_function)
 ADD_BINARY_FUNCTION_WITH_INCLUDES(
     beta, stan::math::opencl_kernels::beta_device_function)
 ADD_BINARY_FUNCTION_WITH_INCLUDES(

stan/math/opencl/kernels/device_functions/Phi.hpp

@@ -21,15 +21,31 @@ static constexpr const char* phi_device_function
            * @return Phi(x)
            */
           inline double Phi(double x) {
-            if (x < -37.5) {
-              return 0;
-            } else if (x < -5.0) {
-              return 0.5 * erfc(-M_SQRT1_2 * x);
-            } else if (x > 8.25) {
-              return 1;
-            } else {
-              return 0.5 * (1.0 + erf(M_SQRT1_2 * x));
+            // erfcx_cody_middle and erfcx_cody_tail come from the erfcx
+            // device function, which the registration includes first.
+            const double u = M_SQRT1_2 * x;
+            const double a = fabs(u);
+            if (a < 0.46875) {
+              const double s = u * u;
+              const double num
+                  = (((1.85777706184603153e-1 * s + 3.16112374387056560) * s
+                      + 1.13864154151050156e2)
+                         * s
+                     + 3.77485237685302021e2)
+                        * s
+                    + 3.20937758913846947e3;
+              const double den
+                  = (((s + 2.36012909523441209e1) * s + 2.44024637934444173e2)
+                         * s
+                     + 1.28261652607737228e3)
+                        * s
+                    + 2.84423683343917062e3;
+              return 0.5 + 0.5 * (u * num / den);
             }
+            const double tail
+                = 0.5 * exp(-(0.5 * x) * x)
+                  * (a >= 4.0 ? erfcx_cody_tail(a) : erfcx_cody_middle(a));
+            return x < 0.0 ? tail : 1.0 - tail;
           }
           // \cond
           ) "\n#endif\n";  // NOLINT