[Bug-gsl] gsl_fit_mul results are WRONG

[blavy pierre]

The doc of gsl_fit_mul states that this functions computes the model  : 
Y = c_1 *X
Therefore, if run with x=1,2,3 and y=100,101,102, it should return 
something around c_1 = 50.
But it returns exactly 1. In fact 1 is the result of the gsl_fit_linear 
function, computed with a not null intercept.

The bug is explained because your implementation assume that c_1's in 
the two following models are equals, which is of course FALSE.
model 1 : y = c_0 +  c_1*x (whant gsl_fit_linear computes)
model 2 : y= 0       + c_1*x (whant gsl_fit_mul SHOULD compute, in fact 
it computes model 1 and drop c_0).

Can you please fix it, or patch the doc.



--- gsl_fit_mul  documentation ---
Function: int gsl_fit_mul (const double * x, const size_t xstride, const 
double * y, const size_t ystride, size_t n, double * c1, double * cov11, 
double * sumsq)
    This function computes the best-fit linear regression coefficient 
c1 of the model Y = c_1 X for the datasets (x, y), two vectors of length 
n with strides xstride and ystride. The errors on y are assumed unknown 
so the variance of the parameter c1 is estimated from the scatter of the 
points around the best-fit line and returned via the parameter cov11. 
The sum of squares of the residuals from the best-fit line is returned 
in sumsq.