home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / STPLRG.C < prev    next >
Encoding:
Text File  |  1984-08-27  |  1.4 KB  |  69 lines
  1.    stplrg(m,x,y,n,c,iflag,a)
  2.  
  3.       /* this function will fit a polynomial of degree n    */
  4.       /* to m observations, by the method of least squares. */
  5.  
  6.       int *iflag,m,n;
  7.       float a[],c[],x[],y[];
  8.  
  9.    {
  10.       int i,ij,in1,n1,n2,nc,nnpn,n1n1,n1n2,j,l,k,jnp2,im1jp1,
  11.           iiflag;
  12.       float xpwr;
  13.  
  14.       n1 = n + 1;
  15.       n2 = n1 + 1;
  16.       nnpn = n-1 + n * n2;
  17.       n1n1 = n + n * n2;
  18.       n1n2 = n1 * n2 -1;
  19.  
  20.       ij = (n * n2);
  21.       for(j = 0; j <= ij; j+=n2)
  22.       {
  23.        for(i = 0; i <= n1; i++)
  24.        {
  25.         k = j + i;
  26.         a[k] = 0.;
  27.        }
  28.       }
  29.  
  30.       for(k = 1; k <= m; k++)
  31.       {
  32.        xpwr = x[k-1];
  33.        for(j = 1; j <= n; j++)
  34.        {
  35.         jnp2 = n1 + j * n2;
  36.         a[j] = a[j] + xpwr;
  37.         a[jnp2] = a[jnp2] + xpwr * y[k-1];
  38.         xpwr = xpwr * x[k-1];
  39.        }
  40.        for(i = 1; i <= n; i++)
  41.        {
  42.         in1 = n + i * n2;
  43.         a[in1] = a[in1] + xpwr;
  44.         xpwr = xpwr * x[k-1];
  45.        }
  46.        a[n1] = a[n1] + y[k-1];
  47.       }
  48.       a[0] = m;
  49.       for(i = n2; i <= n*n2; i+=n2)
  50.       {
  51.        ij = i;
  52.        for(j = 0; j <= n-1; j++)
  53.        {
  54.         im1jp1 = ij + 1 - n2;
  55.         a[ij] = a[im1jp1];
  56.         ij = ij + 1;
  57.        }
  58.       }
  59.       nc = -n2;
  60.       matinv(&iiflag,&iiflag,n1,nc,a,n2,i,xpwr);
  61.       *iflag = iiflag;
  62.       k = n1n2;
  63.       for(i = 0; i <= n; i++)
  64.       {
  65.        c[i] = a[k];
  66.        k = k - n2;
  67.       }
  68.    }
  69.