#include <cmath>
#include <octave/oct.h>
#include <octave/lo-ieee.h>
#include <octave/pager.h>
//using namespace std;

// Author : Etienne Grossmann <address@hidden>, 2005

DEFUN_DLD (loop_add, args, ,\
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} address@hidden =} loop_add (@var{a}, @var{indices}, @var{what})\n\
\n\
This function is equivalent to, but is quicker than :\n\
\n\
   @address@hidden; for @var{j}=1:rows(@var{indices}(:)), @var{b}(@var{indices}(@var{j})) += @var{what}(@var{j}); end\n\
\n\
@var{what} is optional and is 1 by default.\n\
\n\
If  @var{what} has less elements than @var{indices}, values of @var{what} are wrapped. I.e. @var{what}(1 + rem(@var{j}-1,rows(@var{what}(:))) is used instead of @var{what}(@var{j}).\n\
@end deftypefn")
{
  // octave_value_list retval;
  if (args.length()<2){
    error ("loop_add: called w/ %i argument(s); 2-3 are needed.",args.length());
    return octave_value ();
  }
  Matrix res(args(0).matrix_value());
  Matrix inds(args(1).matrix_value());
  //Matrix what(args(2).matrix_value());
  Matrix what(1,1,1.0);
  if (args.length()>2) {
    what = args(2).matrix_value();
  }

  int norig = res.columns() * res.rows() - 1;
  int ninds = inds.columns() * inds.rows();
  int nwhat = what.columns() * what.rows();
  //printf ("Size of A:%i, of INDICES:%i, of WHAT:%i\n",norig+1,ninds,nwhat);

   for (int i=0; i<ninds; i++) {

     int j = ((int)inds(i)) - 1;

     if ((j < 0) || (j > norig)){ // Check index
       if (j<0)
	 error ("loop_add: Can't add to a(%i) : index is not positive.",j+1);
       else
	 error ("loop_add: Can't add to a(%i) because a has only %i elements.",j+1,norig+1);
       return octave_value ();
     }

     res(j) =  res(j) + what(i % nwhat);
  }
  return octave_value (res);
}