1.1 --- a/include/vigra/mathutil.hxx Fri Jun 17 17:30:36 2011 +0200
1.2 +++ b/include/vigra/mathutil.hxx Wed Jun 22 15:57:47 2011 +0200
1.3 @@ -1207,6 +1207,199 @@
1.4 return ga;
1.5 }
1.6
1.7 +/*
1.8 + * the following code is derived from lgamma_r() by Sun
1.9 + *
1.10 + * ====================================================
1.11 + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1.12 + *
1.13 + * Developed at SunPro, a Sun Microsystems, Inc. business.
1.14 + * Permission to use, copy, modify, and distribute this
1.15 + * software is freely granted, provided that this notice
1.16 + * is preserved.
1.17 + * ====================================================
1.18 + *
1.19 + */
1.20 +template <class REAL>
1.21 +REAL loggammaImpl(REAL x)
1.22 +{
1.23 + vigra_precondition(x > 0.0,
1.24 + "loggamma(): argument must be positive.");
1.25 +
1.26 + vigra_precondition(x <= 1.0e307,
1.27 + "loggamma(): argument must not exceed 1e307.");
1.28 +
1.29 + double res;
1.30 +
1.31 + if (x < 4.2351647362715017e-22)
1.32 + {
1.33 + res = -std::log(x);
1.34 + }
1.35 + else if ((x == 2.0) || (x == 1.0))
1.36 + {
1.37 + res = 0.0;
1.38 + }
1.39 + else if (x < 2.0)
1.40 + {
1.41 + static const double a[] = { 7.72156649015328655494e-02,
1.42 + 3.22467033424113591611e-01,
1.43 + 6.73523010531292681824e-02,
1.44 + 2.05808084325167332806e-02,
1.45 + 7.38555086081402883957e-03,
1.46 + 2.89051383673415629091e-03,
1.47 + 1.19270763183362067845e-03,
1.48 + 5.10069792153511336608e-04,
1.49 + 2.20862790713908385557e-04,
1.50 + 1.08011567247583939954e-04,
1.51 + 2.52144565451257326939e-05,
1.52 + 4.48640949618915160150e-05 };
1.53 + static const double t[] = { 4.83836122723810047042e-01,
1.54 + -1.47587722994593911752e-01,
1.55 + 6.46249402391333854778e-02,
1.56 + -3.27885410759859649565e-02,
1.57 + 1.79706750811820387126e-02,
1.58 + -1.03142241298341437450e-02,
1.59 + 6.10053870246291332635e-03,
1.60 + -3.68452016781138256760e-03,
1.61 + 2.25964780900612472250e-03,
1.62 + -1.40346469989232843813e-03,
1.63 + 8.81081882437654011382e-04,
1.64 + -5.38595305356740546715e-04,
1.65 + 3.15632070903625950361e-04,
1.66 + -3.12754168375120860518e-04,
1.67 + 3.35529192635519073543e-04 };
1.68 + static const double u[] = { -7.72156649015328655494e-02,
1.69 + 6.32827064025093366517e-01,
1.70 + 1.45492250137234768737e+00,
1.71 + 9.77717527963372745603e-01,
1.72 + 2.28963728064692451092e-01,
1.73 + 1.33810918536787660377e-02 };
1.74 + static const double v[] = { 0.0,
1.75 + 2.45597793713041134822e+00,
1.76 + 2.12848976379893395361e+00,
1.77 + 7.69285150456672783825e-01,
1.78 + 1.04222645593369134254e-01,
1.79 + 3.21709242282423911810e-03 };
1.80 + static const double tc = 1.46163214496836224576e+00;
1.81 + static const double tf = -1.21486290535849611461e-01;
1.82 + static const double tt = -3.63867699703950536541e-18;
1.83 + if (x <= 0.9)
1.84 + {
1.85 + res = -std::log(x);
1.86 + if (x >= 0.7316)
1.87 + {
1.88 + double y = 1.0-x;
1.89 + double z = y*y;
1.90 + double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1.91 + double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1.92 + double p = y*p1+p2;
1.93 + res += (p-0.5*y);
1.94 + }
1.95 + else if (x >= 0.23164)
1.96 + {
1.97 + double y = x-(tc-1.0);
1.98 + double z = y*y;
1.99 + double w = z*y;
1.100 + double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1.101 + double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1.102 + double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1.103 + double p = z*p1-(tt-w*(p2+y*p3));
1.104 + res += (tf + p);
1.105 + }
1.106 + else
1.107 + {
1.108 + double y = x;
1.109 + double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1.110 + double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1.111 + res += (-0.5*y + p1/p2);
1.112 + }
1.113 + }
1.114 + else
1.115 + {
1.116 + res = 0.0;
1.117 + if (x >= 1.7316)
1.118 + {
1.119 + double y = 2.0-x;
1.120 + double z = y*y;
1.121 + double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1.122 + double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1.123 + double p = y*p1+p2;
1.124 + res += (p-0.5*y);
1.125 + }
1.126 + else if(x >= 1.23164)
1.127 + {
1.128 + double y = x-tc;
1.129 + double z = y*y;
1.130 + double w = z*y;
1.131 + double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1.132 + double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1.133 + double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1.134 + double p = z*p1-(tt-w*(p2+y*p3));
1.135 + res += (tf + p);
1.136 + }
1.137 + else
1.138 + {
1.139 + double y = x-1.0;
1.140 + double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1.141 + double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1.142 + res += (-0.5*y + p1/p2);
1.143 + }
1.144 + }
1.145 + }
1.146 + else if(x < 8.0)
1.147 + {
1.148 + static const double s[] = { -7.72156649015328655494e-02,
1.149 + 2.14982415960608852501e-01,
1.150 + 3.25778796408930981787e-01,
1.151 + 1.46350472652464452805e-01,
1.152 + 2.66422703033638609560e-02,
1.153 + 1.84028451407337715652e-03,
1.154 + 3.19475326584100867617e-05 };
1.155 + static const double r[] = { 0.0,
1.156 + 1.39200533467621045958e+00,
1.157 + 7.21935547567138069525e-01,
1.158 + 1.71933865632803078993e-01,
1.159 + 1.86459191715652901344e-02,
1.160 + 7.77942496381893596434e-04,
1.161 + 7.32668430744625636189e-06 };
1.162 + double i = std::floor(x);
1.163 + double t = 0.0;
1.164 + double y = x-i;
1.165 + double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6]))))));
1.166 + double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6])))));
1.167 + res = 0.5*y+p/q;
1.168 + double z = 1.0;
1.169 + while (i > 2.0)
1.170 + {
1.171 + --i;
1.172 + z *= (y+i);
1.173 + }
1.174 + res += std::log(z);
1.175 + }
1.176 + else if (x < 2.8823037615171174e+17)
1.177 + {
1.178 + static const double w[] = { 4.18938533204672725052e-01,
1.179 + 8.33333333333329678849e-02,
1.180 + -2.77777777728775536470e-03,
1.181 + 7.93650558643019558500e-04,
1.182 + -5.95187557450339963135e-04,
1.183 + 8.36339918996282139126e-04,
1.184 + -1.63092934096575273989e-03 };
1.185 + double t = std::log(x);
1.186 + double z = 1.0/x;
1.187 + double y = z*z;
1.188 + double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6])))));
1.189 + res = (x-0.5)*(t-1.0)+yy;
1.190 + }
1.191 + else
1.192 + {
1.193 + res = x*(std::log(x) - 1.0);
1.194 + }
1.195 +
1.196 + return res;
1.197 +}
1.198 +
1.199 +
1.200 } // namespace detail
1.201
1.202 /*! The gamma function.
1.203 @@ -1225,6 +1418,22 @@
1.204 return detail::gammaImpl(x);
1.205 }
1.206
1.207 + /*! The natural logarithm of the gamma function.
1.208 +
1.209 + This function is based on a free implementation by Sun Microsystems, Inc., see
1.210 + <a href="http://www.sourceware.org/cgi-bin/cvsweb.cgi/~checkout~/src/newlib/libm/mathfp/er_lgamma.c?rev=1.6&content-type=text/plain&cvsroot=src">sourceware.org</a> archive. It can be removed once all compilers support the new C99
1.211 + math functions.
1.212 +
1.213 + The argument must be positive and < 1e30. An exception is thrown when these conditions are violated.
1.214 +
1.215 + <b>\#include</b> \<vigra/mathutil.hxx\><br>
1.216 + Namespace: vigra
1.217 + */
1.218 +inline double loggamma(double x)
1.219 +{
1.220 + return detail::loggammaImpl(x);
1.221 +}
1.222 +
1.223
1.224 namespace detail {
1.225
2.1 --- a/test/math/test.cxx Fri Jun 17 17:30:36 2011 +0200
2.2 +++ b/test/math/test.cxx Wed Jun 22 15:57:47 2011 +0200
2.3 @@ -431,7 +431,28 @@
2.4 try { vigra::gamma(0.0); failTest("No exception thrown"); } catch(vigra::PreconditionViolation &) {}
2.5 try { vigra::gamma(-1.0); failTest("No exception thrown"); } catch(vigra::PreconditionViolation &) {}
2.6
2.7 - double args[5] = {0.0, 1.0, 0.7, -0.7, -1.0};
2.8 + shouldEqual(vigra::loggamma(1.0), 0.0);
2.9 + shouldEqual(vigra::loggamma(2.0), 0.0);
2.10 + shouldEqualTolerance(vigra::loggamma(4.0e-22), 49.2705776847491144296, 1e-15);
2.11 + shouldEqualTolerance(vigra::loggamma(0.1), 2.2527126517342055401, 1e-15);
2.12 + shouldEqualTolerance(vigra::loggamma(0.3), 1.0957979948180756047, 1e-15);
2.13 + shouldEqualTolerance(vigra::loggamma(0.8), 0.15205967839983755563, 1e-15);
2.14 + shouldEqualTolerance(vigra::loggamma(1.1), -0.049872441259839757344, 1e-15);
2.15 + shouldEqualTolerance(vigra::loggamma(1.3), -0.10817480950786048655, 1e-15);
2.16 + shouldEqualTolerance(vigra::loggamma(1.8), -0.071083872914372153717, 1e-15);
2.17 + shouldEqualTolerance(vigra::loggamma(3.0), 0.69314718055994528623, 1e-15);
2.18 + shouldEqualTolerance(vigra::loggamma(3.1), 0.78737508327386251938, 1e-15);
2.19 + shouldEqualTolerance(vigra::loggamma(4.0), 1.79175946922805500081, 1e-15);
2.20 + shouldEqualTolerance(vigra::loggamma(8.0), 8.5251613610654143002, 1e-15);
2.21 + shouldEqualTolerance(vigra::loggamma(1000.0), 5905.2204232091812118261, 1e-15);
2.22 + shouldEqualTolerance(vigra::loggamma(1000.2), 5906.6018942569799037, 1e-15);
2.23 + shouldEqualTolerance(vigra::loggamma(2.8e+17), 1.096859847946237952e+19, 1e-15);
2.24 + shouldEqualTolerance(vigra::loggamma(2.9e+17), 1.1370510622188449792e+19, 1e-15);
2.25 + shouldEqualTolerance(vigra::loggamma(5.7646075230342349e+17), 2.2998295812288974848e+19, 1e-15);
2.26 + try { vigra::loggamma(0.0); failTest("No exception thrown"); } catch(vigra::PreconditionViolation &) {}
2.27 + try { vigra::loggamma(-1.0); failTest("No exception thrown"); } catch(vigra::PreconditionViolation &) {}
2.28 +
2.29 + double args[5] = {0.0, 1.0, 0.7, -0.7, -1.0};
2.30 for(int i=0; i<5; ++i)
2.31 {
2.32 double x = args[i], x2 = x*x;