I am using a least square polynomial approximation to find smooth curves approximating data coming from image analysis. More specifically I've blobs of pixels that I know are a reasonably smooth almost horizontal slightly fat line in an image and I want to find the best mathematical curve $y=p(x)$ where $p$ is a polynomial that minimizes the total "error" defined as
$$ E = \sum_i\left(y_i - p(x_i)\right)^2 $$
The solution is found by setting to zero the partial derivatives $\partial E \over a_i$ for all coefficients $a_i$ of the polynomial. This produces a linear system of $n+1$ equations in $n+1$ unknowns for a degree $n$ polynomial that is then solved explicitly with a Gauss-Jordan solver.
The code works and I get reasonable results when using low-degree polynomials (up to 5). What I find strange is that the numerical stability increases a lot (no problems even with degree 15) if I do a very simple shifting and scaling transformation:
$$ x' = 2{x - \bar{x} \over x_{\rm max} - x_{\rm min}} $$
Transforming also $y$ doesn't seem to have any effect (and thinking to how $y$ values are used that's not a big surprise).
Does this happen because my solver is bad (I'm using a Gauss-Jordan solver with pivoting I wrote myself so that's not to exclude) or it's a intrinsic problem of floating point math?
The image is about 4000 pixels wide so un-scaled data computation has to deal with numbers rather large (e.g. $4000^{30} \approx 10^{108}$ for a degree 15 polynomial) but I thought that FP was designed exactly to handle this case and a double-precision value can handle up to $10^{300}$ (I also found no mention of this scaling/shifting issue and problems begin to appear even at just degree 7).
Sorry if this is a trivial question but I know very little about numerical analysis.