#include <Python.h>
#include <numpy/arrayobject.h>



static unsigned int myrand16(unsigned int next) {
  next = next * 1103515245 + 12345;
  return((unsigned int)(next/65536) % 32768);
}


static unsigned int myrand32(unsigned int prev) {
  unsigned int low = myrand16(prev); 
  unsigned int high = myrand16(low); 

  return (high << 16) | low; 
}


static PyObject *
fill_knapsack(PyObject *self, PyObject *args)
{
  PyObject *array_weights, *array_values; 
  double maxweight; 
  int seed;
  
  /* parse inputs */

  if (!PyArg_ParseTuple(args, "O!O!di",
                        &PyArray_Type, &array_weights, 
                        &PyArray_Type, &array_values, 
                        &maxweight, &seed))
    return NULL;

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

  int n = PyArray_SIZE(array_weights); 
  
  if(!(PyArray_NDIM(array_values) == 1 &&
       PyArray_TYPE(array_values) == NPY_FLOAT64 &&
       PyArray_SIZE(array_values) == n && 
       (PyArray_FLAGS(array_values) & NPY_C_CONTIGUOUS))) {
    PyErr_SetString(PyExc_TypeError, "values should be c-contiguous 1D float array, same size as weights");
    return NULL;
  }
  
  const double *weights =  PyArray_DATA(array_weights);
  const double *values =  PyArray_DATA(array_values);
 
  npy_intp outsize;   
  int *perm = malloc(n * sizeof(int)); 
  double sum_weight = 0, sum_values = 0;
 
  {
    int i;
    for(i = 0; i < n; i++) perm[i] = i;
    
    /* permute until weight is too big */
    unsigned int rand_state = seed; 
    
    for(i = 0; i < n; i++) {
      rand_state = myrand32(rand_state);
      
      {  /* generate next permutation element */
        int j = i + rand_state % (n - i); 
        int tmp = perm[j]; 
        perm[j] = perm[i];
        perm[i] = tmp; 
      }
      
      /* accumulate weights */
      int k = perm[i];
      sum_weight += weights[k];
      
      if(sum_weight > maxweight) break;    
      sum_values += values[k];
    }
    outsize = i;
  }

  /* prepare outputs */
  PyObject *array_out = PyArray_SimpleNew(1, &outsize, NPY_INT32);
  int *output =  PyArray_DATA(array_out);

  int i; 
  for(i = 0; i < outsize; i++) output[i] = perm[i]; 
  
  PyObject *result = Py_BuildValue("(dN)", sum_values, array_out); 

  return result;
}



static PyMethodDef method_list[] = {
  {"fill_knapsack",  fill_knapsack, METH_VARARGS, ""},
  {NULL, NULL, 0, NULL}        /* Sentinel */
};




PyMODINIT_FUNC initc_knapsack() {
  Py_InitModule("c_knapsack", method_list);
  /* required, otherwise numpy functions do not work */
  import_array();
}