diff --git a/src/functions-reference/real-valued_basic_functions.qmd b/src/functions-reference/real-valued_basic_functions.qmd
index e574d7c2e..ec3b29457 100644
--- a/src/functions-reference/real-valued_basic_functions.qmd
+++ b/src/functions-reference/real-valued_basic_functions.qmd
@@ -1107,7 +1107,11 @@ Quantile arguments below 1e-16 are untested; quantiles above 0.999999999 result
 `R` **`Phi_approx`**`(T x)`<br>\newline
 The fast approximation of the unit (may replace `Phi` for probit
 regression with maximum absolute error of 0.00014, see
-[@BowlingEtAl:2009] for details)
+[@BowlingEtAl:2009] for details). Its relative error is large in the
+tails. `Phi_approx` returns 0 for x below -21.6 and 1 for x above 7.1.
+For tail probabilities, use `Phi(x)` and `Phi(-x)`. For the logarithms of
+tail probabilities, use `std_normal_lcdf` and `std_normal_lccdf`, which
+stay accurate further into the tails.
 {{< since "2.0, vectorized in 2.13" >}}
 
 ### Other probability-related functions
diff --git a/src/functions-reference/unbounded_continuous_distributions.qmd b/src/functions-reference/unbounded_continuous_distributions.qmd
index dc2ff8715..c0c5b92ad 100644
--- a/src/functions-reference/unbounded_continuous_distributions.qmd
+++ b/src/functions-reference/unbounded_continuous_distributions.qmd
@@ -52,10 +52,10 @@ constant additive terms.
 
 `real` **`normal_cdf`**`(reals y | reals mu, reals sigma)`<br>\newline
 The cumulative normal distribution of y given location mu and scale
-sigma; normal_cdf will underflow to 0 for $\frac{{y}-{\mu}}{{\sigma}}$
-below -37.5 and overflow to 1 for $\frac{{y}-{\mu}}{{\sigma}}$ above
-8.25; the function `Phi_approx` is more robust in the tails, but must
-be scaled and translated for anything other than a standard normal.
+sigma. It returns 0 for $\frac{{y}-{\mu}}{{\sigma}}$ below about -38 and 1 for
+$\frac{{y}-{\mu}}{{\sigma}}$ above about 8.3, because double precision cannot
+represent probabilities closer to 0 or 1. For more accuracy in the tails, use
+`normal_lcdf` or `normal_lccdf`.
 {{< since 2.0 >}}
 
 <!-- real; normal_lcdf; (reals y | reals mu, reals sigma); -->
@@ -63,11 +63,9 @@ be scaled and translated for anything other than a standard normal.
 
 `real` **`normal_lcdf`**`(reals y | reals mu, reals sigma)`<br>\newline
 The log of the cumulative normal distribution of y given location mu
-and scale sigma; normal_lcdf will underflow to $-\infty$ for
-$\frac{{y}-{\mu}}{{\sigma}}$ below -37.5 and overflow to 0 for
-$\frac{{y}-{\mu}}{{\sigma}}$ above 8.25; `log(Phi_approx(...))` is more
-robust in the tails, but must be scaled and translated for anything other
-than a standard normal.
+and scale sigma. It is accurate in both tails. It returns 0 only for
+$\frac{{y}-{\mu}}{{\sigma}}$ above about 38.5, where the result is smaller in
+magnitude than the smallest double.
 {{< since 2.12 >}}
 
 <!-- real; normal_lccdf; (reals y | reals mu, reals sigma); -->
@@ -75,11 +73,9 @@ than a standard normal.
 
 `real` **`normal_lccdf`**`(reals y | reals mu, reals sigma)`<br>\newline
 The log of the complementary cumulative normal distribution of y given
-location mu and scale sigma; normal_lccdf will overflow to 0 for
-$\frac{{y}-{\mu}}{{\sigma}}$ below -37.5 and underflow to $-\infty$
-for $\frac{{y}-{\mu}}{{\sigma}}$ above 8.25; `log1m(Phi_approx(...))` is
-more robust in the tails, but must be scaled and translated for anything
-other than a standard normal.
+location mu and scale sigma. It is accurate in both tails. It returns 0
+only for $\frac{{y}-{\mu}}{{\sigma}}$ below about -38.5, where the result is
+smaller in magnitude than the smallest double.
 {{< since 2.15 >}}
 
 <!-- R; normal_rng; (reals mu, reals sigma); -->
@@ -141,28 +137,28 @@ of y dropping constant additive terms.
 \index{{\tt \bfseries std\_normal\_cdf }!{\tt (reals y): real}|hyperpage}
 
 `real` **`std_normal_cdf`**`(reals y)`<br>\newline
-The cumulative standard normal distribution of y; std_normal_cdf will
-underflow to 0 for $y$ below -37.5 and overflow to 1 for $y$ above 8.25;
-the function `Phi_approx` is more robust in the tails.
+The cumulative standard normal distribution of y. It returns 0 for $y$
+below about -38 and 1 for $y$ above about 8.3, because double precision
+cannot represent probabilities closer to 0 or 1. For more accuracy in the
+tails, use `std_normal_lcdf` or `std_normal_lccdf`.
 {{< since 2.21 >}}
 
 <!-- real; std_normal_lcdf; (reals y); -->
 \index{{\tt \bfseries std\_normal\_lcdf }!{\tt (reals y): real}|hyperpage}
 
 `real` **`std_normal_lcdf`**`(reals y)`<br>\newline
-The log of the cumulative standard normal distribution of y; std_normal_lcdf
-will underflow to $-\infty$ for $y$ below -37.5 and overflow to 0 for $y$
-above 8.25; `log(Phi_approx(...))` is more robust in the tails.
+The log of the cumulative standard normal distribution of y. It is
+accurate in both tails. It returns 0 only for $y$ above about 38.5, where
+the result is smaller in magnitude than the smallest double.
 {{< since 2.21 >}}
 
 <!-- real; std_normal_lccdf; (reals y); -->
 \index{{\tt \bfseries std\_normal\_lccdf }!{\tt (reals y): real}|hyperpage}
 
 `real` **`std_normal_lccdf`**`(reals y)`<br>\newline
-The log of the complementary cumulative standard normal distribution of y;
-std_normal_lccdf will overflow to 0 for $y$ below -37.5 and underflow to
-$-\infty$ for $y$ above 8.25; `log1m(Phi_approx(...))` is more robust in the
-tails.
+The log of the complementary cumulative standard normal distribution of
+y. It is accurate in both tails. It returns 0 only for $y$ below about
+-38.5, where the result is smaller in magnitude than the smallest double.
 {{< since 2.21 >}}
 
 <!-- R; std_normal_qf; (T x); -->