[Tutor] Extensions in C - numeric arrays & pointers

[Ray, Scott (Dow AgroSciences)]

As part of my Python education, I decided to turn a handful of simple numerical C functions into a small Python extension.  I first selected a simple C function that takes one scalar input (a float or a double) and returns one scalar output (same).  I followed the guidance in Guido and Drake's  "Extending and Embedding" guide and, much to my delight, the extension worked on the first try.

Then I moved to C functions that work with numeric arrays and their associated pointers.  This has not gone so smoothly.  After much experimentation, re-reading of documentation, searching mailing list archives, and scanning through other people's source code, I'm still stuck and in need of some guidance.  So, I now turn to the experts.

Here's a simple C function that is representative of what I'm trying to work with.

	double average(double *xvec, int numpts) {
	  int n;
	  double xbar;

	  xbar = 0.0;
	  for (n = 0; n < numpts; ++n) {
	    xbar += xvec[n];
	  }
	  return xbar/numpts;
	}

Given C's argument passing mechanism (pass by value) , this routine expects its first argument to be a pointer to an array of doubles.

What I want to be able to do is include this C function in a Python extension module (call it "ezmath") and use it in Python as:

	import ezmath
	z = [3.1416, 2.71828, 42]
	ezmath.average(z,3)

I'm thinking that the C wrapper function should look something like,

	static PyObject *
	ezmath_average(PyObject *self, PyObject *args) {
		PyObject *x;
		int num;
		double *xa;
		double y;

		if (!PyArg_ParseTuple(args, "Oi", &x, &num)) {
			return NULL;
		}

		/* big mystery here */

		y = average(xa, num);
		return Py_BuildValue("d",y);
	}

But I haven't figured out what goes in the middle to connect or otherwise convert the Python list + PyObject pointer to the C array and its pointer.  Or, perhaps, my mental model is all wrong and I need to go about this in some other way.

What I have tried is various combinations of:
- using the "array" extension module (under the assumption that Python array objects are somehow "closer" to C arrays)
    e.g. 
	import ezmath
	import array
	z = [3.1416, 2.71828, 42]
	cz = array.array('d',z)
	ezmath.average(cz,3)

- C pointer casts (e.g. xa = (double*) x )
- using PySequence_Size(x) to test whether the argument parsing operation worked (seems to have)
- using PySequence_GetItem(x,i) to attempt to get at the i'th element in the list or array (though the outcome is not making sense to me)

and lots of other more hare-brained ideas.  Nothing seems to work and it feels like I must be missing something simple.

Any help / guidance / pointers / redirection that more experienced Python/C developers can provide will be greatly appreciated. Please note that
- I'm new to Python and to OOP.  My advanced C skills are rusty at best.  Advising me to read the NumPy or SciPy source code probably won't help.  I'd like to work up to this, but I'm not there now.
- My real goal here is learning, I know I could compute a numeric average in a multitude of ways.
- I'm aware of the existence of NumPy and SWIG and of their advantages (for numerical work and extension development, respectively).   If my current Python experiments work out positively, I'll be headed in those directions.  Perhaps the right answer is to go there sooner rather than later.  Right now though I'd prefer, if possible, to solve this puzzle without additional modules and objects types.  It seems to me like there ought to be a solution with Python's lists and/or arrays.  Then again, perhaps I'm being stubborn.
- Finally, I notice in scanning through the mailing list archives that questions on this issue seem to come up periodically.  Maybe I haven't dug far enough, but I have not found the answers provided to be fully illuminating.  If I find instructive ways to solve this problem, I could maybe write a tutorial and post it on one of the Python websites, perhaps tacking on some additional stuff on like converting C result arrays back out to Python, doing it all over with NumPy, etc.

Thanks,

Scott Ray
slray@dow.com