? polyfit.patch
Index: polyfit.m
===================================================================
RCS file: /cvs/octave/scripts/polynomial/polyfit.m,v
retrieving revision 1.23
diff -u -r1.23 polyfit.m
--- polyfit.m	23 Aug 2005 18:38:28 -0000	1.23
+++ polyfit.m	9 Aug 2006 19:24:50 -0000
@@ -78,10 +78,16 @@
   x = reshape (x, l, 1);
   y = reshape (y, l, 1);
 
+  ## scale the x values
+  xmax = max(abs(x));
+  x = x/xmax;
+
+
   X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0));
 
   p = X \ y;
 
+
   if (nargout > 1)
 
     yf = X*p;
@@ -99,8 +105,29 @@
 
   endif
 
-  ## Return value should be a row vector.
-
-  p = p.';
+  ## Scale the polynomial back and return a row vector
+  power = length(p)-1:-1:0;
+  power = (xmax .^ power)(:);
+  p = (p ./ power).';
 
 endfunction
+
+## test simple polynomial; here we check the returned polynomial
+%!test				
+%!  x = [1,2,3];  
+%!  p = poly(x);  
+%!  y = polyval(p,x);
+%!  pnew = polyfit([0,1,2,3]', [-6,0,0,0]', 3);
+%!  assert(p, pnew, 10^-6);
+
+## test difficult case where scaling is really needed
+## here we check the values after the evaluation of the new polynomial
+%!test
+%!  x = [ -1196.4, -1195.2, -1194, -1192.8, -1191.6, -1190.4, -1189.2, -1188, \
+%!        -1186.8, -1185.6, -1184.4, -1183.2, -1182 ];
+%!  y = [ 315571.7086, 315575.9618, 315579.4195, 315582.6206, 315585.4966,    \
+%!        315588.3172, 315590.9326, 315593.5934, 315596.0455, 315598.4201,    \
+%!        315600.7143, 315602.9508, 315605.1765 ];
+%!  p = polyfit(x,y,10);
+%!  ynew = polyval(p,x);
+%!  assert(norm(y-ynew), 0.08,  0.01);