diff -r b280628de44e scripts/polynomial/polyfit.m
--- a/scripts/polynomial/polyfit.m	Mon Jan 28 21:18:03 2008 +0100
+++ b/scripts/polynomial/polyfit.m	Wed Jan 30 18:13:39 2008 +0100
@@ -51,17 +51,31 @@
 ## @item yf
 ## The values of the polynomial for each value of @var{x}.
 ## @end table
+## 
+## If a fourth bool parameter @var{scale} is given and set to
+## @code{true}, the values of @var{x} are scaled by their maximum absolute
+## value before doing the actual calculation. This helps in numercial difficult
+## cases. The returned coefficients @var{p} are valid for the unscaled
+## values of @var{x}, while the values in the structure @var{S} are
+## valid for the scaled polynomial used for the actual calculation.
 ## @end deftypefn
 
 ## Author: KH <address@hidden>
 ## Created: 13 December 1994
 ## Adapted-By: jwe
 
-function [p, s, mu] = polyfit (x, y, n)
+function [p, s, mu] = polyfit (x, y, n, scale)
 
+  if (nargin < 3 || nargin > 4)
+    print_usage ();
+  endif
 
-  if (nargin != 3)
-    print_usage ();
+  if (nargin == 3)
+    scale = false;
+  else
+    if ( ! isbool(scale) || ! isscalar(scale))
+      error ("polyfit: scale must be a bool scalar");
+    endif
   endif
 
   if (! (isvector (x) && isvector (y) && size_equal (x, y)))
@@ -77,8 +91,16 @@ function [p, s, mu] = polyfit (x, y, n)
   l = length (x);
   x = reshape (x, l, 1);
   y = reshape (y, l, 1);
+  ## FIXME: is reshape() faster than (:) ?
+
+  ## scale the x values
+  if (scale)
+    xmax = max(abs(x));
+    x = x/xmax;
+  endif
 
   X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0));
+  # FIXME: maybe use repmat() for the ones() part, but this here seems faster
 
   p = X \ y;
 
@@ -93,16 +115,22 @@ function [p, s, mu] = polyfit (x, y, n)
     endif
 
     [s.R, dummy] = chol (X'*X);
+    ## FIXME: what if dummy is !0 ?
     s.X = X;
     s.df = l - n - 1;
     s.normr = norm (yf - y);
 
   endif
 
-  ## Return value should be a row vector.
+  ## Scale the polynomial back and return a row vector
+  if (scale)
+    power = length(p)-1:-1:0;
+    power = (xmax .^ power)(:);
+    p = (p ./ power);
+  endif
 
-  p = p.';
-
+  ## Return a row vector
+  p = p(:)';
 endfunction
 
 %!test
@@ -123,3 +151,20 @@ endfunction
 %! x = [-2, -1, 0, 1, 2];
 %! fail("polyfit (x, x.^2+x+1, [])");
 
+%% test difficult case where scaling is really needed
+%% checks also usage of second output argument
+%!shared s1, s2
+%!  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 ];
+%!  [p1,s1] = polyfit(x,y,10);
+%!  [p2,s2] = polyfit(x,y,10, true);
+%!
+%!  # fails without scaling (s1.norm ~6.33)
+%!error assert(s1.normr, 0.08,  0.01);
+%!
+%!  # but works with scaling
+%!assert(s2.normr, 0.08,  0.01);
+