added reintroduced (cleaned) erf,erfc,polevl implementations

Src/MMTK_DCD.c

@@ -131,7 +131,7 @@ readDCD(PyObject *dummy, PyObject *args)
     goto error;
   }
   if( atoms != dcdAtoms ){
-    snprintf(buffer, sizeof(buffer),
+    PyOS_snprintf(buffer, sizeof(buffer),
              "number of atoms in DCD file (%d) doesn't "
 	    "match universe (%d)", dcdAtoms, atoms);
     PyErr_SetString(PyExc_ValueError, buffer);

Src/MMTK_surface.c

@@ -29,7 +29,7 @@ add_point(const TPoint p, TPoint *r, int r_index, PyObject *pt_dict)
 {
   PyObject *p_orig;
   char key[512];
-  snprintf(key, sizeof(key), "%.6f,%.6f,%.6f", p[0], p[1], p[2]);
+  PyOS_snprintf(key, sizeof(key), "%.6f,%.6f,%.6f", p[0], p[1], p[2]);
   p_orig = PyDict_GetItemString(pt_dict, key);
   if(p_orig)
     Py_DECREF(p_orig);
@@ -178,7 +178,7 @@ nbor_data_1_atom(PyObject *nbors, int i, PyObject *atom_data,
     {
       char key[128];
       PyObject *alist;
-      snprintf(key, sizeof(key), "%d %d %d", key0+nbor_boxes[boxn][0],
+      PyOS_snprintf(key, sizeof(key), "%d %d %d", key0+nbor_boxes[boxn][0],
                key1+nbor_boxes[boxn][1], key2+nbor_boxes[boxn][2]);
       alist = PyDict_GetItemString(boxes, key);
       if(!alist && i == -1) printf("none in list at %s, atom %d\n", key, i);
@@ -421,7 +421,7 @@ FindNeighbors(PyObject *dummy, PyObject *args)
     PyObject *v, *tmp;
     char key[128];
     PyObject *pos = PyList_GetItem(atom_data, (Py_ssize_t)i);
-    snprintf(key, sizeof(key), "%d %d %d",
+    PyOS_snprintf(key, sizeof(key), "%d %d %d",
              (int)floor(PyFloat_AsDouble(PyTuple_GetItem(pos,
 						     (Py_ssize_t)0))/box_size),
              (int)floor(PyFloat_AsDouble(PyTuple_GetItem(pos,
@@ -484,7 +484,7 @@ FindNeighbors(PyObject *dummy, PyObject *args)
     {
       char key[128];
       PyObject *alist;
-      snprintf(key, sizeof(key), "%d %d %d", key0+nbor_boxes[boxn][0],
+      PyOS_snprintf(key, sizeof(key), "%d %d %d", key0+nbor_boxes[boxn][0],
                key1+nbor_boxes[boxn][1], key2+nbor_boxes[boxn][2]);
       alist = PyDict_GetItemString(boxes, key);
       if(!alist && i == -1) printf("none in list at %s\n", key);
@@ -590,7 +590,7 @@ FindNeighborsOfAtom(PyObject *dummy, PyObject *args)
     {
       char key[128];
       PyObject *alist;
-      snprintf(key, sizeof(key),
+      PyOS_snprintf(key, sizeof(key),
                "%d %d %d", key0+nbor_boxes[boxn][0],
                key1+nbor_boxes[boxn][1], key2+nbor_boxes[boxn][2]);
       alist = PyDict_GetItemString(boxes, key);

Src/MMTK_trajectory.c

@@ -769,7 +769,7 @@ PyTrajectory_TimeStamp(PyTrajectoryObject *self, char *text)
 {
   time_t now = time(NULL);
   static char time_stamp[200];
-  snprintf(time_stamp, sizeof(time_stamp), text, ctime(&now));
+  PyOS_snprintf(time_stamp, sizeof(time_stamp), text, ctime(&now));
   time_stamp[strlen(time_stamp)-1] = '\0';
   return PyNetCDFFile_AddHistoryLine(self->file, time_stamp);
 }
@@ -1397,14 +1397,14 @@ PyTrajectory_Output(PyTrajectoryOutputSpec *spec, int step,
 	  PyEval_RestoreThread(*thread);
 	if (step > 0)
 	  PyFile_WriteString("\n", spec->destination);
-	snprintf(buffer, sizeof(buffer), "Step %d\n", step);
+	PyOS_snprintf(buffer, sizeof(buffer), "Step %d\n", step);
 	PyFile_WriteString(buffer, spec->destination);
 
 	for (var = data; var->name != NULL; var++)
 	  if (spec->what & var->class) {
 	    switch (var->type) {
 	    case PyTrajectory_Scalar:
-	      snprintf(buffer, sizeof(buffer), var->text, *var->value.dp);
+	      PyOS_snprintf(buffer, sizeof(buffer), var->text, *var->value.dp);
 	      PyFile_WriteString(buffer, spec->destination);
 	      break;
 	    case PyTrajectory_ParticleScalar:
@@ -1412,7 +1412,7 @@ PyTrajectory_Output(PyTrajectoryOutputSpec *spec, int step,
 	      PyFile_WriteString(buffer, spec->destination);
 	      array_data = (double *)var->value.array->data;
 	      for (i = 0; i < var->value.array->dimensions[0]; i++) {
-		snprintf(buffer, sizeof(buffer),
+		PyOS_snprintf(buffer, sizeof(buffer),
                          "  %5d: %g\n", i, *array_data++);
 		PyFile_WriteString(buffer, spec->destination);
 	      }
@@ -1422,7 +1422,7 @@ PyTrajectory_Output(PyTrajectoryOutputSpec *spec, int step,
 	      PyFile_WriteString(buffer, spec->destination);
 	      array_data = (double *)var->value.array->data;
 	      for (i = 0; i < var->value.array->dimensions[0]; i++) {
-		snprintf(buffer, sizeof(buffer),"  %5d: %g, %g, %g\n", i,
+		PyOS_snprintf(buffer, sizeof(buffer),"  %5d: %g, %g, %g\n", i,
                          array_data[0], array_data[1], array_data[2]);
 		PyFile_WriteString(buffer, spec->destination);
 		array_data += 3;
@@ -1433,12 +1433,12 @@ PyTrajectory_Output(PyTrajectoryOutputSpec *spec, int step,
 	      PyFile_WriteString(buffer, spec->destination);
 	      array_data = var->value.dp;
 	      for (i = 0; i < var->length; i++) {
-		snprintf(buffer, sizeof(buffer), " %g", *array_data++);
+		PyOS_snprintf(buffer, sizeof(buffer), " %g", *array_data++);
 		PyFile_WriteString(buffer, spec->destination);
 		if (i == var->length-1)
-		  snprintf(buffer, sizeof(buffer), "\n");
+		  PyOS_snprintf(buffer, sizeof(buffer), "\n");
 		else
-		  snprintf(buffer, sizeof(buffer), ",");
+		  PyOS_snprintf(buffer, sizeof(buffer), ",");
 		PyFile_WriteString(buffer, spec->destination);
 	      }
 	      break;

Src/ReadDCD.c

@@ -607,7 +607,7 @@ int write_dcdheader(FILE *fd, char *filename, int N, int NSET, int ISTART,
 	out_integer = 2;
 	fwrite((char *) & out_integer, sizeof(int), 1, fd);
 
-	snprintf(title_string, sizeof(title_string),
+	PyOS_snprintf(title_string, sizeof(title_string),
                  "REMARKS FILENAME=%s CREATED BY VMD", filename);
 	pad(title_string, 80);
 	fwrite(title_string, sizeof(char), 80, fd);
@@ -616,7 +616,7 @@ int write_dcdheader(FILE *fd, char *filename, int N, int NSET, int ISTART,
 	tmbuf=localtime(&cur_time);
 	strftime(time_str, 10, "%m/%d/%y", tmbuf);
 
-	snprintf(title_string, sizeof(title_string),
+	PyOS_snprintf(title_string, sizeof(title_string),
                  "REMARKS DATE: %s CREATED BY MMTK.", time_str);
 	pad(title_string, 80);
 	fwrite(title_string, sizeof(char), 80, fd);

Src/ewald.c

@@ -17,8 +17,7 @@
  * Include erfc code, unless user wants to use an erfc from libm.
  */
 #ifndef LIBM_HAS_ERFC
-#include "polevl.c"
-#include "ndtr.c"
+#include "ndtr.h"
 #endif
 #ifndef M_PI
 #define M_PI   3.14159265358979323846264338327950288

Src/mconf.h

@@ -161,11 +161,11 @@ typedef struct
 /* ANSI function prototypes */
 #define ANSIPROT
 
-extern double polevl ( double x, double *P, int N );
-extern double p1evl ( double x, double *P, int N );
-extern double ndtr ( double a );
-extern double erfc ( double a );
-extern double erf ( double x );
+//extern double polevl ( double x, double *P, int N );
+//extern double p1evl ( double x, double *P, int N );
+//extern double ndtr ( double a );
+//extern double erfc ( double a );
+//extern double erf ( double x );
 
 /* No error reporting */
 #define mtherr(function, type)

Src/ndtr.c

@@ -145,7 +145,7 @@ Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */
 
 
-#include "mconf.h"
+#include "polevl.h"
 
 
 #ifdef UNK
@@ -378,13 +378,14 @@ static unsigned short U[] = {
 #define UTHRESH 37.519379347
 #endif
 
+/*
 #ifndef ANSIPROT
 double polevl(), p1evl(), exp(), log(), fabs();
 double erf(), erfc();
 #endif
+*/
 
-double ndtr(a)
-double a;
+double ndtr(double a)
 {
 double x, y, z;
 
@@ -406,8 +407,7 @@ return(y);
 }
 
 
-double erfc(a)
-double a;
+double erfc(double a)
 {
 double p,q,x,y,z;
 
@@ -457,8 +457,7 @@ return(y);
 
 
 
-double erf(x)
-double x;
+double erf(double x)
 {
 double y, z;
 

Src/ndtr.h

+/*							ndtr.c
+ *
+ *	Normal distribution function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, ndtr();
+ *
+ * y = ndtr( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the area under the Gaussian probability density
+ * function, integrated from minus infinity to x:
+ *
+ *                            x
+ *                             -
+ *                   1        | |          2
+ *    ndtr(x)  = ---------    |    exp( - t /2 ) dt
+ *               sqrt(2pi)  | |
+ *                           -
+ *                          -inf.
+ *
+ *             =  ( 1 + erf(z) ) / 2
+ *             =  erfc(z) / 2
+ *
+ * where z = x/sqrt(2). Computation is via the functions
+ * erf and erfc.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC      -13,0         8000       2.1e-15     4.8e-16
+ *    IEEE     -13,0        30000       3.4e-14     6.7e-15
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition         value returned
+ * erfc underflow    x > 37.519379347       0.0
+ *
+ */
+/*							erf.c
+ *
+ *	Error function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, erf();
+ *
+ * y = erf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * The integral is
+ *
+ *                           x 
+ *                            -
+ *                 2         | |          2
+ *   erf(x)  =  --------     |    exp( - t  ) dt.
+ *              sqrt(pi)   | |
+ *                          -
+ *                           0
+ *
+ * The magnitude of x is limited to 9.231948545 for DEC
+ * arithmetic; 1 or -1 is returned outside this range.
+ *
+ * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
+ * erf(x) = 1 - erfc(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0,1         14000       4.7e-17     1.5e-17
+ *    IEEE      0,1         30000       3.7e-16     1.0e-16
+ *
+ */
+/*							erfc.c
+ *
+ *	Complementary error function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, erfc();
+ *
+ * y = erfc( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ *  1 - erf(x) =
+ *
+ *                           inf. 
+ *                             -
+ *                  2         | |          2
+ *   erfc(x)  =  --------     |    exp( - t  ) dt
+ *               sqrt(pi)   | |
+ *                           -
+ *                            x
+ *
+ *
+ * For small x, erfc(x) = 1 - erf(x); otherwise rational
+ * approximations are computed.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0, 9.2319   12000       5.1e-16     1.2e-16
+ *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition              value returned
+ * erfc underflow    x > 9.231948545 (DEC)       0.0
+ *
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.2:  June, 1992
+Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+extern double ndtr(double a);
+extern double erfc(double a);
+extern double erf(double x);
+

Src/nonbonded.c

@@ -7,6 +7,10 @@
 #define _FORCEFIELD_MODULE
 #define PY_ARRAY_UNIQUE_SYMBOL PyArray_MMTKFF_API
 
+#ifndef LIBM_HAS_ERFC
+#include "ndtr.h"
+#endif
+
 #include "MMTK/forcefield.h"
 #include "MMTK/forcefield_private.h"
 

Src/polevl.c

@@ -50,6 +50,8 @@ Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */
 
 
+#include "polevl.h"
+
 double polevl( double x, double coef[], int N )
 {
 double ans;

Src/polevl.h

+/*							polevl.c
+ *							p1evl.c
+ *
+ *	Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * double x, y, coef[N+1], polevl[];
+ *
+ * y = polevl( x, coef, N );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates polynomial of degree N:
+ *
+ *                     2          N
+ * y  =  C  + C x + C x  +...+ C x
+ *        0    1     2          N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C  , ..., coef[N] = C  .
+ *            N                   0
+ *
+ *  The function p1evl() assumes that coef[N] = 1.0 and is
+ * omitted from the array.  Its calling arguments are
+ * otherwise the same as polevl().
+ *
+ *
+ * SPEED:
+ *
+ * In the interest of speed, there are no checks for out
+ * of bounds arithmetic.  This routine is used by most of
+ * the functions in the library.  Depending on available
+ * equipment features, the user may wish to rewrite the
+ * program in microcode or assembly language.
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.1:  December, 1988
+Copyright 1984, 1987, 1988 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+extern double polevl( double x, double coef[], int N );
+
+extern double p1evl( double x, double coef[], int N );
+

setup.py

 #!/usr/bin/env python
-from click.exceptions import FileError
 
 package_name = "MMTK"