/* 
   Fuzzy substraction of two curves. Substracts the Y values of _op_ to
   the current Function, by making sure that the Y value substracted to
   a given point corresponds to the closest X_ value of the point in _op_.
   This function somehow assumes that the data is reasonably organised,
   and will never go backwards to find a matching X value in _op_.

   In any case, you really should consider using split_monotonic on it first.
 */

static VALUE function_fuzzy_substract(VALUE self, VALUE op)
{
  long ss = function_sanity_check(self);
  const double *xs = Dvector_Data_for_Read(get_x_vector(self),NULL);
  double *ys = Dvector_Data_for_Write(get_y_vector(self),NULL);
  long so = function_sanity_check(op);
  const double *xo = Dvector_Data_for_Read(get_x_vector(op),NULL);
  const double *yo = Dvector_Data_for_Read(get_y_vector(op),NULL);
  long i,j = 0;
  double diff;
  double fuzz = 0;              /* The actual sum of the terms */
  
  for(i = 0; i < ss; i++) 
    {
      /* We first look for the closest point */
      diff = fabs(xs[i] - xo[j]);
      while((j < (so - 1)) && (fabs(xs[i] - xo[j+1]) <  diff))
        diff = fabs(xs[i] - xo[++j]);
      fuzz += diff;
      ys[i] -= yo[j];
    }
  return rb_float_new(fuzz);
}