code: def fitspoly(self, n, x, y, sd=None, wt=1.0, fid=0):
          
          a = max(abs(max(x)), abs(min(x)))
          if a == 0:
              a = 1
          a = 1.0 / a
          b = NUM.power(a, range(n + 1))
          if self.set(n=n + 1, fid=fid):
              self.linear(poly(n), x * a, y, sd, wt, fid)
              self._fitids[fid]["sol"] *= b
              self._fitids[fid]["error"] *= b
              return self.linear(poly(n), x, y, sd, wt, fid)
docs: Create normal equations from the specified condition equations, and
          solve the resulting normal equations. It is in essence a combination.
  
          The method expects that the properties of the fitter to be used have
          been initialized or set (like the number of simultaneous solutions m
          the type; factors). The main reason is to limit the number of
          parameters on the one hand, and on the other hand not to depend
          on the actual array structure to get the variables and type. Before
          fitting the x-range is normalized to values less than 1 to cater for
          large difference in x raised to large powers. Later a shift to make x
          around zero will be added as well.
  
          :param n: the order of the polynomial to solve for
          :param x: the abscissa values
          :param y: the ordinate values
          :param sd: standard deviation of equations (one or more values used
                     cyclically)
          :param wt: an optional alternate for `sd`
          :param fid: the id of the sub-fitter (numerical)
  
          Example::
  
              fit = fitserver()
              x = N.arange(1,11) # we have values at 10 'x' values
              y = 2. + 0.5*x - 0.1*x**2 # which are 2 +0.5x -0.1x^2
              fit.fitspoly(3, x, y) # fit a 3-degree polynomial
              print fit.solution(), fit.error() #  show solution and their errors