diff --git a/doc/distributions/dist_tutorial.qbk b/doc/distributions/dist_tutorial.qbk
index ee28dbdca5..80e8e2e458 100644
--- a/doc/distributions/dist_tutorial.qbk
+++ b/doc/distributions/dist_tutorial.qbk
@@ -128,12 +128,49 @@ And quantiles are just the same:
    quantile(my_dist, p);  // Returns the value of the random variable x
                           // such that cdf(my_dist, x) == p.
 
+As are the logcdf (Natural log of the Cumulative Distribution Function):
+
+   logcdf(my_dist, x); // Returns logcdf at at point x of distribution my_dist.
+
+And the logpdf (Natural log of the Probability Density Function):
+
+   logpdf(my_dist, x); // Returns logpdf at point x of distribution my_dist.
+
 If you're wondering why these aren't member functions, it's to
 make the library more easily extensible: if you want to add additional
 generic operations - let's say the /n'th moment/ - then all you have to
 do is add the appropriate non-member functions, overloaded for each
 implemented distribution type.
 
+The logcdf and logpdf functions are minimally calculated with log(cdf(my_dist, x)),
+and log(pdf(my_dist, x)) respectively. The following distributions have specialized
+implementations of the logcdf:
+
+* Exponential
+* Extreme Value
+* Geometric
+* Laplace
+* Logistic
+* Pareto
+* Rayleigh
+* Weibull
+
+And the following distributions have specialized implementations of logpdf:
+
+* Exponential
+* Extreme Value
+* Gamma
+* Inverse Gamma
+* Inverse Gaussian
+* Laplace
+* Normal
+* Poisson
+* Rayleigh
+* Weibull
+
+These above listed specialized implementations allow a higher degree of precision
+than can be obtained through the naive generic method.
+
 [tip
 
 [*Random numbers that approximate Quantiles of Distributions]
diff --git a/doc/distributions/exponential.qbk b/doc/distributions/exponential.qbk
index 043818b4a4..b47f0b8bcd 100644
--- a/doc/distributions/exponential.qbk
+++ b/doc/distributions/exponential.qbk
@@ -57,6 +57,9 @@ that are generic to all distributions are supported: __usual_accessors.
 
 The domain of the random variable is \[0, +[infin]\].
 
+In this distribution the implementation of both `logcdf`, and `logpdf` are specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The exponential distribution is implemented in terms of the 
@@ -71,7 +74,9 @@ In the following table [lambda] is the parameter lambda of the distribution,
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = [lambda] * exp(-[lambda] * x) ]]
+[[logpdf][log(pdf) = -expm1(-x * [lambda]) ]]
 [[cdf][Using the relation: p = 1 - exp(-x * [lambda]) = -expm1(-x * [lambda]) ]]
+[[logcdf][log(cdf) = log1p(-exp(-x * [lambda])) ]]
 [[cdf complement][Using the relation: q = exp(-x * [lambda]) ]]
 [[quantile][Using the relation: x = -log(1-p) / [lambda] = -log1p(-p) / [lambda]]]
 [[quantile from the complement][Using the relation: x = -log(q) / [lambda]]]
diff --git a/doc/distributions/extreme_value.qbk b/doc/distributions/extreme_value.qbk
index 314917ebc1..43f68546f9 100644
--- a/doc/distributions/extreme_value.qbk
+++ b/doc/distributions/extreme_value.qbk
@@ -81,6 +81,9 @@ that are generic to all distributions are supported: __usual_accessors.
 
 The domain of the random parameter is \[-[infin], +[infin]\].
 
+In this distribution the implementation of both `logcdf`, and `logpdf` are specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The extreme value distribution is implemented in terms of the 
@@ -96,7 +99,9 @@ In the following table:
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = exp((a-x)/b) * exp(-exp((a-x)/b)) / b ]]
+[[logpdf][log(pdf) = log(1/b) + e - exp(e) ]]
 [[cdf][Using the relation: p = exp(-exp((a-x)/b)) ]]
+[[logcdf][log(cdf) = -exp((a-x)/b) ]]
 [[cdf complement][Using the relation: q = -expm1(-exp((a-x)/b)) ]]
 [[quantile][Using the relation: a - log(-log(p)) * b]]
 [[quantile from the complement][Using the relation: a - log(-log1p(-q)) * b]]
