This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
Patch to define GSL_SF_FACT_NMAX and GSL_SF_DOUBLEFACT_NMAX
- From: "Jonathan G. Underwood" <j dot underwood at open dot ac dot uk>
- To: gsl-discuss at sources dot redhat dot com
- Date: Mon, 10 Oct 2005 11:13:54 +0100
- Subject: Patch to define GSL_SF_FACT_NMAX and GSL_SF_DOUBLEFACT_NMAX
Hi
The attached patch exports the following two macros:
GSL_SF_FACT_NMAX
GSL_SF_DOUBLEFACT_NMAX
which define the maximum values of n such that gsl_sf_fact(n) and
gsl_sf_doublefact(n) do not over flow. Macros were internally defined
for this previously, but not exported or exported or named consistently
with other macros. Since gsl has the equivalent GSL_SF_GAMMA_XMAX for
gsl_sf_gamma, I think we should define these macros as they are useful
from a user perspective, and it makes for consistency for the functions
related to the gamma function. Also, this patch adds an #undef
LogRootTwoPi_ for gamma.c, as it was never meant to be exported
(cosmetic). The patch also documents the changes.
Jonathan.
--
------------------------------------------------
Dr Jonathan Underwood
The Department of Physics and Astronomy
The Open University
Walton Hall
Milton Keynes
MK7 6AA
UK
Tel: +44 (0) 1908 652514
Fax: +44 (0) 1908 654192
Index: specfunc/gsl_sf_gamma.h
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/gsl_sf_gamma.h,v
retrieving revision 1.38
diff -u -r1.38 gsl_sf_gamma.h
--- specfunc/gsl_sf_gamma.h 26 Jun 2005 13:27:08 -0000 1.38
+++ specfunc/gsl_sf_gamma.h 9 Oct 2005 23:40:45 -0000
@@ -280,6 +280,11 @@
*/
#define GSL_SF_GAMMA_XMAX 171.0
+/* The maximum n such that gsl_sf_fact(n) does not give an overflow. */
+#define GSL_SF_FACT_NMAX 170
+
+/* The maximum n such that gsl_sf_doublefact(n) does not give an overflow. */
+#define GSL_SF_DOUBLEFACT_NMAX 297
__END_DECLS
Index: specfunc/gamma.c
===================================================================
RCS file: /cvs/gsl/gsl/specfunc/gamma.c,v
retrieving revision 1.72
diff -u -r1.72 gamma.c
--- specfunc/gamma.c 26 Jun 2005 13:27:07 -0000 1.72
+++ specfunc/gamma.c 9 Oct 2005 23:40:45 -0000
@@ -39,9 +39,7 @@
/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
-#define FACT_TABLE_MAX 170
-#define FACT_TABLE_SIZE (FACT_TABLE_MAX+1)
-static struct {int n; double f; long i; } fact_table[FACT_TABLE_SIZE] = {
+static struct {int n; double f; long i; } fact_table[GSL_SF_FACT_NMAX + 1] = {
{ 0, 1.0, 1L },
{ 1, 1.0, 1L },
{ 2, 2.0, 2L },
@@ -250,9 +248,7 @@
*/
};
-#define DOUB_FACT_TABLE_MAX 297
-#define DOUB_FACT_TABLE_SIZE (DOUB_FACT_TABLE_MAX+1)
-static struct {int n; double f; long i; } doub_fact_table[DOUB_FACT_TABLE_SIZE] = {
+static struct {int n; double f; long i; } doub_fact_table[GSL_SF_DOUBLEFACT_NMAX + 1] = {
{ 0, 1.000000000000000000000000000, 1L },
{ 1, 1.000000000000000000000000000, 1L },
{ 2, 2.000000000000000000000000000, 2L },
@@ -725,7 +721,6 @@
return GSL_SUCCESS;
}
-
/* x = eps near zero
* gives double-precision for |eps| < 0.02
*/
@@ -1029,7 +1024,7 @@
result->val = 1.77245385090551602729817;
result->err = GSL_DBL_EPSILON * result->val;
return GSL_SUCCESS;
- } else if (x <= (FACT_TABLE_MAX + 1.0) && x == floor(x)) {
+ } else if (x <= (GSL_SF_FACT_NMAX + 1.0) && x == floor(x)) {
int n = (int) floor (x);
result->val = fact_table[n - 1].f;
result->err = GSL_DBL_EPSILON * result->val;
@@ -1478,7 +1473,7 @@
result->err = 0.0;
return GSL_SUCCESS;
}
- else if(n <= FACT_TABLE_MAX){
+ else if(n <= GSL_SF_FACT_NMAX){
result->val = fact_table[n].f;
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
@@ -1498,7 +1493,7 @@
result->err = 0.0;
return GSL_SUCCESS;
}
- else if(n <= DOUB_FACT_TABLE_MAX){
+ else if(n <= GSL_SF_DOUBLEFACT_NMAX){
result->val = doub_fact_table[n].f;
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
@@ -1513,7 +1508,7 @@
{
/* CHECK_POINTER(result) */
- if(n <= FACT_TABLE_MAX){
+ if(n <= GSL_SF_FACT_NMAX){
result->val = log(fact_table[n].f);
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
@@ -1529,7 +1524,7 @@
{
/* CHECK_POINTER(result) */
- if(n <= DOUB_FACT_TABLE_MAX){
+ if(n <= GSL_SF_DOUBLEFACT_NMAX){
result->val = log(doub_fact_table[n].f);
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
@@ -1589,7 +1584,7 @@
result->err = 0.0;
return GSL_SUCCESS;
}
- else if (n <= FACT_TABLE_MAX) {
+ else if (n <= GSL_SF_FACT_NMAX) {
result->val = (fact_table[n].f / fact_table[m].f) / fact_table[n-m].f;
result->err = 6.0 * GSL_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
@@ -1681,3 +1676,6 @@
{
EVAL_RESULT(gsl_sf_lnchoose_e(n, m, &result));
}
+
+#undef LogRootTwoPi_
+
Index: doc/specfunc-gamma.texi
===================================================================
RCS file: /cvs/gsl/gsl/doc/specfunc-gamma.texi,v
retrieving revision 1.27
diff -u -r1.27 specfunc-gamma.texi
--- doc/specfunc-gamma.texi 13 Sep 2005 09:43:24 -0000 1.27
+++ doc/specfunc-gamma.texi 9 Oct 2005 23:40:46 -0000
@@ -117,6 +117,9 @@
@cindex factorial
These routines compute the factorial @math{n!}. The factorial is
related to the Gamma function by @math{n! = \Gamma(n+1)}.
+The maximum value of @math{n} such that @math{n!} is not
+considered an overflow is given by the macro @code{GSL_SF_FACT_NMAX}
+and is 170.
@comment exceptions: GSL_EDOM, GSL_OVRFLW
@end deftypefun
@@ -124,6 +127,9 @@
@deftypefunx int gsl_sf_doublefact_e (unsigned int @var{n}, gsl_sf_result * @var{result})
@cindex double factorial
These routines compute the double factorial @math{n!! = n(n-2)(n-4) \dots}.
+The maximum value of @math{n} such that @math{n!!} is not
+considered an overflow is given by the macro @code{GSL_SF_DOUBLEFACT_NMAX}
+and is 297.
@comment exceptions: GSL_EDOM, GSL_OVRFLW
@end deftypefun