// polynomial least square fit

#include <cmath>
#include <iostream>
#include "band.h"

using namespace std;

double power(double x, int k)
{  double t = 1.0;
   for (int i = 0; i < k; ++i)
   t *= x;
   return t;
}

double poly(const vector<double>& c, double x )
{// evaluate the poly c[0] + c[1]*pow(x,1) + ... + c[n-1]*pow(x,n-1)
  // using Horner's method
  int i, n=c.size();
  double val=0.0;

  for( i=n-1; i>=0; i-- )
    val = c[i] + x*val;

  return val;
}

double msqdev(const vector<double> &x, const vector<double> &y, const vector<double> &a)
{   int  n = x.size();
    double sqrtdev = 0;
    for (int i = 0;  i < n; ++i)
      {    double  diff = y[i] - poly(a, x[i]);
      sqrtdev  += diff * diff;
      }
    return sqrtdev/n;
}


double rootmeansquare(const vector<double> &x, const vector<double> &y, vector<double> &v)
{   int s = v.size();
    int n = x.size();

    bandMatrix X(s,s-1); // full matrix
    for (int j = 0; j < s; ++j)
      {  for (int k = j; k < s; ++k)
	  { double Xjk = 0;
	    for (int i = 0; i < n; ++i)
	      Xjk +=power(x[i], j+k);
            X(j,k) = Xjk;
	    if (k > j) X(k,j) = Xjk;
	}
      }
    for (int k = 0; k < s; ++k)
	{   double bk = 0;
	for (int i = 0; i < n; ++i)
	  bk += y[i] * power(x[i], k);
        v[k] = bk;
	}
     v = solve(X, v);
     double msd = msqdev(x, y, v);
     return sqrt(msd);
}

	
int main()
{
  int m, j;
  int numout = 201;

  // read poly degree
  cin >> m;
 
  vector<double> tin(11);
  vector<double> xin(11);
  vector<double> yin(11);
  double t, dt = 100./(numout-1), val;
  
  for(j = 0; j < 11; j ++)
    { 
      tin[j] = 1900. + 10*j;
      xin[j] = (tin[j] -1950.)/50.;
    }

  yin[0] = 76;
  yin[1] = 91.97;
  yin[2] = 105.71;
  yin[3] = 123.2;
  yin[4] = 131.67;
  yin[5] = 150.7;
  yin[6] = 179.32;
  yin[7] = 203.21;
  yin[8] = 226.51;
  yin[9] = 249.63;
  yin[10] = 281.42;


  vector<double> b(m+1);
  double rms = rootmeansquare(xin, yin, b);
  cout << "root mean square deviation is  " << rms << endl;
  for (j = 0; j <= m; j ++)
    cout << b[j] << endl;
  cout<<"the fitted population at 2006 is " << poly(b,(2006-1950.)/50.)<<endl;
  cout<<"the fitted population at 2010 will be " << poly(b,(2010-1950.)/50.)<<endl;

  for(j = 0, t = 1900; j<numout; t+=dt, j++){
    // evaluate interpolant at t
    val = poly(b,(t-1950.)/50.);
    cout << t <<' ' << val << endl;
  }

  
  return 0;
}