diff --git a/doc/distributions/gamma.qbk b/doc/distributions/gamma.qbk
index eefcc84a0c..56be002194 100644
--- a/doc/distributions/gamma.qbk
+++ b/doc/distributions/gamma.qbk
@@ -99,6 +99,9 @@ distributions are supported: __usual_accessors.
 
 The domain of the random variable is \[0,+[infin]\].
 
+In this distribution the implementation of `logpdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The gamma distribution is implemented in terms of the 
@@ -115,6 +118,7 @@ and /q = 1-p/.
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = __gamma_p_derivative(k, x / [theta]) / [theta] ]]
+[[logpdf][log(pdf) = -k*log([theta]) + (k-1)*log(x) - lgamma(k) - (x/[theta]) ]]
 [[cdf][Using the relation: p = __gamma_p(k, x / [theta]) ]]
 [[cdf complement][Using the relation: q = __gamma_q(k, x / [theta]) ]]
 [[quantile][Using the relation: x = [theta] * __gamma_p_inv(k, p) ]]
diff --git a/doc/distributions/geometric.qbk b/doc/distributions/geometric.qbk
index 7aa1a33439..cf1e2bff08 100644
--- a/doc/distributions/geometric.qbk
+++ b/doc/distributions/geometric.qbk
@@ -303,6 +303,9 @@ the context of this distribution:
    ``quantile(complement(geometric(p), P))``]]
 ]
 
+In this distribution the implementation of `logcdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 This distribution is implemented using the pow and exp functions, so most results
@@ -322,6 +325,7 @@ the expected number of failures using the quantile.
 [[Function][Implementation Notes]]
 [[pdf][pdf =  p * pow(q, k)]]
 [[cdf][cdf = 1 - q[super k=1]]]
+[[logcdf][log(cdf) = log1p(-exp(log1p(-p) * (k+1)))]]
 [[cdf complement][exp(log1p(-p) * (k+1))]]
 [[quantile][k = log1p(-x) / log1p(-p) -1]]
 [[quantile from the complement][k = log(x) / log1p(-p) -1]]
diff --git a/doc/distributions/inverse_gamma.qbk b/doc/distributions/inverse_gamma.qbk
index 8fccbc19c4..51d6538846 100644
--- a/doc/distributions/inverse_gamma.qbk
+++ b/doc/distributions/inverse_gamma.qbk
@@ -87,6 +87,9 @@ The domain of the random variate is \[0,+[infin]\].
 [note Unlike some definitions, this implementation supports a random variate 
 equal to zero as a special case, returning zero for pdf and cdf.]
 
