3

I was going through the R source code for the function gamma() which points to this C source file.

The snippet of the code is as follows:

#ifdef NOMORE_FOR_THREADS
    static int ngam = 0;
    static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.;

    /* Initialize machine dependent constants, the first time gamma() is called.
    FIXME for threads ! */
    if (ngam == 0) {
    ngam = chebyshev_init(gamcs, 42, DBL_EPSILON/20);/*was .1*d1mach(3)*/
    gammalims(&xmin, &xmax);/*-> ./gammalims.c */
    xsml = exp(fmax2(log(DBL_MIN), -log(DBL_MAX)) + 0.01);
    /*   = exp(.01)*DBL_MIN = 2.247e-308 for IEEE */
    dxrel = sqrt(DBL_EPSILON);/*was sqrt(d1mach(4)) */
    }
#else
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
 * (xmin, xmax) are non-trivial, see ./gammalims.c
 * xsml = exp(.01)*DBL_MIN
 * dxrel = sqrt(DBL_EPSILON) = 2 ^ -26
*/
# define ngam 22
# define xmin -170.5674972726612
# define xmax  171.61447887182298
# define xsml 2.2474362225598545e-308
# define dxrel 1.490116119384765696e-8
#endif

Now I understand that this piece of code is initializing values of variables like ngam, xsml, dxrel etc. But I don't understand what is NOMORE_FOR_THREADS.

I could not find where that is defined so I am confused as to when would it enter the #ifdef clause or else just go in the #else clause and define the variables according to IEEE standard.

If anyone can shed some light on what is NOMORE_FOR_THREADS or where it is defined, it would be really helpful.

Thanks for the help.

Sid
  • 65
  • 8

1 Answers1

2

Looks like that #ifdef was added in revision 13433:

> svn log -vr13433
------------------------------------------------------------------------
r13433 | maechler | 2001-04-10 10:05:54 -0500 (Tue, 10 Apr 2001) | 2 lines
Changed paths:
   M /trunk/src/nmath/beta.c
   M /trunk/src/nmath/gamma.c
   M /trunk/src/nmath/lgamma.c
   M /trunk/src/nmath/lgammacor.c

globals now #ifdef NOMORE_FOR_THREADS

------------------------------------------------------------------------

And here's the diff with the previous revision:

> svn diff -r13432:13433 src/nmath/gamma.c
Index: src/nmath/gamma.c
===================================================================
--- src/nmath/gamma.c   (revision 13432)
+++ src/nmath/gamma.c   (revision 13433)
@@ -1,7 +1,7 @@
 /*
  *  Mathlib : A C Library of Special Functions
  *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000 The R Development Core Team
+ *  Copyright (C) 2000-2001 The R Development Core Team
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -40,7 +40,7 @@

 double gammafn(double x)
 {
-    static /* const */ double gamcs[42] = {
+    const double gamcs[42] = {
    +.8571195590989331421920062399942e-2,
    +.4415381324841006757191315771652e-2,
    +.5685043681599363378632664588789e-1,
@@ -89,6 +89,7 @@
     double y;
     double sinpiy, value;

+#ifdef NOMORE_FOR_THREADS
     static int ngam = 0;
     static double xmin = 0, xmax = 0., xsml = 0., dxrel = 0.;

@@ -101,6 +102,18 @@
    /*   = exp(.01)*DBL_MIN = 2.247e-308 for IEEE */
    dxrel = sqrt(1/DBL_EPSILON);/*was (1/d1mach(4)) */
     }
+#else
+/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
+ * (xmin, xmax) are non-trivial, see ./gammalims.c 
+ * xsml = exp(.01)*DBL_MIN 
+ * dxrel = sqrt(1/DBL_EPSILON) = 2 ^ 26
+*/
+# define ngam 22
+# define xmin -170.5674972726612
+# define xmax  171.61447887182298
+# define xsml 2.2474362225598545e-308
+# define dxrel 67108864.
+#endif

     if(ISNAN(x)) return x;

I didn't see NOMORE_FOR_THREADS defined anywhere in the source code in that revision, the next couple hundred revisions after it was added, or in current revisions. It's possible (though seems unlikely) that it could be defined in some system header file. It seems more likely that it would be specified via a command line argument.

Joshua Ulrich
  • 173,410
  • 32
  • 338
  • 418
  • 1
    is it possible that that piece of code does nothing as R is single threaded? – Sid Nov 10 '15 at 18:36
  • @Sid: possible, and may be to support multi-threaded uses of the R math library from other applications. But that's just speculation. – Joshua Ulrich Nov 10 '15 at 19:37
  • Maybe it's a holdout from the translation from Fortan or something that was just never finished ... A C program called from R could certainly use OpenMP or pthreads so perhaps they never finished implementation. At the same time it looks more like they are checking for precision on a specific platform that may define NOMORE_FOR_THREADS but I realize that's not a definitive answer. – JimLohse Jan 15 '18 at 18:40
  • 1
    Actually found this [here](http://en.cppreference.com/w/cpp/numeric/math/lgamma): "The POSIX version of lgamma is not thread-safe: each execution of the function stores the sign of the gamma function of arg in the static external variable signgam. Some implementations provide lgamma_r, which takes a pointer to user-provided storage for singgam as the second parameter, and is thread-safe. " FWIW – JimLohse Jan 15 '18 at 18:52