diff --git a/src/docs/stan-reference/functions.tex b/src/docs/stan-reference/functions.tex
index 289158277b..467eaa32ed 100644
--- a/src/docs/stan-reference/functions.tex
+++ b/src/docs/stan-reference/functions.tex
@@ -616,7 +616,11 @@ \section{Probability-Related Functions}\label{Phi-function.section}
 The complementary error function of \farg{x}}
 %
 \fitemindexsort{real}{Phi}{real \farg{x}}{
-The cumulative unit normal density function of \farg{x}}{Phi}{phi}
+The cumulative unit normal density function of \farg{x};
+\code{Phi(\farg{x})} will underflow to 0 for \farg{x} below -37.5 and overflow
+to 1 for \farg{x} above 8.25;  derivatives will underflow to 0 below
+-27.5 and overflow to 1 above 27.5.
+}{Phi}{phi}
 %
 \fitemindexsort{real}{Phi\_approx}{real \farg{x}}{
 Fast approximation of the cumulative unit normal density
diff --git a/src/stan/math/functions/Phi.hpp b/src/stan/math/functions/Phi.hpp
index c2d0019bb1..46e8b2c0b9 100644
--- a/src/stan/math/functions/Phi.hpp
+++ b/src/stan/math/functions/Phi.hpp
@@ -20,13 +20,24 @@ namespace stan {
      * This function can be used to implement the inverse link function
      * for probit regression.  
      *
+     * Phi will underflow to 0 below -37.5 and overflow to 1 above 8
+     *
      * @param x Argument.
      * @return Probability random sample is less than or equal to argument. 
      */
     template <typename T>
     inline typename boost::math::tools::promote_args<T>::type
     Phi(const T x) {
-      return 0.5 * (1.0 + boost::math::erf(INV_SQRT_2 * x));
+      // overridden in fvar and var, so can hard-code boost versions
+      // here for scalars only
+      if (x < -37.5)
+        return 0;
+      else if (x < -5.0)
+        return 0.5 * boost::math::erfc(-INV_SQRT_2 * x);
+      else if (x > 8.25)
+        return 1;
+      else
+        return 0.5 * (1.0 + boost::math::erf(INV_SQRT_2 * x));
     }
 
   }
diff --git a/src/test/agrad/fwd/Phi_test.cpp b/src/test/agrad/fwd/Phi_test.cpp
index 95fc0b4e80..ac40f0c3a8 100644
--- a/src/test/agrad/fwd/Phi_test.cpp
+++ b/src/test/agrad/fwd/Phi_test.cpp
@@ -17,6 +17,18 @@ TEST(AgradFwdPhi,Fvar) {
   EXPECT_FLOAT_EQ(exp(stan::prob::normal_log<false>(1.0,0.0,1.0)),
                   Phi_x.d_);
 }
+TEST(AgradFwdPhi, FvarDerivUnderOverFlow) {
+  using stan::agrad::fvar;
+  fvar<double> x = -27.5;
+  x.d_ = 1.0;
+  fvar<double> Phi_x = Phi(x);
+  EXPECT_FLOAT_EQ(0, Phi_x.d_);
+
+  fvar<double> y = 27.5;
+  y.d_ = 1.0;
+  fvar<double> Phi_y = Phi(y);
+  EXPECT_FLOAT_EQ(0, Phi_y.d_);
+}
 TEST(AgradFwdPhi, FvarVar_1stDeriv) {
   using stan::agrad::fvar;
   using stan::agrad::var;
@@ -154,3 +166,122 @@ TEST(AgradFwdPhi, FvarFvarVar_3rdDeriv) {
   a.d_.d_.grad(p,g);
   EXPECT_FLOAT_EQ(0, g[0]);
 }
+
+double fun_value_of(double x) { return x; }
+// inefficient calls by value because neither const& nor & compile
+double fun_value_of(stan::agrad::fvar<double> x) {
+  return x.val();
+}
+double fun_value_of(stan::agrad::fvar<stan::agrad::fvar<double> > x) {
+  return x.val().val();
+}
+
+// tests calculating using R 3.0.2 Snow Leopard build (6558)
+template <typename T>
+void test_tails() {
+  using stan::agrad::Phi;
+
+  EXPECT_EQ(0, fun_value_of(Phi(T(-40))));
+
+  EXPECT_FLOAT_EQ(1, 4.60535300958196e-308 / fun_value_of(Phi(T(-37.5))));
+  EXPECT_FLOAT_EQ(1, 5.72557122252458e-300 / fun_value_of(Phi(T(-37))));
+  EXPECT_FLOAT_EQ(1, 5.54472571307484e-292 / fun_value_of(Phi(T(-36.5))));
+  EXPECT_FLOAT_EQ(1, 4.18262406579728e-284 / fun_value_of(Phi(T(-36))));
+  EXPECT_FLOAT_EQ(1, 2.45769154066194e-276 / fun_value_of(Phi(T(-35.5))));
+  EXPECT_FLOAT_EQ(1, 1.12491070647241e-268 / fun_value_of(Phi(T(-35))));
+  EXPECT_FLOAT_EQ(1, 4.01072896657726e-261 / fun_value_of(Phi(T(-34.5))));
+  EXPECT_FLOAT_EQ(1, 1.11389878557438e-253 / fun_value_of(Phi(T(-34))));
+  EXPECT_FLOAT_EQ(1, 2.40983869512039e-246 / fun_value_of(Phi(T(-33.5))));
+  EXPECT_FLOAT_EQ(1, 4.06118562091586e-239 / fun_value_of(Phi(T(-33))));
+  EXPECT_FLOAT_EQ(1, 5.33142435967881e-232 / fun_value_of(Phi(T(-32.5))));
+  EXPECT_FLOAT_EQ(1, 5.4520806035124e-225 / fun_value_of(Phi(T(-32))));
+  EXPECT_FLOAT_EQ(1, 4.34323260103177e-218 / fun_value_of(Phi(T(-31.5))));
+  EXPECT_FLOAT_EQ(1, 2.6952500812005e-211 / fun_value_of(Phi(T(-31))));
+  EXPECT_FLOAT_EQ(1, 1.30293791317808e-204 / fun_value_of(Phi(T(-30.5))));
+  EXPECT_FLOAT_EQ(1, 4.90671392714819e-198 / fun_value_of(Phi(T(-30))));
+  EXPECT_FLOAT_EQ(1, 1.43947455222918e-191 / fun_value_of(Phi(T(-29.5))));
+  EXPECT_FLOAT_EQ(1, 3.28978526670438e-185 / fun_value_of(Phi(T(-29))));
+  EXPECT_FLOAT_EQ(1, 5.85714125380634e-179 / fun_value_of(Phi(T(-28.5))));
+  EXPECT_FLOAT_EQ(1, 8.12386946965943e-173 / fun_value_of(Phi(T(-28))));
+  EXPECT_FLOAT_EQ(1, 8.77817055687808e-167 / fun_value_of(Phi(T(-27.5))));
+  EXPECT_FLOAT_EQ(1, 7.38948100688502e-161 / fun_value_of(Phi(T(-27))));
+  EXPECT_FLOAT_EQ(1, 4.84616266030332e-155 / fun_value_of(Phi(T(-26.5))));
+  EXPECT_FLOAT_EQ(1, 2.47606331550339e-149 / fun_value_of(Phi(T(-26))));
+  EXPECT_FLOAT_EQ(1, 9.85623651896393e-144 / fun_value_of(Phi(T(-25.5))));
+  EXPECT_FLOAT_EQ(1, 3.05669670638256e-138 / fun_value_of(Phi(T(-25))));
+  EXPECT_FLOAT_EQ(1, 7.38570686148941e-133 / fun_value_of(Phi(T(-24.5))));
+  EXPECT_FLOAT_EQ(1, 1.3903921185497e-127 / fun_value_of(Phi(T(-24))));
+  EXPECT_FLOAT_EQ(1, 2.03936756324998e-122 / fun_value_of(Phi(T(-23.5))));
+  EXPECT_FLOAT_EQ(1, 2.33063700622065e-117 / fun_value_of(Phi(T(-23))));
+  EXPECT_FLOAT_EQ(1, 2.07531079906636e-112 / fun_value_of(Phi(T(-22.5))));
+  EXPECT_FLOAT_EQ(1, 1.43989243514508e-107 / fun_value_of(Phi(T(-22))));
+  EXPECT_FLOAT_EQ(1, 7.78439707718263e-103 / fun_value_of(Phi(T(-21.5))));
+  EXPECT_FLOAT_EQ(1, 3.27927801897904e-98 / fun_value_of(Phi(T(-21))));
+  EXPECT_FLOAT_EQ(1, 1.0764673258791e-93 / fun_value_of(Phi(T(-20.5))));
+  EXPECT_FLOAT_EQ(1, 2.75362411860623e-89 / fun_value_of(Phi(T(-20))));
+  EXPECT_FLOAT_EQ(1, 5.48911547566041e-85 / fun_value_of(Phi(T(-19.5))));
+  EXPECT_FLOAT_EQ(1, 8.52722395263098e-81 / fun_value_of(Phi(T(-19))));
+  EXPECT_FLOAT_EQ(1, 1.03236986895633e-76 / fun_value_of(Phi(T(-18.5))));
+  EXPECT_FLOAT_EQ(1, 9.74094891893715e-73 / fun_value_of(Phi(T(-18))));
+  EXPECT_FLOAT_EQ(1, 7.16345876623504e-69 / fun_value_of(Phi(T(-17.5))));
+  EXPECT_FLOAT_EQ(1, 4.10599620209891e-65 / fun_value_of(Phi(T(-17))));
+  EXPECT_FLOAT_EQ(1, 1.83446300316473e-61 / fun_value_of(Phi(T(-16.5))));
+  EXPECT_FLOAT_EQ(1, 6.38875440053809e-58 / fun_value_of(Phi(T(-16))));
+  EXPECT_FLOAT_EQ(1, 1.73446079179387e-54 / fun_value_of(Phi(T(-15.5))));
+  EXPECT_FLOAT_EQ(1, 3.67096619931275e-51 / fun_value_of(Phi(T(-15))));
+  EXPECT_FLOAT_EQ(1, 6.05749476441522e-48 / fun_value_of(Phi(T(-14.5))));
+  EXPECT_FLOAT_EQ(1, 7.7935368191928e-45 / fun_value_of(Phi(T(-14))));
+  EXPECT_FLOAT_EQ(1, 7.81880730565789e-42 / fun_value_of(Phi(T(-13.5))));
+  EXPECT_FLOAT_EQ(1, 6.11716439954988e-39 / fun_value_of(Phi(T(-13))));
+  EXPECT_FLOAT_EQ(1, 3.73256429887771e-36 / fun_value_of(Phi(T(-12.5))));
+  EXPECT_FLOAT_EQ(1, 1.77648211207768e-33 / fun_value_of(Phi(T(-12))));
+  EXPECT_FLOAT_EQ(1, 6.59577144611367e-31 / fun_value_of(Phi(T(-11.5))));
+  EXPECT_FLOAT_EQ(1, 1.91065957449868e-28 / fun_value_of(Phi(T(-11))));
+  EXPECT_FLOAT_EQ(1, 4.31900631780923e-26 / fun_value_of(Phi(T(-10.5))));
+  EXPECT_FLOAT_EQ(1, 7.61985302416053e-24 / fun_value_of(Phi(T(-10))));
+  EXPECT_FLOAT_EQ(1, 1.04945150753626e-21 / fun_value_of(Phi(T(-9.5))));
+  EXPECT_FLOAT_EQ(1, 1.12858840595384e-19 / fun_value_of(Phi(T(-9))));
+  EXPECT_FLOAT_EQ(1, 9.47953482220332e-18 / fun_value_of(Phi(T(-8.5))));
+  EXPECT_FLOAT_EQ(1, 6.22096057427178e-16 / fun_value_of(Phi(T(-8))));
+  EXPECT_FLOAT_EQ(1, 3.1908916729109e-14 / fun_value_of(Phi(T(-7.5))));
+  EXPECT_FLOAT_EQ(1, 1.27981254388584e-12 / fun_value_of(Phi(T(-7))));
+  EXPECT_FLOAT_EQ(1, 4.01600058385912e-11 / fun_value_of(Phi(T(-6.5))));
+  EXPECT_FLOAT_EQ(1, 9.86587645037698e-10 / fun_value_of(Phi(T(-6))));
+  EXPECT_FLOAT_EQ(1, 1.89895624658877e-08 / fun_value_of(Phi(T(-5.5))));
+  EXPECT_FLOAT_EQ(1, 2.86651571879194e-07 / fun_value_of(Phi(T(-5))));
+  EXPECT_FLOAT_EQ(1, 3.39767312473006e-06 / fun_value_of(Phi(T(-4.5))));
+  EXPECT_FLOAT_EQ(1, 3.16712418331199e-05 / fun_value_of(Phi(T(-4))));
+  EXPECT_FLOAT_EQ(1, 0.000232629079035525 / fun_value_of(Phi(T(-3.5))));
+  EXPECT_FLOAT_EQ(1, 0.00134989803163009 / fun_value_of(Phi(T(-3))));
+  EXPECT_FLOAT_EQ(1, 0.00620966532577613 / fun_value_of(Phi(T(-2.5))));
+  EXPECT_FLOAT_EQ(1, 0.0227501319481792 / fun_value_of(Phi(T(-2))));
+  EXPECT_FLOAT_EQ(1, 0.0668072012688581 / fun_value_of(Phi(T(-1.5))));
+  EXPECT_FLOAT_EQ(1, 0.158655253931457 / fun_value_of(Phi(T(-1))));
+  EXPECT_FLOAT_EQ(1, 0.308537538725987 / fun_value_of(Phi(T(-0.5))));
+  EXPECT_FLOAT_EQ(1, 0.5 / fun_value_of(Phi(T(0))));
+  EXPECT_FLOAT_EQ(1, 0.691462461274013 / fun_value_of(Phi(T(0.5))));
+  EXPECT_FLOAT_EQ(1, 0.841344746068543 / fun_value_of(Phi(T(1))));
+  EXPECT_FLOAT_EQ(1, 0.933192798731142 / fun_value_of(Phi(T(1.5))));
+  EXPECT_FLOAT_EQ(1, 0.977249868051821 / fun_value_of(Phi(T(2))));
+  EXPECT_FLOAT_EQ(1, 0.993790334674224 / fun_value_of(Phi(T(2.5))));
+  EXPECT_FLOAT_EQ(1, 0.99865010196837 / fun_value_of(Phi(T(3))));
+  EXPECT_FLOAT_EQ(1, 0.999767370920964 / fun_value_of(Phi(T(3.5))));
+  EXPECT_FLOAT_EQ(1, 0.999968328758167 / fun_value_of(Phi(T(4))));
+  EXPECT_FLOAT_EQ(1, 0.999996602326875 / fun_value_of(Phi(T(4.5))));
+  EXPECT_FLOAT_EQ(1, 0.999999713348428 / fun_value_of(Phi(T(5))));
+  EXPECT_FLOAT_EQ(1, 0.999999981010438 / fun_value_of(Phi(T(5.5))));
+  EXPECT_FLOAT_EQ(1, 0.999999999013412 / fun_value_of(Phi(T(6))));
+  EXPECT_FLOAT_EQ(1, 0.99999999995984 / fun_value_of(Phi(T(6.5))));
+  EXPECT_FLOAT_EQ(1, 0.99999999999872 / fun_value_of(Phi(T(7))));
+  EXPECT_FLOAT_EQ(1, 0.999999999999968 / fun_value_of(Phi(T(7.5))));
+  EXPECT_FLOAT_EQ(1, 0.999999999999999 / fun_value_of(Phi(T(8))));
+  EXPECT_FLOAT_EQ(1, 1 / fun_value_of(Phi(T(8.5))));
+  EXPECT_FLOAT_EQ(1, 1 / fun_value_of(Phi(T(9))));
+  EXPECT_FLOAT_EQ(1, 1 / fun_value_of(Phi(T(9.5))));
+  EXPECT_FLOAT_EQ(1, 1 / fun_value_of(Phi(T(10))));
+}
+TEST(AgradRev, PhiTails) {
+  using stan::agrad::fvar;
+  test_tails<fvar<double> >();
+  test_tails<fvar<fvar<double> > >();
+}
diff --git a/src/test/agrad/rev/Phi_test.cpp b/src/test/agrad/rev/Phi_test.cpp
index 817cd34f32..4529250b23 100644
--- a/src/test/agrad/rev/Phi_test.cpp
+++ b/src/test/agrad/rev/Phi_test.cpp
@@ -10,12 +10,16 @@ TEST(AgradRev, Phi) {
   y_values.push_back(0.0);
   y_values.push_back(0.9);
   y_values.push_back(-5.0);
+  y_values.push_back(-27.5);
+  y_values.push_back(27.5);
 
   // d/dy = exp(normal_log(value_of(y), 0.0, 1.0))
   std::vector<double> dy_values;
   dy_values.push_back(0.3989423);
   dy_values.push_back(0.2660852);
   dy_values.push_back(1.4867195e-06);
+  dy_values.push_back(0);
+  dy_values.push_back(0);
 
   for (size_t i = 0; i < y_values.size(); i++) {
     var y, phi_y;
@@ -30,3 +34,108 @@ TEST(AgradRev, Phi) {
       << "y = " << y;
   }
 }
+
+// tests calculating using R 3.0.2 Snow Leopard build (6558)
+TEST(AgradRev, PhiTails) {
+  using stan::agrad::Phi;
+  using stan::agrad::var;
+
+  EXPECT_EQ(0, Phi(var(-40)).val());
+
+  EXPECT_FLOAT_EQ(1, 4.60535300958196e-308 / Phi(var(-37.5)).val());
+  EXPECT_FLOAT_EQ(1, 5.72557122252458e-300 / Phi(var(-37)).val());
+  EXPECT_FLOAT_EQ(1, 5.54472571307484e-292 / Phi(var(-36.5)).val());
+  EXPECT_FLOAT_EQ(1, 4.18262406579728e-284 / Phi(var(-36)).val());
+  EXPECT_FLOAT_EQ(1, 2.45769154066194e-276 / Phi(var(-35.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.12491070647241e-268 / Phi(var(-35)).val());
+  EXPECT_FLOAT_EQ(1, 4.01072896657726e-261 / Phi(var(-34.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.11389878557438e-253 / Phi(var(-34)).val());
+  EXPECT_FLOAT_EQ(1, 2.40983869512039e-246 / Phi(var(-33.5)).val());
+  EXPECT_FLOAT_EQ(1, 4.06118562091586e-239 / Phi(var(-33)).val());
+  EXPECT_FLOAT_EQ(1, 5.33142435967881e-232 / Phi(var(-32.5)).val());
+  EXPECT_FLOAT_EQ(1, 5.4520806035124e-225 / Phi(var(-32)).val());
+  EXPECT_FLOAT_EQ(1, 4.34323260103177e-218 / Phi(var(-31.5)).val());
+  EXPECT_FLOAT_EQ(1, 2.6952500812005e-211 / Phi(var(-31)).val());
+  EXPECT_FLOAT_EQ(1, 1.30293791317808e-204 / Phi(var(-30.5)).val());
+  EXPECT_FLOAT_EQ(1, 4.90671392714819e-198 / Phi(var(-30)).val());
+  EXPECT_FLOAT_EQ(1, 1.43947455222918e-191 / Phi(var(-29.5)).val());
+  EXPECT_FLOAT_EQ(1, 3.28978526670438e-185 / Phi(var(-29)).val());
+  EXPECT_FLOAT_EQ(1, 5.85714125380634e-179 / Phi(var(-28.5)).val());
+  EXPECT_FLOAT_EQ(1, 8.12386946965943e-173 / Phi(var(-28)).val());
+  EXPECT_FLOAT_EQ(1, 8.77817055687808e-167 / Phi(var(-27.5)).val());
+  EXPECT_FLOAT_EQ(1, 7.38948100688502e-161 / Phi(var(-27)).val());
+  EXPECT_FLOAT_EQ(1, 4.84616266030332e-155 / Phi(var(-26.5)).val());
+  EXPECT_FLOAT_EQ(1, 2.47606331550339e-149 / Phi(var(-26)).val());
+  EXPECT_FLOAT_EQ(1, 9.85623651896393e-144 / Phi(var(-25.5)).val());
+  EXPECT_FLOAT_EQ(1, 3.05669670638256e-138 / Phi(var(-25)).val());
+  EXPECT_FLOAT_EQ(1, 7.38570686148941e-133 / Phi(var(-24.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.3903921185497e-127 / Phi(var(-24)).val());
+  EXPECT_FLOAT_EQ(1, 2.03936756324998e-122 / Phi(var(-23.5)).val());
+  EXPECT_FLOAT_EQ(1, 2.33063700622065e-117 / Phi(var(-23)).val());
+  EXPECT_FLOAT_EQ(1, 2.07531079906636e-112 / Phi(var(-22.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.43989243514508e-107 / Phi(var(-22)).val());
+  EXPECT_FLOAT_EQ(1, 7.78439707718263e-103 / Phi(var(-21.5)).val());
+  EXPECT_FLOAT_EQ(1, 3.27927801897904e-98 / Phi(var(-21)).val());
+  EXPECT_FLOAT_EQ(1, 1.0764673258791e-93 / Phi(var(-20.5)).val());
+  EXPECT_FLOAT_EQ(1, 2.75362411860623e-89 / Phi(var(-20)).val());
+  EXPECT_FLOAT_EQ(1, 5.48911547566041e-85 / Phi(var(-19.5)).val());
+  EXPECT_FLOAT_EQ(1, 8.52722395263098e-81 / Phi(var(-19)).val());
+  EXPECT_FLOAT_EQ(1, 1.03236986895633e-76 / Phi(var(-18.5)).val());
+  EXPECT_FLOAT_EQ(1, 9.74094891893715e-73 / Phi(var(-18)).val());
+  EXPECT_FLOAT_EQ(1, 7.16345876623504e-69 / Phi(var(-17.5)).val());
+  EXPECT_FLOAT_EQ(1, 4.10599620209891e-65 / Phi(var(-17)).val());
+  EXPECT_FLOAT_EQ(1, 1.83446300316473e-61 / Phi(var(-16.5)).val());
+  EXPECT_FLOAT_EQ(1, 6.38875440053809e-58 / Phi(var(-16)).val());
+  EXPECT_FLOAT_EQ(1, 1.73446079179387e-54 / Phi(var(-15.5)).val());
+  EXPECT_FLOAT_EQ(1, 3.67096619931275e-51 / Phi(var(-15)).val());
+  EXPECT_FLOAT_EQ(1, 6.05749476441522e-48 / Phi(var(-14.5)).val());
+  EXPECT_FLOAT_EQ(1, 7.7935368191928e-45 / Phi(var(-14)).val());
+  EXPECT_FLOAT_EQ(1, 7.81880730565789e-42 / Phi(var(-13.5)).val());
+  EXPECT_FLOAT_EQ(1, 6.11716439954988e-39 / Phi(var(-13)).val());
+  EXPECT_FLOAT_EQ(1, 3.73256429887771e-36 / Phi(var(-12.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.77648211207768e-33 / Phi(var(-12)).val());
+  EXPECT_FLOAT_EQ(1, 6.59577144611367e-31 / Phi(var(-11.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.91065957449868e-28 / Phi(var(-11)).val());
+  EXPECT_FLOAT_EQ(1, 4.31900631780923e-26 / Phi(var(-10.5)).val());
+  EXPECT_FLOAT_EQ(1, 7.61985302416053e-24 / Phi(var(-10)).val());
+  EXPECT_FLOAT_EQ(1, 1.04945150753626e-21 / Phi(var(-9.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.12858840595384e-19 / Phi(var(-9)).val());
+  EXPECT_FLOAT_EQ(1, 9.47953482220332e-18 / Phi(var(-8.5)).val());
+  EXPECT_FLOAT_EQ(1, 6.22096057427178e-16 / Phi(var(-8)).val());
+  EXPECT_FLOAT_EQ(1, 3.1908916729109e-14 / Phi(var(-7.5)).val());
+  EXPECT_FLOAT_EQ(1, 1.27981254388584e-12 / Phi(var(-7)).val());
+  EXPECT_FLOAT_EQ(1, 4.01600058385912e-11 / Phi(var(-6.5)).val());
+  EXPECT_FLOAT_EQ(1, 9.86587645037698e-10 / Phi(var(-6)).val());
+  EXPECT_FLOAT_EQ(1, 1.89895624658877e-08 / Phi(var(-5.5)).val());
+  EXPECT_FLOAT_EQ(1, 2.86651571879194e-07 / Phi(var(-5)).val());
+  EXPECT_FLOAT_EQ(1, 3.39767312473006e-06 / Phi(var(-4.5)).val());
+  EXPECT_FLOAT_EQ(1, 3.16712418331199e-05 / Phi(var(-4)).val());
+  EXPECT_FLOAT_EQ(1, 0.000232629079035525 / Phi(var(-3.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.00134989803163009 / Phi(var(-3)).val());
+  EXPECT_FLOAT_EQ(1, 0.00620966532577613 / Phi(var(-2.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.0227501319481792 / Phi(var(-2)).val());
+  EXPECT_FLOAT_EQ(1, 0.0668072012688581 / Phi(var(-1.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.158655253931457 / Phi(var(-1)).val());
+  EXPECT_FLOAT_EQ(1, 0.308537538725987 / Phi(var(-0.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.5 / Phi(var(0)).val());
+  EXPECT_FLOAT_EQ(1, 0.691462461274013 / Phi(var(0.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.841344746068543 / Phi(var(1)).val());
+  EXPECT_FLOAT_EQ(1, 0.933192798731142 / Phi(var(1.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.977249868051821 / Phi(var(2)).val());
+  EXPECT_FLOAT_EQ(1, 0.993790334674224 / Phi(var(2.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.99865010196837 / Phi(var(3)).val());
+  EXPECT_FLOAT_EQ(1, 0.999767370920964 / Phi(var(3.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.999968328758167 / Phi(var(4)).val());
+  EXPECT_FLOAT_EQ(1, 0.999996602326875 / Phi(var(4.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.999999713348428 / Phi(var(5)).val());
+  EXPECT_FLOAT_EQ(1, 0.999999981010438 / Phi(var(5.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.999999999013412 / Phi(var(6)).val());
+  EXPECT_FLOAT_EQ(1, 0.99999999995984 / Phi(var(6.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.99999999999872 / Phi(var(7)).val());
+  EXPECT_FLOAT_EQ(1, 0.999999999999968 / Phi(var(7.5)).val());
+  EXPECT_FLOAT_EQ(1, 0.999999999999999 / Phi(var(8)).val());
+  EXPECT_FLOAT_EQ(1, 1 / Phi(var(8.5)).val());
+  EXPECT_FLOAT_EQ(1, 1 / Phi(var(9)).val());
+  EXPECT_FLOAT_EQ(1, 1 / Phi(var(9.5)).val());
+  EXPECT_FLOAT_EQ(1, 1 / Phi(var(10)).val());
+}
diff --git a/src/test/math/functions/Phi_test.cpp b/src/test/math/functions/Phi_test.cpp
index 50c8270f61..e93306f4d1 100644
--- a/src/test/math/functions/Phi_test.cpp
+++ b/src/test/math/functions/Phi_test.cpp
@@ -6,3 +6,107 @@ TEST(MathFunctions, Phi) {
   EXPECT_FLOAT_EQ(0.5 + 0.5 * boost::math::erf(0.9/std::sqrt(2.0)), stan::math::Phi(0.9));
   EXPECT_EQ(0.5 + 0.5 * boost::math::erf(-5.0/std::sqrt(2.0)), stan::math::Phi(-5.0));
 }
+
+// tests calculating using R 3.0.2 Snow Leopard build (6558)
+TEST(MathFunctions, PhiTails) {
+  using stan::math::Phi;
+
+  EXPECT_EQ(0, Phi(-40));
+
+  EXPECT_FLOAT_EQ(1, 4.60535300958196e-308 / Phi(-37.5));
+  EXPECT_FLOAT_EQ(1, 5.72557122252458e-300 / Phi(-37));
+  EXPECT_FLOAT_EQ(1, 5.54472571307484e-292 / Phi(-36.5));
+  EXPECT_FLOAT_EQ(1, 4.18262406579728e-284 / Phi(-36));
+  EXPECT_FLOAT_EQ(1, 2.45769154066194e-276 / Phi(-35.5));
+  EXPECT_FLOAT_EQ(1, 1.12491070647241e-268 / Phi(-35));
+  EXPECT_FLOAT_EQ(1, 4.01072896657726e-261 / Phi(-34.5));
+  EXPECT_FLOAT_EQ(1, 1.11389878557438e-253 / Phi(-34));
+  EXPECT_FLOAT_EQ(1, 2.40983869512039e-246 / Phi(-33.5));
+  EXPECT_FLOAT_EQ(1, 4.06118562091586e-239 / Phi(-33));
+  EXPECT_FLOAT_EQ(1, 5.33142435967881e-232 / Phi(-32.5));
+  EXPECT_FLOAT_EQ(1, 5.4520806035124e-225 / Phi(-32));
+  EXPECT_FLOAT_EQ(1, 4.34323260103177e-218 / Phi(-31.5));
+  EXPECT_FLOAT_EQ(1, 2.6952500812005e-211 / Phi(-31));
+  EXPECT_FLOAT_EQ(1, 1.30293791317808e-204 / Phi(-30.5));
+  EXPECT_FLOAT_EQ(1, 4.90671392714819e-198 / Phi(-30));
+  EXPECT_FLOAT_EQ(1, 1.43947455222918e-191 / Phi(-29.5));
+  EXPECT_FLOAT_EQ(1, 3.28978526670438e-185 / Phi(-29));
+  EXPECT_FLOAT_EQ(1, 5.85714125380634e-179 / Phi(-28.5));
+  EXPECT_FLOAT_EQ(1, 8.12386946965943e-173 / Phi(-28));
+  EXPECT_FLOAT_EQ(1, 8.77817055687808e-167 / Phi(-27.5));
+  EXPECT_FLOAT_EQ(1, 7.38948100688502e-161 / Phi(-27));
+  EXPECT_FLOAT_EQ(1, 4.84616266030332e-155 / Phi(-26.5));
+  EXPECT_FLOAT_EQ(1, 2.47606331550339e-149 / Phi(-26));
+  EXPECT_FLOAT_EQ(1, 9.85623651896393e-144 / Phi(-25.5));
+  EXPECT_FLOAT_EQ(1, 3.05669670638256e-138 / Phi(-25));
+  EXPECT_FLOAT_EQ(1, 7.38570686148941e-133 / Phi(-24.5));
+  EXPECT_FLOAT_EQ(1, 1.3903921185497e-127 / Phi(-24));
+  EXPECT_FLOAT_EQ(1, 2.03936756324998e-122 / Phi(-23.5));
+  EXPECT_FLOAT_EQ(1, 2.33063700622065e-117 / Phi(-23));
+  EXPECT_FLOAT_EQ(1, 2.07531079906636e-112 / Phi(-22.5));
+  EXPECT_FLOAT_EQ(1, 1.43989243514508e-107 / Phi(-22));
+  EXPECT_FLOAT_EQ(1, 7.78439707718263e-103 / Phi(-21.5));
+  EXPECT_FLOAT_EQ(1, 3.27927801897904e-98 / Phi(-21));
+  EXPECT_FLOAT_EQ(1, 1.0764673258791e-93 / Phi(-20.5));
+  EXPECT_FLOAT_EQ(1, 2.75362411860623e-89 / Phi(-20));
+  EXPECT_FLOAT_EQ(1, 5.48911547566041e-85 / Phi(-19.5));
+  EXPECT_FLOAT_EQ(1, 8.52722395263098e-81 / Phi(-19));
+  EXPECT_FLOAT_EQ(1, 1.03236986895633e-76 / Phi(-18.5));
+  EXPECT_FLOAT_EQ(1, 9.74094891893715e-73 / Phi(-18));
+  EXPECT_FLOAT_EQ(1, 7.16345876623504e-69 / Phi(-17.5));
+  EXPECT_FLOAT_EQ(1, 4.10599620209891e-65 / Phi(-17));
+  EXPECT_FLOAT_EQ(1, 1.83446300316473e-61 / Phi(-16.5));
+  EXPECT_FLOAT_EQ(1, 6.38875440053809e-58 / Phi(-16));
+  EXPECT_FLOAT_EQ(1, 1.73446079179387e-54 / Phi(-15.5));
+  EXPECT_FLOAT_EQ(1, 3.67096619931275e-51 / Phi(-15));
+  EXPECT_FLOAT_EQ(1, 6.05749476441522e-48 / Phi(-14.5));
+  EXPECT_FLOAT_EQ(1, 7.7935368191928e-45 / Phi(-14));
+  EXPECT_FLOAT_EQ(1, 7.81880730565789e-42 / Phi(-13.5));
+  EXPECT_FLOAT_EQ(1, 6.11716439954988e-39 / Phi(-13));
+  EXPECT_FLOAT_EQ(1, 3.73256429887771e-36 / Phi(-12.5));
+  EXPECT_FLOAT_EQ(1, 1.77648211207768e-33 / Phi(-12));
+  EXPECT_FLOAT_EQ(1, 6.59577144611367e-31 / Phi(-11.5));
+  EXPECT_FLOAT_EQ(1, 1.91065957449868e-28 / Phi(-11));
+  EXPECT_FLOAT_EQ(1, 4.31900631780923e-26 / Phi(-10.5));
+  EXPECT_FLOAT_EQ(1, 7.61985302416053e-24 / Phi(-10));
+  EXPECT_FLOAT_EQ(1, 1.04945150753626e-21 / Phi(-9.5));
+  EXPECT_FLOAT_EQ(1, 1.12858840595384e-19 / Phi(-9));
+  EXPECT_FLOAT_EQ(1, 9.47953482220332e-18 / Phi(-8.5));
+  EXPECT_FLOAT_EQ(1, 6.22096057427178e-16 / Phi(-8));
+  EXPECT_FLOAT_EQ(1, 3.1908916729109e-14 / Phi(-7.5));
+  EXPECT_FLOAT_EQ(1, 1.27981254388584e-12 / Phi(-7));
+  EXPECT_FLOAT_EQ(1, 4.01600058385912e-11 / Phi(-6.5));
+  EXPECT_FLOAT_EQ(1, 9.86587645037698e-10 / Phi(-6));
+  EXPECT_FLOAT_EQ(1, 1.89895624658877e-08 / Phi(-5.5));
+  EXPECT_FLOAT_EQ(1, 2.86651571879194e-07 / Phi(-5));
+  EXPECT_FLOAT_EQ(1, 3.39767312473006e-06 / Phi(-4.5));
+  EXPECT_FLOAT_EQ(1, 3.16712418331199e-05 / Phi(-4));
+  EXPECT_FLOAT_EQ(1, 0.000232629079035525 / Phi(-3.5));
+  EXPECT_FLOAT_EQ(1, 0.00134989803163009 / Phi(-3));
+  EXPECT_FLOAT_EQ(1, 0.00620966532577613 / Phi(-2.5));
+  EXPECT_FLOAT_EQ(1, 0.0227501319481792 / Phi(-2));
+  EXPECT_FLOAT_EQ(1, 0.0668072012688581 / Phi(-1.5));
+  EXPECT_FLOAT_EQ(1, 0.158655253931457 / Phi(-1));
+  EXPECT_FLOAT_EQ(1, 0.308537538725987 / Phi(-0.5));
+  EXPECT_FLOAT_EQ(1, 0.5 / Phi(0));
+  EXPECT_FLOAT_EQ(1, 0.691462461274013 / Phi(0.5));
+  EXPECT_FLOAT_EQ(1, 0.841344746068543 / Phi(1));
+  EXPECT_FLOAT_EQ(1, 0.933192798731142 / Phi(1.5));
+  EXPECT_FLOAT_EQ(1, 0.977249868051821 / Phi(2));
+  EXPECT_FLOAT_EQ(1, 0.993790334674224 / Phi(2.5));
+  EXPECT_FLOAT_EQ(1, 0.99865010196837 / Phi(3));
+  EXPECT_FLOAT_EQ(1, 0.999767370920964 / Phi(3.5));
+  EXPECT_FLOAT_EQ(1, 0.999968328758167 / Phi(4));
+  EXPECT_FLOAT_EQ(1, 0.999996602326875 / Phi(4.5));
+  EXPECT_FLOAT_EQ(1, 0.999999713348428 / Phi(5));
+  EXPECT_FLOAT_EQ(1, 0.999999981010438 / Phi(5.5));
+  EXPECT_FLOAT_EQ(1, 0.999999999013412 / Phi(6));
+  EXPECT_FLOAT_EQ(1, 0.99999999995984 / Phi(6.5));
+  EXPECT_FLOAT_EQ(1, 0.99999999999872 / Phi(7));
+  EXPECT_FLOAT_EQ(1, 0.999999999999968 / Phi(7.5));
+  EXPECT_FLOAT_EQ(1, 0.999999999999999 / Phi(8));
+  EXPECT_FLOAT_EQ(1, 1 / Phi(8.5));
+  EXPECT_FLOAT_EQ(1, 1 / Phi(9));
+  EXPECT_FLOAT_EQ(1, 1 / Phi(9.5));
+  EXPECT_FLOAT_EQ(1, 1 / Phi(10));
+}