Index: doc/guide/plugin.rst
===================================================================
--- doc/guide/plugin.rst	(revision 2015f021ed00f8b3e48aa9ae7f690780441be4d8)
+++ doc/guide/plugin.rst	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -428,14 +428,14 @@
 	def random():
 	...
-	
-This function provides a model-specific random parameter set which shows model 
-features in the USANS to SANS range.  For example, core-shell sphere sets the 
-outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q 
-value for the transition from flat to falling.  It then uses a beta distribution 
-to set the percentage of the shape which is shell, giving a preference for very 
-thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp 
-should show a reasonable variety of curves over the default sascomp q range.  
-The parameter set is returned as a dictionary of `{parameter: value, ...}`.  
-Any model parameters not included in the dictionary will default according to 
+
+This function provides a model-specific random parameter set which shows model
+features in the USANS to SANS range.  For example, core-shell sphere sets the
+outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q
+value for the transition from flat to falling.  It then uses a beta distribution
+to set the percentage of the shape which is shell, giving a preference for very
+thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp
+should show a reasonable variety of curves over the default sascomp q range.
+The parameter set is returned as a dictionary of `{parameter: value, ...}`.
+Any model parameters not included in the dictionary will default according to
 the code in the `_randomize_one()` function from sasmodels/compare.py.
 
@@ -701,6 +701,6 @@
     erf, erfc, tgamma, lgamma:  **do not use**
         Special functions that should be part of the standard, but are missing
-        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma
-        instead (see below). Note: lgamma(x) has not yet been tested.
+        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma
+        and sas_lgamma instead (see below).
 
 Some non-standard constants and functions are also provided:
@@ -769,9 +769,27 @@
         Gamma function sas_gamma\ $(x) = \Gamma(x)$.
 
