Though I was partially responsible for all this discussions,
I am now rethinking this approach.
polyfit should be as stable numerically as we can make it,
but I do not think it should try to make a "better"
fit than the one it was asked for.
Polyfit does fit to a*x^n+b*x^(n-1)... polynomial.
If people want it to fit to a*(x-x0)^n+b*(x-x0)^(n-1)
they should pre-condition data accordingly prior to
passing it to polyfit.
So to make long story short -- I think the original
Thomas' suggestion (just to normalize the data to max(x))
was the good one.
On related note -- it appears that wpolifit does a
better job than polyfit on un-conditioned data as well.
It appears that the difference is in polyfit
uses straight "\", while the wpolyfit does QR decomposition.
Perhaps it is worth to port this to polyfit as well.
Frankly I do not see why would not we copy wpolyfit to
octave as polyfit (with an additional change to scale x
unconditionally).