// evaluate a monomial basis interpolant at  many points

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

using namespace std;

typedef double (*FP)(double); // an FP can point to f below.

double f( double x ); // function to be interpolated

void getMonCoef( const vector<double>& y,  vector<double>& c, FP f );

double poly( const vector<double>& c, double x );

double f( double x )
{ return 1.0/(1+25*x*x);}

void getMonCoef( const vector<double>& y,  vector<double>& c, FP f )
{ // return c so that at the points x in y we have
  // c[0] + x*(c[1] + x *(... )) = f(x)
  int n = y.size();
  bandMatrix A(n,n-1); // full matrix
  int i,j;

  // form the Vandermonde matrix
  for( i=0; i<n; i++ )
    A(i,0) = 1.0;
  for( i=0; i<n; i++ )
    for( j=1; j<n; j++ )
      A(i,j) = A(i,j-1)*y[i]; // A(i,j) = pow(y[i],j)

  // for vector of function values
  for( i=0; i<n; i++ )
    c[i] = f(y[i]);

  c = solve( A, c );
  
}

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;
}


 int main()
{
  int n,i,j;
  int numout = 201;

  // read poly degree
  cin >> n;

  vector<double> y(n+1);
  double x, dx = 2.0/(numout-1), h =2.0/n, val;
  double pi = 4*atan(1.0);
  double ph = pi/n;

  // generate the n+1 interpolation points
  for(i=0; i<=n; i++ )
     y[i] = -1.0 + i*h;

  // get monomial basis coefficients
  vector<double> c(n+1);
  getMonCoef(y,c,f);

  for( j=0, x=-1.0; j<numout; x+=dx, j++){
    // evaluate interpolant at x
    val = poly(c,x);
    cout << x << ' ' << f(x) << ' ' << val << ' ' << f(x) - val << endl;
  }


  return 0;
}