More robust polyfit

Hello, during some work I found the results of polyfit unnecessarily
inaccurate.

I have found out, that lsq system matrix is constructed via x.^n
(which supposedly walks through log/exp pow), and also the A\y
inversion did not work to expectations, compared to other numerical
tools.

Below, see link to a test data file.
x=readdlm("polyfit_8th_order_test.dat"); p=polyfit(x[:,1], x[:,2],
8); plot(x[:,2]-p(x[:,1])
-- a huge difference between previous (meaningless result) and
patched version. (And yes, the 8th order is physical in this data
set ;-))
https://github.com/Keno/Polynomials.jl/files/659739/polyfit_8th_order_test.dat.zip

Author: Marek P

src/Polynomials.jl

@@ -541,13 +541,29 @@ polyfit(xs, ys, 2)
 ```
 
 Original by [ggggggggg](https://github.com/Keno/Polynomials.jl/issues/19)
+More robust version by Marek Peca <[email protected]>
+    (1. no exponentiation in system matrix construction, 2. QR least squares)
 """
 function polyfit(x, y, n::Int=length(x)-1, sym::Symbol=:x)
     length(x) == length(y) || throw(DomainError)
     1 <= n <= length(x) - 1 || throw(DomainError)
 
-    A = [ float(x[i])^p for i = 1:length(x), p = 0:n ]
-    Poly(A \ y, sym)
+    #
+    # here unsure, whether zeros(Float64), or similar(x,...)
+    # however similar may yield unwanted surprise in case of e.g. x::Int
+    #
+    A=zeros(Float64, length(x), n+1)
+    #
+    # TODO: add support for poly coef bitmap
+    # (i.e. polynomial with some terms fixed to zero)
+    #
+    A[:,1]=1
+    for i=1:n
+        A[:,i+1]=A[:,i] .* x   # cumulative product more precise than x.^n
+    end
+    Aqr=qrfact(A)   # returns QR object, not a matrix
+    p=Aqr\y         # least squares solution via QR
+    Poly(p, sym)
 end
 polyfit(x,y,sym::Symbol) = polyfit(x,y,length(x)-1, sym)