diff -u gsl-1.15/randist/beta.c gsl-sam/randist/beta.c
--- gsl-1.15/randist/beta.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/beta.c	2012-06-15 15:55:57.000000000 +0100
@@ -42,9 +42,15 @@
 double
 gsl_ran_beta_pdf (const double x, const double a, const double b)
 {
+  return exp(gsl_ran_beta_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_beta_lnpdf (const double x, const double a, const double b)
+{
   if (x < 0 || x > 1)
     {
-      return 0 ;
+      return GSL_NEGINF;
     }
   else 
     {
@@ -54,20 +60,13 @@
       double ga = gsl_sf_lngamma (a);
       double gb = gsl_sf_lngamma (b);
       
-      if (x == 0.0 || x == 1.0) 
-        {
-	  if (a > 1.0 && b > 1.0)
-	    {
-	      p = 0.0;
-	    }
-	  else
-	    {
-	      p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
-	    }
-        }
+      if ((x == 0.0 || x == 1.0) && a > 1.0 && b > 1.0)
+	{
+	  p = GSL_NEGINF;
+	}
       else
         {
-          p = exp (gab - ga - gb + log(x) * (a - 1)  + log1p(-x) * (b - 1));
+          p = gab - ga - gb + log (x) * (a - 1) + log1p (-x) * (b - 1);
         }
 
       return p;
diff -u gsl-1.15/randist/bigauss.c gsl-sam/randist/bigauss.c
--- gsl-1.15/randist/bigauss.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/bigauss.c	2012-06-15 13:38:21.000000000 +0100
@@ -61,10 +61,18 @@
                                 const double sigma_x, const double sigma_y,
                                 const double rho)
 {
+  return exp (gsl_ran_bivariate_gaussian_lnpdf (x, y, sigma_x, sigma_y, rho));
+}
+
+double
+gsl_ran_bivariate_gaussian_lnpdf (const double x, const double y, 
+				  const double sigma_x, const double sigma_y,
+				  const double rho)
+{
   double u = x / sigma_x ;
   double v = y / sigma_y ;
   double c = 1 - rho*rho ;
-  double p = (1 / (2 * M_PI * sigma_x * sigma_y * sqrt(c))) 
-    * exp (-(u * u - 2 * rho * u * v + v * v) / (2 * c));
+  double p = -log(2 * M_PI * sigma_x * sigma_y * sqrt(c))
+    - (u * u - 2 * rho * u * v + v * v) / (2 * c);
   return p;
 }
diff -u gsl-1.15/randist/binomial.c gsl-sam/randist/binomial.c
--- gsl-1.15/randist/binomial.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/binomial.c	2012-06-15 16:21:17.000000000 +0100
@@ -23,6 +23,7 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_math.h>
 
 /* The binomial distribution has the form,
 
@@ -72,9 +73,16 @@
 gsl_ran_binomial_pdf (const unsigned int k, const double p,
                       const unsigned int n)
 {
+  return exp (gsl_ran_binomial_lnpdf (k, p, n));
+}
+
+double
+gsl_ran_binomial_lnpdf (const unsigned int k, const double p,
+			const unsigned int n)
+{
   if (k > n)
     {
-      return 0;
+      return GSL_NEGINF;
     }
   else
     {
@@ -82,17 +90,16 @@
 
       if (p == 0) 
         {
-          P = (k == 0) ? 1 : 0;
+          P = (k == 0) ? 0 : GSL_NEGINF;
         }
       else if (p == 1)
         {
-          P = (k == n) ? 1 : 0;
+          P = (k == n) ? 0 : GSL_NEGINF;
         }
       else
         {
           double ln_Cnk = gsl_sf_lnchoose (n, k);
           P = ln_Cnk + k * log (p) + (n - k) * log1p (-p);
-          P = exp (P);
         }
 
       return P;
diff -u gsl-1.15/randist/cauchy.c gsl-sam/randist/cauchy.c
--- gsl-1.15/randist/cauchy.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/cauchy.c	2012-06-15 15:00:07.000000000 +0100
@@ -45,7 +45,13 @@
 double
 gsl_ran_cauchy_pdf (const double x, const double a)
 {
+  return exp (gsl_ran_cauchy_lnpdf (x, a));
+}
+
+double
+gsl_ran_cauchy_lnpdf (const double x, const double a)
+{
   double u = x / a;
-  double p = (1 / (M_PI * a)) / (1 + u * u);
+  double p = -log(M_PI * a) - log1p(u * u);
   return p;
 }
diff -u gsl-1.15/randist/chisq.c gsl-sam/randist/chisq.c
--- gsl-1.15/randist/chisq.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/chisq.c	2012-06-15 16:21:23.000000000 +0100
@@ -22,6 +22,7 @@
 #include <gsl/gsl_sf_gamma.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The chisq distribution has the form
 
@@ -39,22 +40,28 @@
 double
 gsl_ran_chisq_pdf (const double x, const double nu)
 {
+  return exp (gsl_ran_chisq_lnpdf (x, nu));
+}
+
+double
+gsl_ran_chisq_lnpdf (const double x, const double nu)
+{
   if (x < 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
       if(nu == 2.0)
         {
-          return exp(-x/2.0) / 2.0;
+          return (-x/2.0) - log (2.0);
         }
       else
         {
           double p;
           double lngamma = gsl_sf_lngamma (nu / 2);
 
-          p = exp ((nu / 2 - 1) * log (x/2) - x/2 - lngamma) / 2;
+          p = (nu / 2 - 1) * log (x/2) - x/2 - lngamma - log (2.0);
           return p;
         }
     }
diff -u gsl-1.15/randist/erlang.c gsl-sam/randist/erlang.c
--- gsl-1.15/randist/erlang.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/erlang.c	2012-06-15 16:21:53.000000000 +0100
@@ -22,6 +22,7 @@
 #include <gsl/gsl_sf_gamma.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The sum of N samples from an exponential distribution gives an
    Erlang distribution
@@ -39,16 +40,22 @@
 double
 gsl_ran_erlang_pdf (const double x, const double a, const double n)
 {
+  return exp (gsl_ran_erlang_lnpdf (x, a, n));
+}
+
+double
+gsl_ran_erlang_lnpdf (const double x, const double a, const double n)
+{
   if (x <= 0) 
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
       double p;
       double lngamma = gsl_sf_lngamma (n);
 
-      p = exp ((n - 1) * log (x/a) - x/a - lngamma) / a;
+      p = (n - 1) * log (x/a) - x/a - lngamma - log (a);
       return p;
     }
 }
diff -u gsl-1.15/randist/exponential.c gsl-sam/randist/exponential.c
--- gsl-1.15/randist/exponential.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/exponential.c	2012-06-15 16:02:41.000000000 +0100
@@ -40,13 +40,19 @@
 double
 gsl_ran_exponential_pdf (const double x, const double mu)
 {
+  return exp (gsl_ran_exponential_lnpdf (x, mu));
+}
+
+double
+gsl_ran_exponential_lnpdf (const double x, const double mu)
+{
   if (x < 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
-      double p = exp (-x/mu)/mu;
+      double p = (-x/mu) - log (mu);
       
       return p;
     }
diff -u gsl-1.15/randist/exppow.c gsl-sam/randist/exppow.c
--- gsl-1.15/randist/exppow.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/exppow.c	2012-06-15 14:59:19.000000000 +0100
@@ -115,8 +115,14 @@
 double
 gsl_ran_exppow_pdf (const double x, const double a, const double b)
 {
+  return exp (gsl_ran_exppow_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_exppow_lnpdf (const double x, const double a, const double b)
+{
   double p;
   double lngamma = gsl_sf_lngamma (1 + 1 / b);
-  p = (1 / (2 * a)) * exp (-pow (fabs (x / a), b) - lngamma);
+  p = -log(2 * a) + (-pow (fabs (x / a), b) - lngamma);
   return p;
 }
diff -u gsl-1.15/randist/fdist.c gsl-sam/randist/fdist.c
--- gsl-1.15/randist/fdist.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/fdist.c	2012-06-15 16:22:01.000000000 +0100
@@ -22,6 +22,7 @@
 #include <gsl/gsl_sf_gamma.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The F distribution has the form
 
@@ -46,9 +47,15 @@
 double
 gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2)
 {
+  return exp (gsl_ran_fdist_lnpdf (x, nu1, nu2));
+}
+
+double
+gsl_ran_fdist_lnpdf (const double x, const double nu1, const double nu2)
+{
   if (x < 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
@@ -59,9 +66,8 @@
       double lg1 = gsl_sf_lngamma (nu1 / 2);
       double lg2 = gsl_sf_lngamma (nu2 / 2);
 
-      p =
-	exp (lglg + lg12 - lg1 - lg2 + (nu1 / 2 - 1) * log (x) -
-	     ((nu1 + nu2) / 2) * log (nu2 + nu1 * x));
+      p = (lglg + lg12 - lg1 - lg2 + (nu1 / 2 - 1) * log (x) -
+	   ((nu1 + nu2) / 2) * log (nu2 + nu1 * x));
 
       return p;
     }
diff -u gsl-1.15/randist/flat.c gsl-sam/randist/flat.c
--- gsl-1.15/randist/flat.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/flat.c	2012-06-15 16:22:08.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* This is the uniform distribution in the range [a, b)
 
@@ -42,12 +43,18 @@
 double
 gsl_ran_flat_pdf (double x, const double a, const double b)
 {
+  return exp (gsl_ran_flat_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_flat_lnpdf (double x, const double a, const double b)
+{
   if (x < b && x >= a)
     {
-      return 1 / (b - a);
+      return -log (b - a);
     }
   else
     {
-      return 0;
+      return GSL_NEGINF;
     }
 }
diff -u gsl-1.15/randist/gamma.c gsl-sam/randist/gamma.c
--- gsl-1.15/randist/gamma.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gamma.c	2012-06-15 12:34:34.000000000 +0100
@@ -150,26 +150,32 @@
 double
 gsl_ran_gamma_pdf (const double x, const double a, const double b)
 {
+  return exp (gsl_ran_gamma_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_gamma_lnpdf (const double x, const double a, const double b)
+{
   if (x < 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else if (x == 0)
     {
       if (a == 1)
-        return 1/b ;
+        return -log(b) ;
       else
-        return 0 ;
+        return GSL_NEGINF ;
     }
   else if (a == 1)
     {
-      return exp(-x/b)/b ;
+      return (-x/b) - log(b) ;
     }
   else 
     {
       double p;
       double lngamma = gsl_sf_lngamma (a);
-      p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
+      p = (a - 1) * log (x/b) - x/b - lngamma - log (b);
       return p;
     }
 }
diff -u gsl-1.15/randist/gauss.c gsl-sam/randist/gauss.c
--- gsl-1.15/randist/gauss.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gauss.c	2012-06-15 13:31:48.000000000 +0100
@@ -118,8 +118,14 @@
 double
 gsl_ran_gaussian_pdf (const double x, const double sigma)
 {
+  return exp (gsl_ran_gaussian_lnpdf (x, sigma));
+}
+
+double
+gsl_ran_gaussian_lnpdf (const double x, const double sigma)
+{
   double u = x / fabs (sigma);
-  double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
+  double p = -log(sqrt (2 * M_PI) * fabs (sigma)) + (-u * u / 2);
   return p;
 }
 
@@ -136,8 +142,13 @@
 }
 
 double
-gsl_ran_ugaussian_pdf (const double x)
+gsl_ran_ugaussian_lnpdf (const double x)
 {
-  return gsl_ran_gaussian_pdf (x, 1.0);
+  return gsl_ran_gaussian_lnpdf (x, 1.0);
 }
 
+double
+gsl_ran_ugaussian_pdf (const double x)
+{
+  return exp (gsl_ran_gaussian_lnpdf (x, 1.0));
+}
diff -u gsl-1.15/randist/gausstail.c gsl-sam/randist/gausstail.c
--- gsl-1.15/randist/gausstail.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gausstail.c	2012-06-15 13:34:07.000000000 +0100
@@ -75,9 +75,15 @@
 double
 gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma)
 {
+  return exp (gsl_ran_gaussian_tail_lnpdf (x, a, sigma));
+}
+
+double
+gsl_ran_gaussian_tail_lnpdf (const double x, const double a, const double sigma)
+{
   if (x < a)
     {
-      return 0;
+      return GSL_NEGINF;
     }
   else
     {
@@ -88,7 +94,7 @@
 
       N = 0.5 * f;
 
-      p = (1 / (N * sqrt (2 * M_PI) * sigma)) * exp (-u * u / 2);
+      p = -log(N * sqrt (2 * M_PI) * sigma) + (-u * u / 2);
 
       return p;
     }
@@ -103,5 +109,11 @@
 double
 gsl_ran_ugaussian_tail_pdf (const double x, const double a)
 {
-  return gsl_ran_gaussian_tail_pdf (x, a, 1.0) ;
+  return exp (gsl_ran_gaussian_tail_lnpdf (x, a, 1.0)) ;
+}
+
+double
+gsl_ran_ugaussian_tail_lnpdf (const double x, const double a)
+{
+  return gsl_ran_gaussian_tail_lnpdf (x, a, 1.0) ;
 }
diff -u gsl-1.15/randist/geometric.c gsl-sam/randist/geometric.c
--- gsl-1.15/randist/geometric.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/geometric.c	2012-06-15 16:22:18.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* Geometric distribution (bernoulli trial with probability p) 
 
@@ -51,17 +52,23 @@
 double
 gsl_ran_geometric_pdf (const unsigned int k, const double p)
 {
+  return exp (gsl_ran_geometric_lnpdf (k, p));
+}
+
+double
+gsl_ran_geometric_lnpdf (const unsigned int k, const double p)
+{
   if (k == 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else if (k == 1)
     {
-      return p ;
+      return log (p) ;
     }
   else
     {
-      double P = p * pow (1 - p, k - 1.0);
+      double P = log (p) + log1p (-p) * (k - 1.0);
       return P;
     }
 }
diff -u gsl-1.15/randist/gsl_randist.h gsl-sam/randist/gsl_randist.h
--- gsl-1.15/randist/gsl_randist.h	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gsl_randist.h	2012-06-15 15:49:21.000000000 +0100
@@ -38,23 +38,29 @@
 
 double gsl_ran_beta (const gsl_rng * r, const double a, const double b);
 double gsl_ran_beta_pdf (const double x, const double a, const double b);
+double gsl_ran_beta_lnpdf (const double x, const double a, const double b);
 
 unsigned int gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n);
 unsigned int gsl_ran_binomial_knuth (const gsl_rng * r, double p, unsigned int n);
 unsigned int gsl_ran_binomial_tpe (const gsl_rng * r, double p, unsigned int n);
 double gsl_ran_binomial_pdf (const unsigned int k, const double p, const unsigned int n);
+double gsl_ran_binomial_lnpdf (const unsigned int k, const double p, const unsigned int n);
 
 double gsl_ran_exponential (const gsl_rng * r, const double mu);
 double gsl_ran_exponential_pdf (const double x, const double mu);
+double gsl_ran_exponential_lnpdf (const double x, const double mu);
 
 double gsl_ran_exppow (const gsl_rng * r, const double a, const double b);
 double gsl_ran_exppow_pdf (const double x, const double a, const double b);
+double gsl_ran_exppow_lnpdf (const double x, const double a, const double b);
 
 double gsl_ran_cauchy (const gsl_rng * r, const double a);
 double gsl_ran_cauchy_pdf (const double x, const double a);
+double gsl_ran_cauchy_lnpdf (const double x, const double a);
 
 double gsl_ran_chisq (const gsl_rng * r, const double nu);
 double gsl_ran_chisq_pdf (const double x, const double nu);
+double gsl_ran_chisq_lnpdf (const double x, const double nu);
 
 void gsl_ran_dirichlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
 double gsl_ran_dirichlet_pdf (const size_t K, const double alpha[], const double theta[]);
@@ -62,16 +68,20 @@
 
 double gsl_ran_erlang (const gsl_rng * r, const double a, const double n);
 double gsl_ran_erlang_pdf (const double x, const double a, const double n);
+double gsl_ran_erlang_lnpdf (const double x, const double a, const double n);
 
 double gsl_ran_fdist (const gsl_rng * r, const double nu1, const double nu2);
 double gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2);
+double gsl_ran_fdist_lnpdf (const double x, const double nu1, const double nu2);
 
 double gsl_ran_flat (const gsl_rng * r, const double a, const double b);
 double gsl_ran_flat_pdf (double x, const double a, const double b);
+double gsl_ran_flat_lnpdf (double x, const double a, const double b);
 
 double gsl_ran_gamma (const gsl_rng * r, const double a, const double b);
 double gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a);
 double gsl_ran_gamma_pdf (const double x, const double a, const double b);
+double gsl_ran_gamma_lnpdf (const double x, const double a, const double b);
 double gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b);
 double gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b);
 
@@ -79,43 +89,55 @@
 double gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma);
 double gsl_ran_gaussian_ziggurat (const gsl_rng * r, const double sigma);
 double gsl_ran_gaussian_pdf (const double x, const double sigma);
+double gsl_ran_gaussian_lnpdf (const double x, const double sigma);
 
 double gsl_ran_ugaussian (const gsl_rng * r);
 double gsl_ran_ugaussian_ratio_method (const gsl_rng * r);
 double gsl_ran_ugaussian_pdf (const double x);
+double gsl_ran_ugaussian_lnpdf (const double x);
 
 double gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma);
 double gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma);
+double gsl_ran_gaussian_tail_lnpdf (const double x, const double a, const double sigma);
 
 double gsl_ran_ugaussian_tail (const gsl_rng * r, const double a);
 double gsl_ran_ugaussian_tail_pdf (const double x, const double a);
+double gsl_ran_ugaussian_tail_lnpdf (const double x, const double a);
 
 void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y);
 double gsl_ran_bivariate_gaussian_pdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);
+double gsl_ran_bivariate_gaussian_lnpdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);
 
 double gsl_ran_landau (const gsl_rng * r);
 double gsl_ran_landau_pdf (const double x);
 
 unsigned int gsl_ran_geometric (const gsl_rng * r, const double p);
 double gsl_ran_geometric_pdf (const unsigned int k, const double p);
+double gsl_ran_geometric_lnpdf (const unsigned int k, const double p);
 
 unsigned int gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2, unsigned int t);
 double gsl_ran_hypergeometric_pdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);
+double gsl_ran_hypergeometric_lnpdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);
 
 double gsl_ran_gumbel1 (const gsl_rng * r, const double a, const double b);
 double gsl_ran_gumbel1_pdf (const double x, const double a, const double b);
+double gsl_ran_gumbel1_lnpdf (const double x, const double a, const double b);
 
 double gsl_ran_gumbel2 (const gsl_rng * r, const double a, const double b);
 double gsl_ran_gumbel2_pdf (const double x, const double a, const double b);
+double gsl_ran_gumbel2_lnpdf (const double x, const double a, const double b);
 
 double gsl_ran_logistic (const gsl_rng * r, const double a);
 double gsl_ran_logistic_pdf (const double x, const double a);
+double gsl_ran_logistic_lnpdf (const double x, const double a);
 
 double gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma);
 double gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma);
+double gsl_ran_lognormal_lnpdf (const double x, const double zeta, const double sigma);
 
 unsigned int gsl_ran_logarithmic (const gsl_rng * r, const double p);
 double gsl_ran_logarithmic_pdf (const unsigned int k, const double p);
+double gsl_ran_logarithmic_lnpdf (const unsigned int k, const double p);
 
 void gsl_ran_multinomial (const gsl_rng * r, const size_t K,
                           const unsigned int N, const double p[],
@@ -128,35 +150,44 @@
 
 unsigned int gsl_ran_negative_binomial (const gsl_rng * r, double p, double n);
 double gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n);
+double gsl_ran_negative_binomial_lnpdf (const unsigned int k, const double p, double n);
 
 unsigned int gsl_ran_pascal (const gsl_rng * r, double p, unsigned int n);
 double gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n);
+double gsl_ran_pascal_lnpdf (const unsigned int k, const double p, unsigned int n);
 
 double gsl_ran_pareto (const gsl_rng * r, double a, const double b);
 double gsl_ran_pareto_pdf (const double x, const double a, const double b);
+double gsl_ran_pareto_lnpdf (const double x, const double a, const double b);
 
 unsigned int gsl_ran_poisson (const gsl_rng * r, double mu);
 void gsl_ran_poisson_array (const gsl_rng * r, size_t n, unsigned int array[],
                             double mu);
 double gsl_ran_poisson_pdf (const unsigned int k, const double mu);
+double gsl_ran_poisson_lnpdf (const unsigned int k, const double mu);
 
 double gsl_ran_rayleigh (const gsl_rng * r, const double sigma);
 double gsl_ran_rayleigh_pdf (const double x, const double sigma);
+double gsl_ran_rayleigh_lnpdf (const double x, const double sigma);
 
 double gsl_ran_rayleigh_tail (const gsl_rng * r, const double a, const double sigma);
 double gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma);
+double gsl_ran_rayleigh_tail_lnpdf (const double x, const double a, const double sigma);
 
 double gsl_ran_tdist (const gsl_rng * r, const double nu);
 double gsl_ran_tdist_pdf (const double x, const double nu);
+double gsl_ran_tdist_lnpdf (const double x, const double nu);
 
 double gsl_ran_laplace (const gsl_rng * r, const double a);
 double gsl_ran_laplace_pdf (const double x, const double a);
+double gsl_ran_laplace_lnpdf (const double x, const double a);
 
 double gsl_ran_levy (const gsl_rng * r, const double c, const double alpha);
 double gsl_ran_levy_skew (const gsl_rng * r, const double c, const double alpha, const double beta);
 
 double gsl_ran_weibull (const gsl_rng * r, const double a, const double b);
 double gsl_ran_weibull_pdf (const double x, const double a, const double b);
+double gsl_ran_weibull_lnpdf (const double x, const double a, const double b);
 
 void gsl_ran_dir_2d (const gsl_rng * r, double * x, double * y);
 void gsl_ran_dir_2d_trig_method (const gsl_rng * r, double * x, double * y);
diff -u gsl-1.15/randist/gumbel.c gsl-sam/randist/gumbel.c
--- gsl-1.15/randist/gumbel.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/gumbel.c	2012-06-15 16:22:29.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The Type I Gumbel distribution has the form,
 
@@ -45,7 +46,13 @@
 double
 gsl_ran_gumbel1_pdf (const double x, const double a, const double b)
 {
-  double p = a * b *  exp (-(b * exp(-a * x) + a * x));
+  return exp (gsl_ran_gumbel1_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_gumbel1_lnpdf (const double x, const double a, const double b)
+{
+  double p = log(a) + log(b) - (b * exp(-a * x) + a * x);
   return p;
 }
 
@@ -62,13 +69,19 @@
 double
 gsl_ran_gumbel2_pdf (const double x, const double a, const double b)
 {
+  return exp (gsl_ran_gumbel2_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_gumbel2_lnpdf (const double x, const double a, const double b)
+{
   if (x <= 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
-      double p = b * a *  pow(x,-(a+1)) * exp (-b * pow(x, -a));
+      double p = log(b) + log(a) - log(x) * (a+1) - (b * pow(x, -a));
       return p;
     }
 }
diff -u gsl-1.15/randist/hyperg.c gsl-sam/randist/hyperg.c
--- gsl-1.15/randist/hyperg.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/hyperg.c	2012-06-15 16:22:39.000000000 +0100
@@ -22,6 +22,7 @@
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_math.h>
 
 /* The hypergeometric distribution has the form,
 
@@ -96,6 +97,15 @@
                             const unsigned int n2, 
                             unsigned int t)
 {
+  return exp (gsl_ran_hypergeometric_lnpdf (k, n1, n2, t));
+}
+
+double
+gsl_ran_hypergeometric_lnpdf (const unsigned int k, 
+			      const unsigned int n1, 
+			      const unsigned int n2, 
+			      unsigned int t)
+{
   if (t > n1 + n2)
     {
       t = n1 + n2 ;
@@ -103,11 +113,11 @@
 
   if (k > n1 || k > t)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else if (t > n2 && k + n2 < t )
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else 
     {
@@ -117,7 +127,7 @@
       double c2 = gsl_sf_lnchoose(n2,t-k);
       double c3 = gsl_sf_lnchoose(n1+n2,t);
 
-      p = exp(c1 + c2 - c3) ;
+      p = c1 + c2 - c3 ;
 
       return p;
     }
diff -u gsl-1.15/randist/laplace.c gsl-sam/randist/laplace.c
--- gsl-1.15/randist/laplace.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/laplace.c	2012-06-15 16:22:59.000000000 +0100
@@ -51,7 +51,14 @@
 double
 gsl_ran_laplace_pdf (const double x, const double a)
 {
-  double p = (1/(2*a)) * exp (-fabs (x)/a);
+  return exp (gsl_ran_laplace_lnpdf (x, a));
+}
+
+double
+gsl_ran_laplace_lnpdf (const double x, const double a)
+{
+  double p = -log(2*a) - (fabs (x)/a);
+
   return p;
 }
 
diff -u gsl-1.15/randist/logarithmic.c gsl-sam/randist/logarithmic.c
--- gsl-1.15/randist/logarithmic.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/logarithmic.c	2012-06-15 16:23:18.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* Logarithmic distribution 
 
@@ -64,13 +65,19 @@
 double
 gsl_ran_logarithmic_pdf (const unsigned int k, const double p)
 {
+  return exp (gsl_ran_logarithmic_lnpdf (k, p));
+}
+
+double
+gsl_ran_logarithmic_lnpdf (const unsigned int k, const double p)
+{
   if (k == 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else 
     {
-      double P = pow(p, (double)k) / (double) k / log(1/(1-p)) ;
+      double P = log(p) * k - log(k) - log(-log1p(-p)) ;
       return P;
     }
 }
diff -u gsl-1.15/randist/logistic.c gsl-sam/randist/logistic.c
--- gsl-1.15/randist/logistic.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/logistic.c	2012-06-15 14:20:03.000000000 +0100
@@ -47,7 +47,13 @@
 double
 gsl_ran_logistic_pdf (const double x, const double a)
 {
-  double u = exp (-fabs(x)/a);
-  double p = u / (fabs(a) * (1 + u) * (1 + u));
+  return exp (gsl_ran_logistic_lnpdf (x, a));
+}
+
+double
+gsl_ran_logistic_lnpdf (const double x, const double a)
+{
+  double u = -fabs(x)/a;
+  double p = u - (log(fabs(a)) + log1p (exp (u)) * 2);
   return p;
 }
diff -u gsl-1.15/randist/lognormal.c gsl-sam/randist/lognormal.c
--- gsl-1.15/randist/lognormal.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/lognormal.c	2012-06-15 16:14:53.000000000 +0100
@@ -57,14 +57,20 @@
 double
 gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma)
 {
+  return exp (gsl_ran_lognormal_lnpdf (x, zeta, sigma));
+}
+
+double
+gsl_ran_lognormal_lnpdf (const double x, const double zeta, const double sigma)
+{
   if (x <= 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
       double u = (log (x) - zeta)/sigma;
-      double p = 1 / (x * fabs(sigma) * sqrt (2 * M_PI)) * exp (-(u * u) /2);
+      double p = -log(x * fabs(sigma) * sqrt (2 * M_PI)) - (u * u) / 2;
       return p;
     }
 }
diff -u gsl-1.15/randist/nbinomial.c gsl-sam/randist/nbinomial.c
--- gsl-1.15/randist/nbinomial.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/nbinomial.c	2012-06-15 14:46:32.000000000 +0100
@@ -42,13 +42,19 @@
 double
 gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n)
 {
+  return exp( gsl_ran_negative_binomial_lnpdf (k, p, n));
+}
+
+double
+gsl_ran_negative_binomial_lnpdf (const unsigned int k, const double p, double n)
+{
   double P;
 
   double f = gsl_sf_lngamma (k + n) ;
   double a = gsl_sf_lngamma (n) ;
   double b = gsl_sf_lngamma (k + 1.0) ;
 
-  P = exp(f-a-b) * pow (p, n) * pow (1 - p, (double)k);
+  P = (f-a-b) + log (p) * n + log1p (-p) * k;
   
   return P;
 }
diff -u gsl-1.15/randist/pareto.c gsl-sam/randist/pareto.c
--- gsl-1.15/randist/pareto.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/pareto.c	2012-06-15 16:23:26.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The Pareto distribution has the form,
 
@@ -41,14 +42,21 @@
 double
 gsl_ran_pareto_pdf (const double x, const double a, const double b)
 {
+  return exp (gsl_ran_pareto_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_pareto_lnpdf (const double x, const double a, const double b)
+{
   if (x >= b)
     {
-      double p = (a/b) / pow (x/b, a + 1);
+      double lb = log(b);
+      double p = (log(a) - lb) - (log (x) - lb) * (a + 1);
       return p;
     }
   else
     {
-      return 0;
+      return GSL_NEGINF;
     }
 }
 
diff -u gsl-1.15/randist/pascal.c gsl-sam/randist/pascal.c
--- gsl-1.15/randist/pascal.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/pascal.c	2012-06-15 14:49:27.000000000 +0100
@@ -45,6 +45,11 @@
 double
 gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n)
 {
-  double P = gsl_ran_negative_binomial_pdf (k, p, (double) n);
-  return P;
+  return exp (gsl_ran_negative_binomial_lnpdf (k, p, (double) n));
+}
+
+double
+gsl_ran_pascal_lnpdf (const unsigned int k, const double p, unsigned int n)
+{
+  return gsl_ran_negative_binomial_lnpdf (k, p, (double) n);
 }
diff -u gsl-1.15/randist/poisson.c gsl-sam/randist/poisson.c
--- gsl-1.15/randist/poisson.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/poisson.c	2012-06-15 14:57:26.000000000 +0100
@@ -85,9 +85,15 @@
 double
 gsl_ran_poisson_pdf (const unsigned int k, const double mu)
 {
+  return exp(gsl_ran_poisson_lnpdf (k, mu));
+}
+
+double
+gsl_ran_poisson_lnpdf (const unsigned int k, const double mu)
+{
   double p;
   double lf = gsl_sf_lnfact (k); 
 
-  p = exp (log (mu) * k - lf - mu);
+  p = log (mu) * k - lf - mu;
   return p;
 }
diff -u gsl-1.15/randist/rayleigh.c gsl-sam/randist/rayleigh.c
--- gsl-1.15/randist/rayleigh.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/rayleigh.c	2012-06-15 16:23:34.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The Rayleigh distribution has the form
 
@@ -39,14 +40,20 @@
 double
 gsl_ran_rayleigh_pdf (const double x, const double sigma)
 {
+  return exp (gsl_ran_rayleigh_lnpdf (x, sigma));
+}
+
+double
+gsl_ran_rayleigh_lnpdf (const double x, const double sigma)
+{
   if (x < 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
       double u = x / sigma ;
-      double p = (u / sigma) * exp(-u * u / 2.0) ;
+      double p = log(u / sigma) - (u * u / 2.0) ;
       
       return p;
     }
@@ -69,17 +76,23 @@
 double
 gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma)
 {
+  return exp (gsl_ran_rayleigh_tail_lnpdf (x, a, sigma));
+}
+
+double
+gsl_ran_rayleigh_tail_lnpdf (const double x, const double a, const double sigma)
+{
   if (x < a)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else
     {
       double u = x / sigma ;
       double v = a / sigma ;
 
-      double p = (u / sigma) * exp((v + u) * (v - u) / 2.0) ;
-      
+      double p = log(u / sigma) + ((v + u) * (v - u) / 2.0) ;
+
       return p;
     }
 }
diff -u gsl-1.15/randist/tdist.c gsl-sam/randist/tdist.c
--- gsl-1.15/randist/tdist.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/tdist.c	2012-06-15 17:18:00.000000000 +0100
@@ -67,13 +67,20 @@
 double
 gsl_ran_tdist_pdf (const double x, const double nu)
 {
+  return exp (gsl_ran_tdist_lnpdf (x, nu));
+}
+
+double
+gsl_ran_tdist_lnpdf (const double x, const double nu)
+{
   double p;
 
   double lg1 = gsl_sf_lngamma (nu / 2);
   double lg2 = gsl_sf_lngamma ((nu + 1) / 2);
+  
+  p = (((lg2 - lg1) - log(sqrt (M_PI * nu)))
+       + log1p (x * x / nu) * ((-1 - nu) / 2));
 
-  p = ((exp (lg2 - lg1) / sqrt (M_PI * nu)) 
-       * pow ((1 + x * x / nu), -(nu + 1) / 2));
   return p;
 }
 
diff -u gsl-1.15/randist/weibull.c gsl-sam/randist/weibull.c
--- gsl-1.15/randist/weibull.c	2010-12-26 17:57:08.000000000 +0000
+++ gsl-sam/randist/weibull.c	2012-06-15 16:23:46.000000000 +0100
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_math.h>
 
 /* The Weibull distribution has the form,
 
@@ -41,24 +42,30 @@
 double
 gsl_ran_weibull_pdf (const double x, const double a, const double b)
 {
+  return exp (gsl_ran_weibull_lnpdf (x, a, b));
+}
+
+double
+gsl_ran_weibull_lnpdf (const double x, const double a, const double b)
+{
   if (x < 0)
     {
-      return 0 ;
+      return GSL_NEGINF ;
     }
   else if (x == 0)
     {
       if (b == 1)
-        return 1/a ;
+        return -log(a) ;
       else
-        return 0 ;
+        return GSL_NEGINF ;
     }
   else if (b == 1)
     {
-      return exp(-x/a)/a ;
+      return (-x/a) - log(a) ;
     }
   else
     {
-      double p = (b/a) * exp (-pow (x/a, b) + (b - 1) * log (x/a));
+      double p = log(b/a) + (-pow (x/a, b) + (b - 1) * log (x/a));
       return p;
     }
 }