/*
  Splits the function into strictly monotonic sub-functions.
  Returns the array of the subfunctions. The returned values are
  necessarily new values.
*/

static VALUE function_split_monotonic(VALUE self)
{
  VALUE ret = rb_ary_new();
  VALUE cur_x = Dvector_Create();
  VALUE cur_y = Dvector_Create();

  long size = function_sanity_check(self);
  long i;
  if(size < 2)
    rb_raise(rb_eRuntimeError, "Function needs to have at least 2 points");

  double *x = Dvector_Data_for_Read(get_x_vector(self),NULL);
  double *y = Dvector_Data_for_Read(get_y_vector(self),NULL);

  double last_x;
  double direction; /* -1 if down, +1 if up, so that the product of 
                       (x - last_x) with direction should always be positive
                    */
  VALUE f;
                     
                       
  /* bootstrap */
  if(x[1] > x[0])
    direction = 1;
  else
    direction = -1;
  last_x = x[1];
  for(i = 0; i < 2; i++)
    {
      Dvector_Push_Double(cur_x, x[i]);
      Dvector_Push_Double(cur_y, y[i]);
    }

  for(i = 2; i < size; i++) 
    {
      if(direction * (x[i] - last_x) <= 0) 
        {
          /* we need to add a new set of Dvectors */
          f = Function_Create(cur_x, cur_y);
          rb_ary_push(ret, f);
          cur_x = Dvector_Create();
          cur_y = Dvector_Create();
          /* We don't store the previous point if 
           the X value is the same*/
          if(x[i] != last_x) 
            {
              Dvector_Push_Double(cur_x, x[i-1]);
              Dvector_Push_Double(cur_y, y[i-1]);
            }
          direction *= -1;
        }
      /* store the current point */
      Dvector_Push_Double(cur_x, x[i]);
      Dvector_Push_Double(cur_y, y[i]);
      last_x = x[i];
    }
  f = Function_Create(cur_x, cur_y);
  rb_ary_push(ret, f);
  return ret;
}