Changeset a062d18 in sasmodels

Timestamp:

Apr 27, 2016 7:27:02 AM (8 years ago)

Author:

Piotr Rozyczko <piotr.rozyczko@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

29464ba

Parents:

2a55a6f

git-author:

Piotr Rozyczko <piotr.rozyczko@…> (04/27/16 07:18:43)

git-committer:

Piotr Rozyczko <piotr.rozyczko@…> (04/27/16 07:27:02)

Message:

Reverted the workaround - needs more analysis.

File:

• sasmodels/models/lib/sas_gamma.c (modified) (1 diff)

sasmodels/models/lib/sas_gamma.c

@@ -7 +7 @@
 */
 
-#if defined(NEED_TGAMMA)
-static double cephes_stirf(double x)
-{
-        const double MAXSTIR=143.01608;
-        const double SQTPI=2.50662827463100050242E0;
-        double y, w, v;
-
-        w = 1.0 / x;
-
-        w = ((((
-                7.87311395793093628397E-4*w
-                - 2.29549961613378126380E-4)*w
-                - 2.68132617805781232825E-3)*w
-                + 3.47222221605458667310E-3)*w
-                + 8.33333333333482257126E-2)*w
-                + 1.0;
-        y = exp(x);
-        if (x > MAXSTIR)
-        { /* Avoid overflow in pow() */
-                v = pow(x, 0.5 * x - 0.25);
-                y = v * (v / y);
-        }
-        else
-        {
-                y = pow(x, x - 0.5) / y;
-        }
-        y = SQTPI * y * w;
-        return(y);
-}
-
-static double tgamma(double x) {
-        double p, q, z;
-        int sgngam;
-        int i;
-
-        sgngam = 1;
-        if (isnan(x))
-                return(x);
-        q = fabs(x);
-
-        if (q > 33.0)
-        {
-                if (x < 0.0)
-                {
-                        p = floor(q);
-                        if (p == q)
-                        {
-                                return (NAN);
-                        }
-                        i = p;
-                        if ((i & 1) == 0)
-                                sgngam = -1;
-                        z = q - p;
-                        if (z > 0.5)
-                        {
-                                p += 1.0;
-                                z = q - p;
-                        }
-                        z = q * sin(M_PI * z);
-                        if (z == 0.0)
-                        {
-                                return(NAN);
-                        }
-                        z = fabs(z);
-                        z = M_PI / (z * cephes_stirf(q));
-                }
-                else
-                {
-                        z = cephes_stirf(x);
-                }
-                return(sgngam * z);
-        }
-
-        z = 1.0;
-        while (x >= 3.0)
-        {
-                x -= 1.0;
-                z *= x;
-        }
-
-        while (x < 0.0)
-        {
-                if (x > -1.E-9)
-                        goto small;
-                z /= x;
-                x += 1.0;
-        }
-
-        while (x < 2.0)
-        {
-                if (x < 1.e-9)
-                        goto small;
-                z /= x;
-                x += 1.0;
-        }
-
-        if (x == 2.0)
-                return(z);
-
-        x -= 2.0;
-        p = (((((
-                +1.60119522476751861407E-4*x
-                + 1.19135147006586384913E-3)*x
-                + 1.04213797561761569935E-2)*x
-                + 4.76367800457137231464E-2)*x
-                + 2.07448227648435975150E-1)*x
-                + 4.94214826801497100753E-1)*x
-                + 9.99999999999999996796E-1;
-        q = ((((((
-                -2.31581873324120129819E-5*x
-                + 5.39605580493303397842E-4)*x
-                - 4.45641913851797240494E-3)*x
-                + 1.18139785222060435552E-2)*x
-                + 3.58236398605498653373E-2)*x
-                - 2.34591795718243348568E-1)*x
-                + 7.14304917030273074085E-2)*x
-                + 1.00000000000000000320E0;
-        return(z * p / q);
-
-small:
-        if (x == 0.0)
-        {
-                return (NAN);
-        }
-        else
-                return(z / ((1.0 + 0.5772156649015329 * x) * x));
-}
-#endif
-
-
 inline double sas_gamma( double x) { return tgamma(x+1)/x; }