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
