--- a/liboctave/numeric/randgamma.cc	Sat Mar 17 07:02:30 2018 -0500
+++ b/liboctave/numeric/randgamma.cc	Sat Mar 17 11:59:05 2018 -0500
@@ -89,37 +89,31 @@
 
 namespace octave
 {
-
-#define INFINITE lo_ieee_isinf
-#define RUNI rand_uniform<double> ()
-#define RNOR rand_normal<double> ()
-#define REXP rand_exponential<double> ()
-
-  template <> void rand_gamma<double> (double a, octave_idx_type n, double *r)
+  template <typename T> void rand_gamma (T a, octave_idx_type n, T *r)
   {
     octave_idx_type i;
     /* If a < 1, start by generating gamma (1+a) */
-    const double d = (a < 1. ? 1.+a : a) - 1./3.;
-    const double c = 1./std::sqrt (9.*d);
+    const T d = (a < 1. ? 1.+a : a) - 1./3.;
+    const T c = 1./std::sqrt (9.*d);
 
     /* Handle invalid cases */
-    if (a <= 0 || INFINITE(a))
+    if (a <= 0 || lo_ieee_isinf (a))
       {
         for (i=0; i < n; i++)
-          r[i] = numeric_limits<double>::NaN ();
+          r[i] = numeric_limits<T>::NaN ();
         return;
       }
 
     for (i=0; i < n; i++)
       {
-        double x, xsq, v, u;
+        T x, xsq, v, u;
       restart:
-        x = RNOR;
+        x = rand_normal<T> ();
         v = (1+c*x);
         v *= v*v;
         if (v <= 0)
           goto restart; /* rare, so don't bother moving up */
-        u = RUNI;
+        u = rand_uniform<T> ();
         xsq = x*x;
         if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v)))
           goto restart;
@@ -130,55 +124,10 @@
         /* Use gamma(a) = gamma(1+a)*U^(1/a) */
         /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
         for (i = 0; i < n; i++)
-          r[i] *= exp (-REXP/a);
+          r[i] *= exp (-rand_exponential<T> () / a);
       }
   }
 
-#undef RUNI
-#undef RNOR
-#undef REXP
-
-#define RUNI rand_uniform<float> ()
-#define RNOR rand_normal<float> ()
-#define REXP rand_exponential<float> ()
-
-  template <> void rand_gamma<float> (float a, octave_idx_type n, float *r)
-  {
-    octave_idx_type i;
-    /* If a < 1, start by generating gamma(1+a) */
-    const float d = (a < 1. ? 1.+a : a) - 1./3.;
-    const float c = 1./std::sqrt (9.*d);
-
-    /* Handle invalid cases */
-    if (a <= 0 || INFINITE(a))
-      {
-        for (i=0; i < n; i++)
-          r[i] = numeric_limits<float>::NaN ();
-        return;
-      }
-
-    for (i=0; i < n; i++)
-      {
-        float x, xsq, v, u;
-      frestart:
-        x = RNOR;
-        v = (1+c*x);
-        v *= v*v;
-        if (v <= 0)
-          goto frestart; /* rare, so don't bother moving up */
-        u = RUNI;
-        xsq = x*x;
-        if (u >= 1.-0.0331*xsq*xsq && std::log (u) >= 0.5*xsq + d*(1-v+std::log (v)))
-          goto frestart;
-        r[i] = d*v;
-      }
-    if (a < 1)
-      {
-        /* Use gamma(a) = gamma(1+a)*U^(1/a) */
-        /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
-        for (i = 0; i < n; i++)
-          r[i] *= exp (-REXP/a);
-      }
-  }
-
+  template void rand_gamma (double, octave_idx_type, double *);
+  template void rand_gamma (float, octave_idx_type, float *);
 }