+In this distribution the implementation of `logpdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The inverse gamma distribution is implemented in terms of the 
@@ -99,12 +102,13 @@ But in general, inverse_gamma results are accurate to a few epsilon,
 [h4 Implementation]
 
 In the following table [alpha] is the shape parameter of the distribution, 
-[alpha] is its scale parameter, /x/ is the random variate, /p/ is the probability
+[beta] is its scale parameter, /x/ is the random variate, /p/ is the probability
 and /q = 1-p/.
 
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = __gamma_p_derivative([alpha], [beta]/ x, [beta]) / x * x ]]
+[[logpdf][log(pdf) = [alpha] * log([beta]) + (-[alpha]-1)*log(x) - lgamma([alpha]) - ([beta]/x) ]]
 [[cdf][Using the relation: p = __gamma_q([alpha], [beta] / x) ]]
 [[cdf complement][Using the relation: q = __gamma_p([alpha], [beta] / x) ]]
 [[quantile][Using the relation: x = [beta]/ __gamma_q_inv([alpha], p) ]]
diff --git a/doc/distributions/inverse_gaussian.qbk b/doc/distributions/inverse_gaussian.qbk
index c5b824385f..6187a3835b 100644
--- a/doc/distributions/inverse_gaussian.qbk
+++ b/doc/distributions/inverse_gaussian.qbk
@@ -115,6 +115,9 @@ The domain of the random variate is \[0,+[infin]).
 [note Unlike some definitions, this implementation supports a random variate 
 equal to zero as a special case, returning zero for both pdf and cdf.]
 
+In this distribution the implementation of `logpdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The inverse_gaussian distribution is implemented in terms of the 
@@ -134,6 +137,7 @@ are used for the inverse gaussian function.
 [table
 [[Function] [Implementation Notes] ]
 [[pdf] [ [sqrt]([lambda]/ 2[pi]x[super 3]) e[super -[lambda](x - [mu])[sup2]/ 2[mu][sup2]x]]]
+[[logpdf] [log(pdf) = (-[lambda]*pow([mu]-x, 2)/(x*[mu][super 2]) + log([lambda]) - 3*log(x) - log(2*[pi])) / 2 ]]
 [[cdf][ [Phi]{[sqrt]([lambda]/x) (x/[mu]-1)} + e[super 2[mu]/[lambda]] [Phi]{-[sqrt]([lambda]/[mu]) (1+x/[mu])} ]]
 [[cdf complement] [using complement of [Phi] above.] ]
 [[quantile][No closed form known. Estimated using a guess refined by Newton-Raphson iteration.]]
diff --git a/doc/distributions/laplace.qbk b/doc/distributions/laplace.qbk
index 93327e0228..1afa7f4e1e 100644
--- a/doc/distributions/laplace.qbk
+++ b/doc/distributions/laplace.qbk
@@ -76,6 +76,9 @@ distributions are supported: __usual_accessors.
 
 The domain of the random variable is \[-[infin],+[infin]\].
 
+In this distribution the implementation of both `logcdf`, and `logpdf` are specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The laplace distribution is implemented in terms of the
@@ -90,11 +93,19 @@ and its complement /q = 1-p/.
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = e[super -abs(x-[mu]) \/ [sigma]] \/ (2 * [sigma]) ]]
+[[logpdf][log(pdf) = -abs(x-[mu])/[sigma] - log([sigma]) - log(2) ]]
 [[cdf][Using the relations:
 
 x <  [mu] : p =  e[super (x-[mu])/[sigma] ] \/ [sigma]
 
 x >= [mu] : p =  1 - e[super ([mu]-x)/[sigma] ] \/ [sigma]
+]]
+[[logcdf][log(cdf) = 
+
+x <  [mu] : p = ((x - [mu]) / [sigma]) - ln(2)
+
+x >= [mu] : p = log1p(-exp(([mu]-x) / [sigma]) / 2)
+
 ]]
 [[cdf complement][Using the relation:
 
diff --git a/doc/distributions/logistic.qbk b/doc/distributions/logistic.qbk
index 0a22b48d42..84b3a24f1d 100644
--- a/doc/distributions/logistic.qbk
+++ b/doc/distributions/logistic.qbk
@@ -67,6 +67,9 @@ At `p=1` and `p=0`, the quantile function returns the result of
 quantile function returns the result of -__overflow_error and 
 +__overflow_error respectively. 
 
+In this distribution the implementation of `logcdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The logistic distribution is implemented in terms of the `std::exp` 
@@ -82,6 +85,7 @@ in such cases, only a low /absolute error/ can be guaranteed.
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = e[super -(x-u)/s] / (s*(1+e[super -(x-u)/s])[super 2])]]
 [[cdf][Using the relation: p = 1/(1+e[super -(x-u)/s])]]
+[[logcdf][log(cdf) = -log1p(exp((u-x)/s)) ]]
 [[cdf complement][Using the relation: q = 1/(1+e[super (x-u)/s])]]
 [[quantile][Using the relation: x = u - s*log(1/p-1)]]
 [[quantile from the complement][Using the relation: x = u + s*log(p/1-p)]]
diff --git a/doc/distributions/non_members.qbk b/doc/distributions/non_members.qbk
index 9c900bf127..99728b8700 100644
--- a/doc/distributions/non_members.qbk
+++ b/doc/distributions/non_members.qbk
@@ -17,6 +17,8 @@ to go straight to the function you want if you already know its name.
 * __hazard.
 * __kurtosis.
 * __kurtosis_excess
+* __logcdf.
+* __logpdf.
 * __mean.
 * __median.
 * __mode.
@@ -349,6 +351,20 @@ Kurtosis excess can have a value from -2 to + infinity.
 
 The kurtosis excess of a normal distribution is zero.
 
+[h4:logcdf Natural Log of the Cumulative Distribution Function]
+
+   template <class RealType, ``__Policy``>
+   RealType logcdf(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
+
+Returns the natural log of the CDF of distribution /dist/.
+
+[h4:logpdf Natural Log of the Probability Density Function]
+
+   template <class RealType, ``__Policy``>
+   RealType logcdf(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
+
+Returns the natural log of the CDF of distribution /dist/.
+
 [h4:cdfPQ P and Q]
 
 The terms P and Q are sometimes used to refer to the __cdf
diff --git a/doc/distributions/normal.qbk b/doc/distributions/normal.qbk
index 52ac44e96b..eefe4e967b 100644
--- a/doc/distributions/normal.qbk
+++ b/doc/distributions/normal.qbk
@@ -97,6 +97,7 @@ and /s/ is its standard deviation.
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = e[super -(x-m)[super 2]\/(2s[super 2])] \/ (s * sqrt(2*pi)) ]]
+[[logpdf][log(pdf) = -log(s) - log(2*[pi])/2 - (x-mean)[super 2]/(2*s[super 2]) ]]
 [[cdf][Using the relation: p = 0.5 * __erfc(-(x-m)/(s*sqrt(2))) ]]
 [[cdf complement][Using the relation: q = 0.5 * __erfc((x-m)/(s*sqrt(2))) ]]
 [[quantile][Using the relation: x = m - s * sqrt(2) * __erfc_inv(2*p)]]
diff --git a/doc/distributions/pareto.qbk b/doc/distributions/pareto.qbk
index fcc7eee425..6a89b0dbc5 100644
--- a/doc/distributions/pareto.qbk
+++ b/doc/distributions/pareto.qbk
@@ -73,6 +73,9 @@ distributions are supported: __usual_accessors.
 
 The supported domain of the random variable is \[scale, [infin]\].
 
+In this distribution the implementation of `logcdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The Pareto distribution is implemented in terms of the
@@ -91,6 +94,7 @@ and its complement /q = 1-p/.
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf p = [alpha][beta][super [alpha]]/x[super [alpha] +1] ]]
 [[cdf][Using the relation: cdf p = 1 - ([beta] / x)[super [alpha]] ]]