-        The standard math function, tgamma(x) is unstable for $x < 1$
+        The standard math function, tgamma(x), is unstable for $x < 1$
         on some platforms.
 
         :code:`source = ["lib/sas_gamma.c", ...]`
         (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
+
+    sas_gammaln(x):
+        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
+
+        The standard math function, lgamma(x), is incorrect for single
+        precision on some platforms.
+
+        :code:`source = ["lib/sas_gammainc.c", ...]`
+        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
+
+    sas_gammainc(a, x), sas_gammaincc(a, x):
+        Incomplete gamma function
+        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
+        and complementary incomplete gamma function
+        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
+
+        :code:`source = ["lib/sas_gammainc.c", ...]`
+        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
 
     sas_erf(x), sas_erfc(x):
Index: explore/precision.py
===================================================================
--- explore/precision.py	(revision 2a7e20e9c4ba72da5c6088e743a64f3345ee554d)
+++ explore/precision.py	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -95,4 +95,7 @@
             neg:    [-100,100]
 
+        For arbitrary range use "start:stop:steps:scale" where scale is
+        one of log, lin, or linear.
+
         *diff* is "relative", "absolute" or "none"
 
@@ -102,14 +105,21 @@
         linear = not xrange.startswith("log")
         if xrange == "zoom":
-            lin_min, lin_max, lin_steps = 1000, 1010, 2000
+            start, stop, steps = 1000, 1010, 2000
         elif xrange == "neg":
-            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000
+            start, stop, steps = -100.1, 100.1, 2000
         elif xrange == "linear":
-            lin_min, lin_max, lin_steps = 1, 1000, 2000
-            lin_min, lin_max, lin_steps = 0.001, 2, 2000
+            start, stop, steps = 1, 1000, 2000
+            start, stop, steps = 0.001, 2, 2000
         elif xrange == "log":
-            log_min, log_max, log_steps = -3, 5, 400
+            start, stop, steps = -3, 5, 400
         elif xrange == "logq":
-            log_min, log_max, log_steps = -4, 1, 400
+            start, stop, steps = -4, 1, 400
+        elif ':' in xrange:
+            parts = xrange.split(':')
+            linear = parts[3] != "log" if len(parts) == 4 else True
+            steps = int(parts[2]) if len(parts) > 2 else 400
+            start = float(parts[0])
+            stop = float(parts[1])
+
         else:
             raise ValueError("unknown range "+xrange)
@@ -121,17 +131,17 @@
             # value to x in the given precision.
             if linear:
-                lin_min = max(lin_min, self.limits[0])
-                lin_max = min(lin_max, self.limits[1])
-                qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single')
-                #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double')
+                start = max(start, self.limits[0])
+                stop = min(stop, self.limits[1])
+                qrf = np.linspace(start, stop, steps, dtype='single')
+                #qrf = np.linspace(start, stop, steps, dtype='double')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = mp.linspace(lin_min, lin_max, lin_steps)
+                #qr = mp.linspace(start, stop, steps)
             else:
-                log_min = np.log10(max(10**log_min, self.limits[0]))
-                log_max = np.log10(min(10**log_max, self.limits[1]))
-                qrf = np.logspace(log_min, log_max, log_steps, dtype='single')
-                #qrf = np.logspace(log_min, log_max, log_steps, dtype='double')
+                start = np.log10(max(10**start, self.limits[0]))
+                stop = np.log10(min(10**stop, self.limits[1]))
+                qrf = np.logspace(start, stop, steps, dtype='single')
+                #qrf = np.logspace(start, stop, steps, dtype='double')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
+                #qr = [10**v for v in mp.linspace(start, stop, steps)]
 
         target = self.call_mpmath(qr, bits=500)
@@ -176,5 +186,5 @@
     """
     if diff == "relative":
-        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
+        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
         pylab.loglog(x, err, '-', label=label)
@@ -197,4 +207,21 @@
     return model_info
 
+# Hack to allow second parameter A in two parameter functions
+A = 1
+def parse_extra_pars():
+    global A
+
+    A_str = str(A)
+    pop = []
+    for k, v in enumerate(sys.argv[1:]):
+        if v.startswith("A="):
+            A_str = v[2:]
+            pop.append(k+1)
+    if pop:
+        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop]
+        A = float(A_str)
+
+parse_extra_pars()
+
 
 # =============== FUNCTION DEFINITIONS ================
@@ -297,4 +324,23 @@
     ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
     limits=(-3.1, 10),
+)
+add_function(
+    name="gammaln(x)",
+    mp_function=mp.loggamma,
+    np_function=scipy.special.gammaln,
+    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]),
+    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"),
+)
+add_function(
+    name="gammainc(x)",
+    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a),
+    np_function=lambda x, a=A: scipy.special.gammainc(a, x),
+    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]),
+)
+add_function(
+    name="gammaincc(x)",
+    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a),
+    np_function=lambda x, a=A: scipy.special.gammaincc(a, x),
+    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]),
 )
 add_function(
@@ -463,6 +509,6 @@
 lanczos_gamma = """\
     const double coeff[] = {
-            76.18009172947146,     -86.50532032941677,
-            24.01409824083091,     -1.231739572450155,
+            76.18009172947146, -86.50532032941677,
+            24.01409824083091, -1.231739572450155,
             0.1208650973866179e-2,-0.5395239384953e-5
             };
@@ -475,5 +521,5 @@
 """
 add_function(
-    name="log gamma(x)",
+    name="loggamma(x)",
     mp_function=mp.loggamma,
     np_function=scipy.special.gammaln,
@@ -599,5 +645,5 @@
 
 ALL_FUNCTIONS = set(FUNCTIONS.keys())
-ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
+ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead
 ALL_FUNCTIONS.discard("3j1/x:taylor")
 ALL_FUNCTIONS.discard("3j1/x:trig")
@@ -615,10 +661,15 @@
     -r indicates that the relative error should be plotted (default),
     -x<range> indicates the steps in x, where <range> is one of the following
-      log indicates log stepping in [10^-3, 10^5] (default)
-      logq indicates log stepping in [10^-4, 10^1]
-      linear indicates linear stepping in [1, 1000]
-      zoom indicates linear stepping in [1000, 1010]
-      neg indicates linear stepping in [-100.1, 100.1]
-and name is "all" or one of:
+        log indicates log stepping in [10^-3, 10^5] (default)
+        logq indicates log stepping in [10^-4, 10^1]
+        linear indicates linear stepping in [1, 1000]
+        zoom indicates linear stepping in [1000, 1010]
+        neg indicates linear stepping in [-100.1, 100.1]
+        start:stop:n[:stepping] indicates an n-step plot in [start, stop]
+            or [10^start, 10^stop] if stepping is "log" (default n=400)
+Some functions (notably gammainc/gammaincc) have an additional parameter A
+which can be set from the command line as A=value.  Default is A=1.
+
+Name is one of:
     """+names)
     sys.exit(1)
