/* 5.1. DEFINITIONS, TYPES, DECLARATIONS */
/* ================================ */

#define degree_plus_one 56
#define number_of_calculated_column_vectors 6
#define radius_of_domain 3.5e+00
#define radius_of_ball 5.0e-09
#define scaling_factor 8.956e+00
#define approx_linear_coeff 6.989740395952895e-01

/* types and external variables */
#include "Section5.2"

/* reps-macros and rup */
#include "Section5.8"

/* More declarations */
/* ----------------- */
main()
{
  extern int
    n,n1;

  extern reps
    rrho,rup(),rmaxabs();

  extern scalar
    srho;

  extern vector
    varg,vfact;

  int
    j,jmax;

  reps
    repsilon,rbeta,rtheta;

  scalar st1;

  vector
    v0,v0i,v0t1,v0t2;

/* Initialisations */
/* --------------- */
  n1=degree_plus_one;
  n=n1-1;
  read_approx_fixed_point(&v0);
  rrho=radius_of_domain;
  srconst(rrho,&srho);
  init_pack_and_unpack();
  init_multiplication_by(&vfact);
  init_composition_with(&varg);

/* The approximate fixed point */
/* --------------------------- */
/* calculate image */
  m(&v0,&v0i);

/* calculate norm of difference */
  printf("NORM OF DIFFERENCE :\n\n");
  vequal(&v0i,&v0t1);
  vminus(&v0,&v0t1);
  snorm(&v0t1,&st1);
  repsilon=rmaxabs(&st1);
  printf("epsilon %32.16e %30.16e\n\n",st1.l,st1.u);

/* Bound on the tangent map DM */
/* --------------------------- */
/* prepare ball around approx fixed point */
  rbeta=radius_of_ball;
  sscalar(rzero,rbeta,&v0.g);
  m(&v0,&v0t1);

/* norm of polynomial part */
  printf("BOUND ON DM :\n\n");
  jmax=number_of_calculated_column_vectors;
  rtheta=rzero;
  for(j=2;j<=jmax;j++){
    vsetzo(&v0t1);
    sequal(&sone,v0t1.p+j);
    dm(&v0t1,&v0t2);
    snorm(&v0t2,&st1);
    printf("column %d %31.16e %30.16e\n",j,st1.l,st1.u);
    rtheta=rmax2(rtheta,rmaxabs(&st1));
  }

/* higher order term */
  vsetzo(&v0t1);
  for(j=jmax+1;j<=n;j++)
    sscalar(rzero,rone,v0t1.p+j);
  sscalar(rzero,rone,&v0t1.h);
  dm(&v0t1,&v0t2);
  snorm(&v0t2,&st1);
  printf("higher order %27.16e %30.16e\n\n",st1.l,st1.u);
  rtheta=rmax2(rtheta,rmaxabs(&st1));

/* check if contraction */
 if(repsilon<rdcprod(rddiff(rone,rtheta),rbeta)){
    printf("***********************\n");
    printf("* contraction on ball *\n");
    printf("***********************\n");
  }
}

/* Include the subroutines */
/* ----------------------- */
/* initialisation routines */
#include "Section5.3"

/* the contraction M */
#include "Section5.4"

/* routines involving borel-functions */
#include "Section5.5"

/* vector operations */
#include "Section5.6"

/* routines acting on scalars */
#include "Section5.7"
/*---------------------------------------------*/