/* Computes the derivative of the Function and returns it as a new Function. The newly created function shares the X vector with the previous one. WARNING: this is a very naive 3-points algorithm; you should consider using diff_5p */ static VALUE function_derivative(VALUE self) { long size = function_sanity_check(self); const double *x = Dvector_Data_for_Read(get_x_vector(self),NULL); const double *y = Dvector_Data_for_Read(get_y_vector(self),NULL); VALUE derivative = Dvector_Create(); long i = 0; /* First value */ Dvector_Push_Double(derivative, (y[i+1] - y[i]) /(x[i+1] - x[i])); i++; while(i < (size - 1)) { Dvector_Push_Double(derivative, .5 * ( (y[i+1] - y[i]) /(x[i+1] - x[i]) + (y[i] - y[i-1]) /(x[i] - x[i-1]) )); i++; } Dvector_Push_Double(derivative, (y[i] - y[i-1]) /(x[i] - x[i-1])); return Function_Create(get_x_vector(self), derivative); }