#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include "nr.h"
using namespace std;

#define M 25
#define KMAX 200

// global variables
double mu=0.8;
double eps=1.e-5;

double recursion( double y1, double y0, double c );

int main()
{
  int i, j, k, nn[M];
  double c[M], cg[M];
  double y00[M], y11[M], y22[M], y[M];
  double norm;
  
  // put random vectors (as normal distributed)
  for(i=0;i<M;i++) {
     nn[i] = rand( ) %100;
     y00[i] = NR::gasdev(nn[i]);
     y11[i] = NR::gasdev(nn[i]);
     c[i] = 1.0;
     cg[i] = c[i]/(1.0-mu);  
  }

  k = 0;
  do {
     k++;

     for (j=0; j<M; j++)
         y22[j] = recursion( y11[j], y00[j], c[j] );
     norm = 0.0; 
     for(i=0;i<M;i++) {
       norm = norm + abs(y22[i]-cg[i]);
     }

     for(i=0;i<M;i++) {
       y00[i] = y11[i];
       y11[i] = y22[i];
     }
   }
   while ( norm > eps && k < KMAX);
   
   cout << "#         iteration analysis" << endl;
   cout << "# dim M: " << M << ", mu= " << mu << ", eps = " << eps << endl;
   cout << "# norm:  " << norm << ", # no. of its k = " << k << endl << endl;
   cout << "#    result (vector)      limit " << endl;
   cout << "#============================================= " << endl;
   for(j=0;j<M;j++) { 
     printf("\t %13.10f \t %6.4f\n",  y22[j], cg[j] ); 
   }
   return 0;
}

// ========================================
 
double recursion( double y1, double y0, double c )
{
     return 2.0*mu*y1 - mu*y0 + c;
}
// =========================================