+[[logcdf][log(cdf) = log1p(-pow([beta]/x, [alpha])) ]]
 [[cdf complement][Using the relation: q = 1 - p = -([beta] / x)[super [alpha]] ]]
 [[quantile][Using the relation: x = [beta] / (1 - p)[super 1/[alpha]] ]]
 [[quantile from the complement][Using the relation: x =  [beta] / (q)[super 1/[alpha]] ]]
diff --git a/doc/distributions/poisson.qbk b/doc/distributions/poisson.qbk
index 50a3b12cea..01cfb13cc5 100644
--- a/doc/distributions/poisson.qbk
+++ b/doc/distributions/poisson.qbk
@@ -62,6 +62,9 @@ distributions are supported: __usual_accessors.
 
 The domain of the random variable is \[0, [infin]\].
 
+In this distribution the implementation of `logpdf` is specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The Poisson distribution is implemented in terms of the 
@@ -81,6 +84,7 @@ In the following table [lambda] is the mean of the distribution,
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = e[super -[lambda]] [lambda][super k] \/ k! ]]
+[[logpdf][log(pdf) = -lgamma(k+1) + k*log([lambda]) - [lambda] if k > 0 and [lambda] > 0 ]]
 [[cdf][Using the relation: p = [Gamma](k+1, [lambda]) \/ k! = __gamma_q(k+1, [lambda])]]
 [[cdf complement][Using the relation: q = __gamma_p(k+1, [lambda]) ]]
 [[quantile][Using the relation: k = __gamma_q_inva([lambda], p) - 1]]
diff --git a/doc/distributions/rayleigh.qbk b/doc/distributions/rayleigh.qbk
index 5fd6fe44c2..609b654da7 100644
--- a/doc/distributions/rayleigh.qbk
+++ b/doc/distributions/rayleigh.qbk
@@ -77,6 +77,9 @@ distributions are supported: __usual_accessors.
 
 The domain of the random variable is \[0, max_value\].
 
