First pass vectorise

Author: andrjohns

stan/math/prim/fun/inv_Phi.hpp

@@ -43,8 +43,8 @@ inline double inv_Phi_lambda(double p) {
          1.9715909503065514427e+03, 1.3731693765509461125e+04,
          4.5921953931549871457e+04, 6.7265770927008700853e+04,
          3.3430575583588128105e+04, 2.5090809287301226727e+03};
-  static constexpr double b[7]
-      = {4.2313330701600911252e+01, 6.8718700749205790830e+02,
+  static constexpr double b[8]
+      = {1.0, 4.2313330701600911252e+01, 6.8718700749205790830e+02,
          5.3941960214247511077e+03, 2.1213794301586595867e+04,
          3.9307895800092710610e+04, 2.8729085735721942674e+04,
          5.2264952788528545610e+03};
@@ -53,8 +53,8 @@ inline double inv_Phi_lambda(double p) {
          5.76949722146069140550e+00, 3.64784832476320460504e+00,
          1.27045825245236838258e+00, 2.41780725177450611770e-01,
          2.27238449892691845833e-02, 7.74545014278341407640e-04};
-  static constexpr double d[7]
-      = {2.05319162663775882187e+00, 1.67638483018380384940e+00,
+  static constexpr double d[8]
+      = {1.0, 2.05319162663775882187e+00, 1.67638483018380384940e+00,
          6.89767334985100004550e-01, 1.48103976427480074590e-01,
          1.51986665636164571966e-02, 5.47593808499534494600e-04,
          1.05075007164441684324e-09};
@@ -63,72 +63,49 @@ inline double inv_Phi_lambda(double p) {
          1.78482653991729133580e+00, 2.96560571828504891230e-01,
          2.65321895265761230930e-02, 1.24266094738807843860e-03,
          2.71155556874348757815e-05, 2.01033439929228813265e-07};
-  static constexpr double f[7]
-      = {5.99832206555887937690e-01, 1.36929880922735805310e-01,
+  static constexpr double f[8]
+      = {1.0, 5.99832206555887937690e-01, 1.36929880922735805310e-01,
          1.48753612908506148525e-02, 7.86869131145613259100e-04,
          1.84631831751005468180e-05, 1.42151175831644588870e-07,
          2.04426310338993978564e-15};
 
   double q = p - 0.5;
-  double r;
+  double r = q < 0 ? p : 1 - p;
+
+  if (r <= 0) {
+    return 0;
+  }
+
   double val;
+  double inner_r;
+  double pre_mult;
+  using Vector8d = Eigen::Matrix<double, 8, 1>;
+  Eigen::Map<const Vector8d> numerator_map(NULL);
+  Eigen::Map<const Vector8d> denonimator_map(NULL);
 
   if (std::fabs(q) <= .425) {
-    r = .180625 - square(q);
-    return q
-           * (((((((a[7] * r + a[6]) * r + a[5]) * r + a[4]) * r + a[3]) * r
-                + a[2])
-                   * r
-               + a[1])
-                  * r
-              + a[0])
-           / (((((((b[6] * r + b[5]) * r + b[4]) * r + b[3]) * r + b[2]) * r
-                + b[1])
-                   * r
-               + b[0])
-                  * r
-              + 1.0);
+    inner_r = .180625 - square(q);
+    pre_mult = q;
+
+    new (&numerator_map) Eigen::Map<const Vector8d>(a, 8);
+    new (&denonimator_map) Eigen::Map<const Vector8d>(b, 8);
   } else {
-    r = q < 0 ? p : 1 - p;
-
-    if (r <= 0)
-      return 0;
-
-    r = std::sqrt(-std::log(r));
-
-    if (r <= 5.0) {
-      r += -1.6;
-      val = (((((((c[7] * r + c[6]) * r + c[5]) * r + c[4]) * r + c[3]) * r
-               + c[2])
-                  * r
-              + c[1])
-                 * r
-             + c[0])
-            / (((((((d[6] * r + d[5]) * r + d[4]) * r + d[3]) * r + d[2]) * r
-                 + d[1])
-                    * r
-                + d[0])
-                   * r
-               + 1.0);
+
+    double temp_r = std::sqrt(-std::log(r));
+    if (temp_r <= 5.0) {
+      inner_r = temp_r - 1.6;
+      new (&numerator_map) Eigen::Map<const Vector8d>(c, 8);
+      new (&denonimator_map) Eigen::Map<const Vector8d>(d, 8);
     } else {
-      r -= 5.0;
-      val = (((((((e[7] * r + e[6]) * r + e[5]) * r + e[4]) * r + e[3]) * r
-               + e[2])
-                  * r
-              + e[1])
-                 * r
-             + e[0])
-            / (((((((f[6] * r + f[5]) * r + f[4]) * r + f[3]) * r + f[2]) * r
-                 + f[1])
-                    * r
-                + f[0])
-                   * r
-               + 1.0);
+      inner_r = temp_r - 5.0;
+      new (&numerator_map) Eigen::Map<const Vector8d>(e, 8);
+      new (&denonimator_map) Eigen::Map<const Vector8d>(f, 8);
     }
-    if (q < 0.0)
-      return -val;
+    pre_mult = q < 0 ? -1 : 1;
   }
-  return val;
+
+  Eigen::VectorXd r_pow = pow(inner_r, Eigen::ArrayXd::LinSpaced(8, 0, 7)) / 10.0;
+  return pre_mult * (numerator_map.dot(r_pow) * 10.0) / (denonimator_map.dot(r_pow) * 10.0);
 }
 }  // namespace internal