Index: sasmodels/models/lib/sas_gammainc.c
===================================================================
--- sasmodels/models/lib/sas_gammainc.c	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
+++ sasmodels/models/lib/sas_gammainc.c	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -0,0 +1,458 @@
+#if FLOAT_SIZE>4  // double precision
+// based on cephes/double/igam.c
+double gammaln(double x);
+double gammainc(double a, double x);
+double gammaincc(double a, double x);
+
+double cephes_igamc(double a, double x);
+double cephes_igam(double a, double x);
+double cephes_lgam(double x);
+double cephes_lgam2(double x);
+
+double gammainc(double a, double x)
+{
+    if ((x <= 0) || ( a <= 0)) return 0.0;
+    if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamc(a, x);
+    return cephes_igam(a, x);
+}
+
+double gammaincc(double a, double x)
+{
+    if ((x <= 0) || (a <= 0)) return 1.0;
+    if ((x < 1.0) || (x < a)) return 1.0 - cephes_igam(a, x);
+    return cephes_igamc(a, x);
+}
+
+double gammaln(double x)
+{
+    if (isnan(x)) return(x);
+    if (!isfinite(x)) return(INFINITY);
+    return cephes_lgam(x);
+}
+
+double cephes_igamc(double a, double x)
+{
+    const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53
+    const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized
+    const double BIG = 4.503599627370496e15;
+    const double BIGINV = 2.22044604925031308085e-16;
+    double ans, ax, c, yc, r, t, y, z;
+    double pk, pkm1, pkm2, qk, qkm1, qkm2;
+
+    /* Compute  x**a * exp(-x) / gamma(a)  */
+    ax = a * log(x) - x - cephes_lgam(a);
+    if (ax < -MAXLOG) return 0.0;  // underflow
+    ax = exp(ax);
+
+    /* continued fraction */
+    y = 1.0 - a;
+    z = x + y + 1.0;
+    c = 0.0;
+    pkm2 = 1.0;
+    qkm2 = x;
+    pkm1 = x + 1.0;
+    qkm1 = z * x;
+    ans = pkm1/qkm1;
+
+    do {
+        c += 1.0;
+        y += 1.0;
+        z += 2.0;
+        yc = y * c;
+        pk = pkm1 * z  -  pkm2 * yc;
+        qk = qkm1 * z  -  qkm2 * yc;
+        if (qk != 0) {
+            r = pk/qk;
+            t = fabs( (ans - r)/r );
+            ans = r;
+        } else {
+            t = 1.0;
+        }
+        pkm2 = pkm1;
+        pkm1 = pk;
+        qkm2 = qkm1;
+        qkm1 = qk;
+        if (fabs(pk) > BIG) {
+            pkm2 *= BIGINV;
+            pkm1 *= BIGINV;
+            qkm2 *= BIGINV;
+            qkm1 *= BIGINV;
+        }
+    } while( t > MACHEP );
+
+    return( ans * ax );
+}
+
+
+double cephes_igam(double a, double x)
+{
+    const double MACHEP = 1.11022302462515654042E-16; // IEEE 2**-53
+    const double MAXLOG = 7.09782712893383996843E2; // IEEE log(2**1024) denormalized
+    double ans, ax, c, r;
+
+    /* Compute  x**a * exp(-x) / gamma(a)  */
+    ax = a * log(x) - x - cephes_lgam(a);
+    if (ax < -MAXLOG) return 0.0; // underflow
+    ax = exp(ax);
+
+    /* power series */
+    r = a;
+    c = 1.0;
+    ans = 1.0;
+
+    do {
+        r += 1.0;
+        c *= x/r;
+        ans += c;
+    } while (c/ans > MACHEP);
+
+    return ans * ax/a;
+}
+
+/* Logarithm of gamma function */
+
+
+double cephes_lgam(double x)
+{
+    const double LOGPI = 1.14472988584940017414; // log(pi)
+    int sgngam = 1;
+    double p, q, w, z;
+    int i;
+
+    if (x < -34.0) {
+        q = -x;
+        w = cephes_lgam2(q);
+        p = floor(q);
+        if (p == q) return INFINITY;
+        i = p;
+        sgngam = ((i&1) == 0) ? -1 : 1;
+        z = q - p;
+        if (z > 0.5) {
+            p += 1.0;
+            z = p - q;
+        }
+        z = q * sin(M_PI * z);
+        if (z == 0.0) return INFINITY;
+        z = LOGPI - log(z) - w;
+        return z;
+    } else {
+        return cephes_lgam2(x);
+    }
+}
+
+double cephes_lgam2(double x)
+{
+    const double LS2PI = 0.91893853320467274178; // log(sqrt(2*pi))
+    const double MAXLGM = 2.556348e305;
+    int sgngam = 1;
+    double p, q, u, z;
+    double A, B, C;
+
+    if (x < 13.0) {
+        z = 1.0;
+        p = 0.0;
+        u = x;
+        while (u >= 3.0) {
+            p -= 1.0;
+            u = x + p;
+            z *= u;
+        }
+        while (u < 2.0) {
+            if (u == 0.0) return INFINITY;
+            z /= u;
+            p += 1.0;
+            u = x + p;
+        }
+        if (z < 0.0) {
+            sgngam = -1;
+            z = -z;
+        } else {
+            sgngam = 1;
+        }
+        if (u == 2.0) return log(z);
+        p -= 2.0;
+        x = x + p;
+        B = ((((((
+            -1.37825152569120859100E3)*x
+            -3.88016315134637840924E4)*x
+            -3.31612992738871184744E5)*x
+            -1.16237097492762307383E6)*x
+            -1.72173700820839662146E6)*x
+            -8.53555664245765465627E5);
+        C = ((((((
+            /* 1.00000000000000000000E0)* */x
+            -3.51815701436523470549E2)*x
+            -1.70642106651881159223E4)*x
+            -2.20528590553854454839E5)*x
+            -1.13933444367982507207E6)*x
+            -2.53252307177582951285E6)*x
+            -2.01889141433532773231E6);
+        p = x * B / C;
+        return log(z) + p;
+    }
+
+    if (x > MAXLGM) return sgngam * INFINITY;
+
+    q = (x - 0.5) * log(x) - x + LS2PI;
+    if (x > 1.0e8) return q;
+
+    p = 1.0/(x*x);
+    if (x >= 1000.0) {
+        q += ((7.9365079365079365079365e-4 * p
+            - 2.7777777777777777777778e-3) * p
+            + 0.0833333333333333333333) / x;
+    } else {
+        A = (((((
+            +8.11614167470508450300E-4)*p
+            -5.95061904284301438324E-4)*p
+            +7.93650340457716943945E-4)*p
+            -2.77777777730099687205E-3)*p
+            +8.33333333333331927722E-2);
+        q += A / x;
+    }
+    return q;
+}
+
+#else // single precision
+
+// based on cephes/float/igam.c
+float gammalnf(float x);
+float gammaincf(float a, float x);
+float gammainccf(float a, float x);
+
+float cephes_igamcf(float a, float x);
+float cephes_igamf(float a, float x);
+float cephes_lgamf(float x);
+float cephes_lgam2f(float x);
+
+// Note: original uses logf, fabsf, floorf and sinf
+
+float gammaincf(float a, float x)
+{
+    if ((x <= 0) || ( a <= 0)) return 0.0;
+    if ((x > 1.0) && (x > a)) return 1.0 - cephes_igamcf(a, x);
+    return cephes_igamf(a, x);
+}
+
+float gammainccf(float a, float x)
+{
+    if ((x <= 0) || (a <= 0)) return 1.0;
+    if ((x < 1.0) || (x < a)) return 1.0 - cephes_igamf(a, x);
+    return cephes_igamcf(a, x);
+}
+
+float gammalnf(float x)
+{
+    if (isnan(x)) return(x);
+    if (!isfinite(x)) return(INFINITY);
+    return cephes_lgamf(x);
+}
+
+float cephes_igamcf(float a, float x)
+{
+    const float MAXLOGF = 88.72283905206835;
+    const float MACHEPF = 5.9604644775390625E-8;
+    const float BIG = 16777216.;
+    float ans, c, yc, ax, y, z;
+    float pk, pkm1, pkm2, qk, qkm1, qkm2;
+    float r, t;
+
+    ax = a * log(x) - x - cephes_lgamf(a);
+    if (ax < -MAXLOGF) return 0.0; // underflow
+    ax = expf(ax);
+
+    /* continued fraction */
+    y = 1.0 - a;
+    z = x + y + 1.0;
+    c = 0.0;
+    pkm2 = 1.0;
+    qkm2 = x;
+    pkm1 = x + 1.0;
+    qkm1 = z * x;
+    ans = pkm1/qkm1;
+
+    do {
+        c += 1.0;
+        y += 1.0;
+        z += 2.0;
+        yc = y * c;
+        pk = pkm1 * z  -  pkm2 * yc;
+        qk = qkm1 * z  -  qkm2 * yc;
+        if (qk != 0) {
+            r = pk/qk;
+            t = fabs((ans - r)/r);
+            ans = r;
+        } else {
+            t = 1.0;
+        }
+        pkm2 = pkm1;
+        pkm1 = pk;
+        qkm2 = qkm1;
+        qkm1 = qk;
+        if (fabs(pk) > BIG) {
+            pkm2 *= MACHEPF;
+            pkm1 *= MACHEPF;
+            qkm2 *= MACHEPF;
+            qkm1 *= MACHEPF;
+        }
+    } while (t > MACHEPF);
+
+    return ans * ax;
+}
+
+float cephes_igamf(float a, float x)
+{
+    const float MAXLOGF = 88.72283905206835;
+    const float MACHEPF = 5.9604644775390625E-8;
+    float ans, ax, c, r;
+
+    /* Compute  x**a * exp(-x) / gamma(a)  */
+    ax = a * log(x) - x - cephes_lgamf(a);
+    if (ax < -MAXLOGF) return 0.0; // underflow
+    ax = expf(ax);
+
+    /* power series */
+    r = a;
+    c = 1.0;
+    ans = 1.0;
+
+    do {
+        r += 1.0;
+        c *= x/r;
+        ans += c;
+    } while (c/ans > MACHEPF);
+
+    return ans * ax/a;
+}
+
+float cephes_lgamf(float x)
+{
+    const float PIINV =  0.318309886183790671538;
+    float p, q, w, z;
+    int i;
+    int sgngamf;
+
+    if (x < 0.0) {
+        q = -x;
+        w = cephes_lgam2f(q);
+        p = floor(q);
+        if (p == q) return INFINITY;
+        i = p;
+        sgngamf = ((i&1) == 0) ? -1 : 1;
+        z = q - p;
+        if (z > 0.5) {
+            p += 1.0;
+            z = p - q;
+        }
+        z = q * sin(M_PI * z);
+        if (z == 0.0) return sgngamf * INFINITY;
+        z = -log(PIINV*z) - w;
+        return z;
+    } else {
+        return cephes_lgam2f(x);
+    }
+}
+
+float cephes_lgam2f(float x)
+{
+    const float LS2PI = 0.91893853320467274178; // log(sqrt(2*pi))
+    const float MAXLGM = 2.035093e36;
+    float p, q, z;
+    float nx, tx;
+    int direction;
+    int sgngamf = 1;
+
+    if (x < 6.5) {
+        direction = 0;
+        z = 1.0;
+        tx = x;
+        nx = 0.0;
+        if (x >= 1.5) {
+            while (tx > 2.5) {
+                nx -= 1.0;
+                tx = x + nx;
+                z *=tx;
+            }
+            x += nx - 2.0;
+        } else if (x >= 1.25) {
+            z *= x;
+            x -= 1.0; /* x + 1 - 2 */
+            direction = 1;
+        } else if (x >= 0.75) {
+            x -= 1.0;
+            /* log gamma(x+1), -.25 < x < .25 */
+            // p = x * polevlf( x, C, 7 );
+            p = ((((((((
+                +1.369488127325832E-001)*x
+                -1.590086327657347E-001)*x
+                +1.692415923504637E-001)*x
+                -2.067882815621965E-001)*x
+                +2.705806208275915E-001)*x
+                -4.006931650563372E-001)*x
+                +8.224670749082976E-001)*x
+                -5.772156501719101E-001)*x;
+            q = 0.0;
+            return p + q;
+        } else {
+            while (tx < 1.5) {
+                if (tx == 0.0) return sgngamf * INFINITY;
+                z *=tx;
+                nx += 1.0;
+                tx = x + nx;
+            }
+            direction = 1;
+            x += nx - 2.0;
+        }
+
+        /* log gamma(x+2), -.5 < x < .5 */
+        // p = x * polevlf( x, B, 7 );
+        p = ((((((((
+            +6.055172732649237E-004)*x
+            -1.311620815545743E-003)*x
+            +2.863437556468661E-003)*x
+            -7.366775108654962E-003)*x
+            +2.058355474821512E-002)*x
+            -6.735323259371034E-002)*x
+            +3.224669577325661E-001)*x
+            +4.227843421859038E-001)*x;
+
+        if (z < 0.0) {
+            sgngamf = -1;
+            z = -z;
+        } else {
+            sgngamf = 1;
+        }
+        q = log(z);
+        if (direction) q = -q;
+        return p + q;
+    } else if (x > MAXLGM) {
+        return sgngamf * INFINITY;
+    } else {
+        /* Note, though an asymptotic formula could be used for x >= 3,
+        * there is cancellation error in the following if x < 6.5.
+        */
+        q = LS2PI - x;
+        q += (x - 0.5) * log(x);
+
+        if (x <= 1.0e4) {
+            z = 1.0/x;
+            p = z * z;
+            q += ((6.789774945028216E-004 * p
+                - 2.769887652139868E-003 ) * p
+                +  8.333316229807355E-002 ) * z;
+        }
+        return q;
+    }
+}
+
+#endif // !single precision
+
+#if FLOAT_SIZE>4
+#define sas_gammaln gammaln
+#define sas_gammainc gammainc
+#define sas_gammaincc gammaincc
+#else
+#define sas_gammaln gammalnf
+#define sas_gammainc gammaincf
+#define sas_gammaincc gammainccf
+#endif
Index: sasmodels/special.py
===================================================================
--- sasmodels/special.py	(revision df69efa2cdb726d21ce081492d47b44552fc0d19)
+++ sasmodels/special.py	(revision fba9ca0239129ed28284a0b435ba2a9524913f36)
@@ -113,4 +113,16 @@
         The standard math function, tgamma(x) is unstable for $x < 1$
         on some platforms.
+
+    sas_gammaln(x):
+        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
+
+        The standard math function, lgamma(x), is incorrect for single
+        precision on some platforms.
+
+    sas_gammainc(a, x), sas_gammaincc(a, x):
+        Incomplete gamma function
+        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
+        and complementary incomplete gamma function
+        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
 
     sas_erf(x), sas_erfc(x):
@@ -207,4 +219,7 @@
 from numpy import pi, nan, inf
 from scipy.special import gamma as sas_gamma
+from scipy.special import gammaln as sas_gammaln
+from scipy.special import gammainc as sas_gammainc
+from scipy.special import gammaincc as sas_gammaincc
 from scipy.special import erf as sas_erf
 from scipy.special import erfc as sas_erfc