+In this distribution the implementation of both `logcdf`, and `logpdf` are specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The Rayleigh distribution is implemented in terms of the
@@ -92,7 +95,9 @@ In the following table [sigma] is the sigma parameter of the distribution,
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = x * exp(-x[super 2])/2 [sigma][super 2] ]]
+[[logpdf][log(pdf) = -(x[super 2])/(2*[sigma][super 2]) - 2*log([sigma]) + log(x) ]]
 [[cdf][Using the relation: p = 1 - exp(-x[super 2]/2) [sigma][super 2]= -__expm1(-x[super 2]/2) [sigma][super 2]]]
+[[logcdf][log(cdf) = log1p(-exp(-x[super 2] / (2*[sigma][super 2]))) ]]
 [[cdf complement][Using the relation: q =  exp(-x[super 2]/ 2) * [sigma][super 2] ]]
 [[quantile][Using the relation: x = sqrt(-2 * [sigma] [super 2]) * log(1 - p)) = sqrt(-2 * [sigma] [super 2]) * __log1p(-p))]]
 [[quantile from the complement][Using the relation: x = sqrt(-2 * [sigma] [super 2]) * log(q)) ]]
diff --git a/doc/distributions/weibull.qbk b/doc/distributions/weibull.qbk
index 95c9e461e2..3ec18b18d9 100644
--- a/doc/distributions/weibull.qbk
+++ b/doc/distributions/weibull.qbk
@@ -88,6 +88,9 @@ distributions are supported: __usual_accessors.
 
 The domain of the random variable is \[0, [infin]\].
 
+In this distribution the implementation of both `logcdf`, and `logpdf` are specialized
+to improve numerical accuracy.
+
 [h4 Accuracy]
 
 The Weibull distribution is implemented in terms of the 
@@ -104,7 +107,9 @@ and /q = 1-p/.
 [table
 [[Function][Implementation Notes]]
 [[pdf][Using the relation: pdf = [alpha][beta][super -[alpha] ]x[super [alpha] - 1] e[super -(x/beta)[super alpha]] ]]
+[[logpdf][log(pdf) = log([alpha]) - [alpha] * log([beta]) + log(x) * ([alpha]-1) - pow(x/[beta], [alpha]) ]]
 [[cdf][Using the relation: p = -__expm1(-(x\/[beta])[super [alpha]]) ]]
+[[logcdf][log(cdf) = log1p(-exp(-pow(x / [beta], [alpha]))) ]]
 [[cdf complement][Using the relation: q = e[super -(x\/[beta])[super [alpha]]] ]]
 [[quantile][Using the relation: x = [beta] * (-__log1p(-p))[super 1\/[alpha]] ]]
 [[quantile from the complement][Using the relation: x = [beta] * (-log(q))[super 1\/[alpha]] ]]
diff --git a/doc/math.qbk b/doc/math.qbk
index 4b8804dbbb..d6b90efb0b 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -424,7 +424,7 @@ and use the function's name as the link text.]
 
 
 [def __usual_accessors __cdf, __pdf, __quantile, __hazard,
-   __chf, __mean, __median, __mode, __variance, __sd, __skewness,
+   __chf, __logcdf, __logpdf, __mean, __median, __mode, __variance, __sd, __skewness,
    __kurtosis, __kurtosis_excess, __range and __support]
 
 [def __real_concept [link math_toolkit.real_concepts  real concept]]
diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp
index 55e2236f99..4d1fc46569 100644
--- a/include/boost/math/distributions/laplace.hpp
+++ b/include/boost/math/distributions/laplace.hpp
@@ -182,7 +182,7 @@ BOOST_MATH_GPU_ENABLED inline RealType logpdf(const laplace_distribution<RealTyp
    const RealType mu = dist.scale();
    const RealType b = dist.location();
 
-   // if b is 0 avoid divde by 0 error
+   // if b is 0 avoid divide by 0 error
    if(abs(b) < boost::math::numeric_limits<RealType>::epsilon())
    {
       result = log(pdf(dist, x));