#include "math.h"

#include <Python.h>

#include <numpy/arrayobject.h>

/*
Get the include directories within python using

  import distutils.sysconfig
  print distutils.sysconfig.get_python_inc()
  import numpy as np 
  print np.get_include()

Linux: 

gcc  -fPIC -shared -g -Wall -O3 \
   -I /usr/include/python2.7 -I /usr/lib64/python2.7/site-packages/numpy/core/include \
  mymath.c -o mymath.so 


Mac: 

gcc -Wl,-F. -bundle -undefined dynamic_lookup -g -Wall -O3 \
   -I /usr/include/python2.7 -I /System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/core/include \
   mymath.c -o mymath.so

Windows (windows also needs a .lib and the dll must be renamed) :

cl /LD /O2 mymath.c ^
   /I c:\Users\pret\Anaconda\include ^
   /I C:\Users\pret\Anaconda\lib\site-packages\\numpy\\core\\include ^
   c:\Users\pret\Anaconda\libs\python27.lib 

copy mymath.dll mymath.pyd



*/

static PyObject *
array_sqrt(PyObject *self, PyObject *args)
{
  PyObject * array; 

  /* parse inputs */

  if (!PyArg_ParseTuple(args, "O!",
                        &PyArray_Type, &array))
    return NULL;

  if(!(PyArray_NDIM(array) == 1 &&
       PyArray_TYPE(array) == NPY_FLOAT32 &&
       (PyArray_FLAGS(array) & NPY_C_CONTIGUOUS))) {
    PyErr_SetString(PyExc_TypeError, "input should be c-contiguous 1D float array");
    return NULL;
  }

  const float *input =  PyArray_DATA(array);
  const npy_intp *shape = PyArray_DIMS(array);
  npy_intp n = shape[0];

  /* prepare outputs */
  PyObject *array_out = PyArray_SimpleNew(1, &n, NPY_FLOAT32);
  float *output =  PyArray_DATA(array_out);
  
  /* actual computation */
  npy_intp i; 
  for(i = 0; i < n; i++) output[i] = sqrtf(input[i]);

  return array_out;
}



static PyMethodDef method_list[] = {
  {"sqrt",  array_sqrt, METH_VARARGS, "Computes square root of an array."},
  {NULL, NULL, 0, NULL}        /* Sentinel */
};




PyMODINIT_FUNC initmymath() {
  Py_InitModule3("mymath", method_list, "this module contains specific math functions");
  /* required, otherwise numpy functions do not work */
  import_array();
}