5. The Computer Program

/∗ 5.1. DEFINITIONS, TYPES, DECLARATIONS ∗/  printf(”BOUND ON DM :\n\n”);
/∗ ================================ ∗/  jmax=number of calculated column vectors;

rtheta=rzero;

#define degree plus one 56  for(j=2;j<=jmax;j++){
#define number of calculated column vectors 6  vsetzo(&v0t1);
#define radius of domain 3.5e+00  sequal(&sone,v0t1.p+j);
#define radius of ball 5.0e−09  dm(&v0t1,&v0t2);
#define scaling factor 8.956e+00  snorm(&v0t2,&st1);
#define approx linear coeff 6.989740395952895e−01  printf(”column %d %31.16e %30.16e\n”,j,st1.l,st1.u);

rtheta=rmax2(rtheta,rmaxabs(&st1));

/∗ types and external variables ∗/  }

#include ”Section5.2”

/∗ higher order term ∗/

/∗ reps−macros and rup ∗/  vsetzo(&v0t1);
#include ”Section5.8”  for(j=jmax+1;j<=n;j++)

sscalar(rzero,rone,v0t1.p+j);

/∗ More declarations ∗/  sscalar(rzero,rone,&v0t1.h);

dm(&v0t1,&v0t2);

main()  snorm(&v0t2,&st1);
{  printf(”higher order %27.16e %30.16e\n\n”,st1.l,st1.u);
extern int  rtheta=rmax2(rtheta,rmaxabs(&st1));

n,n1;

/∗ check if contraction ∗/

extern reps  if(repsilon<rdcprod(rddiff(rone,rtheta),rbeta)){
rrho,rup(),rmaxabs();  printf(”∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗\n”);

printf(”∗ contraction on ball ∗\n”);

extern scalar  printf(”∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗\n”);
srho;  }

}

extern vector

varg,vfact;  /∗ Include the subroutines ∗/
int  /∗ initialisation routines ∗/
j,jmax;  #include ”Section5.3”
reps  /∗ the contraction M ∗/
repsilon,rbeta,rtheta;  #include ”Section5.4”
scalar st1;  /∗ routines involving borel−functions ∗/

#include ”Section5.5”

vector

v0,v0i,v0t1,v0t2;  /∗ vector operations ∗/

#include ”Section5.6”

/∗ Initialisations ∗/

/∗ routines acting on scalars ∗/
n1=degree plus one;  #include ”Section5.7”

n=n1−1; read approx fixed point(&v0);

rrho=radius of domain;  /∗ 5.2. TYPES AND EXTERNAL VARIABLES ∗/
srconst(rrho,&srho);  /∗ ============================== ∗/

init pack and unpack();

init multiplication by(&vfact);  #include <stdio.h>

init composition with(&varg);

typedef int logical;

/∗ The approximate fixed point ∗/  typedef double reps;
/∗ calculate image ∗/  typedef struct s{
m(&v0,&v0i);  reps l;

reps u;

/∗ calculate norm of difference ∗/  }
printf(”NORM OF DIFFERENCE :\n\n”);  scalar;

vequal(&v0i,&v0t1);

vminus(&v0,&v0t1);  typedef struct v{
snorm(&v0t1,&st1);  scalar p[degree plus one];
repsilon=rmaxabs(&st1);  scalar g;
printf(”epsilon %32.16e %30.16e\n\n”,st1.l,st1.u);  scalar h;

}

/∗ Bound on the tangent map DM ∗/  vector;
/∗ prepare ball around approx fixed point ∗/  typedef struct b{
rbeta=radius of ball;  scalar p[degree plus one];
sscalar(rzero,rbeta,&v0.g);  scalar g;
m(&v0,&v0t1);  scalar h;

reps e;

/∗ norm of polynomial part ∗/  }

borel function;  siconst(2,&s4);

sinv(&s2,&si2);

int n,n1;  sinv(&s3,&si3);

sinv(&s4,&si4);

reps  }

rrho,r0save,r1save,r2save,r3save,

rzero =0.0e+00,  init multiplication by(pv)
rone =1.0e+00,  vector ∗pv;
rtwo =2.0e+00,  {
rthree=3.0e+00,  extern scalar
rfour =4.0e+00,  srho;

rhalf =0.5e+00,

r3half=1.5e+00;  scalar s;
scalar  /∗ produce v(z)=(1+z∗rho/33)∗∗(−1/2) ∗/
sa,srho,sdnu,srarg,snormal,  siconst(33,&s);
s2,s3,s4,si2,si3,si4,  squot(&srho,&s,&s);
szero ={0.0e+00, 0.0e+00},  vsetms(pv,&s);
sone ={1.0e+00, 1.0e+00},  }

stwo ={2.0e+00, 2.0e+00},

sthree={3.0e+00, 3.0e+00},  init composition with(pv)
sfour ={4.0e+00, 4.0e+00},  vector ∗pv;
shalf ={0.5e+00, 0.5e+00},  {
s3half={1.5e+00, 1.5e+00};  extern reps

r0save,r1save,r2save,r3save,

vector  rmaxabs(),rnorm();

ve,varg,vfact,vhhat,vpower[degree plus one];

extern scalar

borel function bhto7;  srho,srarg;

extern vector

/∗ 5.3. INITIALISATION ROUTINES ∗/  vpower[];

/∗ ======================= ∗/

int

/∗ read approx fixed point ∗/  i,j,k,ii,kh,nh;

/∗ init pack and unpack ∗/

/∗ init multiplication by  ∗/  reps
/∗ init composition with  ∗/  r1,r2,r3;

scalar

read approx fixed point(pv)  st1,st2;

vector ∗pv;

{  vector
extern FILE  vt1,vt2,vt3,vz,vs[6];

∗fopen();

srconst(rone/scaling factor,&srarg);

FILE ∗fp;

/∗ produce the vector ∗/

int i;  /∗ v(z)=((3/11+z∗rho/22)/(1+z∗rho/33))/(r∗rho) ∗/

siconst(33,&st1);

double d;  squot(&srho,&st1,&st1);

vsetin(&vz,&st1);

fp=fopen(”Section5.9”,”r”);  siconst(11,&st1);
fscanf(fp,”%∗86c”);  sprod(&st1,&srho,&st1);
for(i=0;i<=n;i++){  squot(&sthree,&st1,&st1);
fscanf(fp,”%lf”,&d);  siconst(22,&st2);
srconst(d,pv−>p+i);  sinv(&st2,&st2);
}  vsetli(&vt1,&st1,&st2);
fclose(fp);  vmult(&vt1,&vz,&vz);

sinv(&srarg,&st1);

sequal(&szero,&pv−>g);  vsmult(&vz,&st1,&vz);
sequal(&szero,&pv−>h);  vequal(&vz,pv);

}

/∗ save powers of vector v ∗/

init pack and unpack()  vmult(&vz,&vz,&vt2);
{  j=2;
extern scalar  i=1;
sa,s2,s3,s4,si2,si3,si4;  while(j<=n1){

vequal(&vt2,&vs[i++]);

extern vector  vmult(&vt2,&vt2,&vt2);
ve;  j=j+j;

}

/∗ make approximate eigenvector ve ∗/  ii=1;
vsetzo(&ve);  vequal(&vz,&vpower[ii++]);
srconst(rone,ve.p+1);  r1=rmax2(rone,rnorm(0,&vz));
srconst(.6e+00,ve.p+2);  r2=rnorm(1,&vz);

for(i=3;i<=n1;i+=2){

/∗ choose coordinates ∗/  k=(i−1)/2;
srconst(approx linear coeff,&sa);  j=1;
siconst(1,&s2);  vsetli(&vt2,&sone,&szero);
siconst(3,&s3);  while(k>0){

kh=k/2;

if(kh+kh!=k){  extern borel function
vequal(&vs[j],&vt3);  bhto7;

vmult(&vt2,&vt3,&vt2);

}  scalar
k=kh;  st1,st2;

++j;

}  vector
vequal(&vt2,&vpower[ii++]);  vt1,vt2;

r1=rmax2(r1,rnorm(0,&vt2));

r2=rmax2(r2,rnorm(1,&vt2));  borel function
if(ii<=n){  bt1,bt2,bt4,bt8;

vmult(&vt2,&vz,&vt3);

vequal(&vt3,&vpower[ii++]);  /∗ take borel transform ∗/
r1=rmax2(r1,rnorm(0,&vt3));  borel(pv,&bt1);

r2=rmax2(r2,rnorm(1,&vt3));

}  /∗ take 8th power of bh ∗/
}  srconst(.02e+00,&st1);
nh=n/2;  bmult(&bt1,&bt1,&st1,&bt2);
if(nh+nh==n)  srconst(.04e+00,&st1);
r3=rnorm(0,&vt3);  bmult(&bt2,&bt2,&st1,&bt4);
else  srconst(.04e+00,&st1);
r3=rnorm(0,&vt2);  bmult(&bt4,&bt4,&st1,&bt8);

if(r3>rone)

printf(”error in vcomp, norm = %30.16e\n”,r3);  /∗ take inverse borel transform ∗/

iborel(&bt8,&vt1);

/∗ save constants ∗/

r0save=rmaxabs(vz.p);  /∗ substitute varg in vh hat ∗/
r1save=r1;  vargcmp(&vt1,&vt1);

r2save=r2;

r3save=r3;  /∗ multiply with vfact ∗/
}  vmult(&vt1,&vfact,&vt1);

/∗ normalize to one ∗/

/∗ 5.4. THE CONTRACTION M ∗/  sinv(vt1.p,&snormal);
/∗ ==================== ∗/  vsmult(&vt1,&snormal,&vhhat);

sequal(&sone,vhhat.p);

/∗ m  : the contraction  ∗/  vequal(&vhhat,pvi);
/∗ dm  : the tangent map of m  ∗/
/∗ r  : the RG transformation  ∗/  /∗ save 7th power of bh ∗/
/∗ dr  : the tangent map of r  ∗/  srconst(.03e+00,&st1);
/∗ unpack : unpack operator used in m  ∗/  bmult(&bt4,&bt2,&st1,&bhto7);
/∗ dunpack : unpack operator used in dm  ∗/  srconst(.03e+00,&st1);
/∗ wzero : part of unpack and dunpack  ∗/  bmult(&bhto7,&bt1,&st1,&bhto7);
/∗ pack  : the left inverse of unpack  ∗/  siconst(8,&st1);
/∗ u param : computes parameter for unpack ∗/  bsmult(&bhto7,&st1,&bhto7);
/∗ w param : computes parameter for dunpack ∗/  }

dr(pv,pvi)

m(pv0,pv0i)  vector ∗pv,∗pvi;
vector ∗pv0,∗pv0i;  {
{  extern scalar
scalar snu;  snormal;
vector vt1;  extern vector

vhhat;

unpack(pv0,&snu,&vt1);

r(&vt1,&vt1);  extern borel function
pack(&snu,&vt1,pv0i);  bhto7;

}

scalar st1;

dm(pv0,pv0i)

vector ∗pv0,∗pv0i;  vector
{  vt1,vt2;

scalar snu0;

borel function bt1;

vector vt1;

/∗ take borel transform ∗/
dunpack(pv0,&snu0,&vt1);  borel(pv,&bt1);

dr(&vt1,&vt1);

pack(&snu0,&vt1,pv0i);  /∗ multiply with 7th power of bh ∗/
}  srconst(.04e+00,&st1);

bmult(&bt1,&bhto7,&st1,&bt1);

r(pv,pvi)

vector ∗pv,∗pvi;  /∗ take inverse borel transform ∗/
{  iborel(&bt1,&vt1);

extern scalar

snormal;  /∗ substitute varg in vh hat ∗/

vargcmp(&vt1,&vt1);

extern vector

vhhat;  /∗ multiply with vfact ∗/

vmult(&vt1,&vfact,&vt1);

sequal(&szero,pv0−>p);

/∗ renormalize ∗/  sequal(&szero,pv0−>p+1);
vsmult(&vhhat,vt1.p,&vt2);  sprod(vt1.p+2,&s2,pv0−>p+2);
vminus(&vt2,&vt1);  sprod(vt1.p+3,&s3,pv0−>p+3);
sequal(&szero,vt1.p);  sprod(vt1.p+4,&s4,pv0−>p+4);

vsmult(&vt1,&snormal,&vt1);

vequal(&vt1,pvi);  for(i=5;i<=n;i++)
}  sequal(vt1.p+i,pv0−>p+i);
unpack(pv0,ps,pv)  ru=rmax2(rmax3(s2.u,s3.u,s4.u),rone);
vector ∗pv0,∗pv;  rl=rmin2(rmin3(s2.l,s3.l,s4.l),rone);
scalar ∗ps;  sscalar(rl,ru,&st1);

{

vector vt1;  sprod(&vt1.g,&st1,&pv0−>g);

sequal(&vt1.h,&pv0−>h);

wzero(pv0,pv);  }

sequal(&sone,pv−>p);

sequal(&sa,pv−>p+1);  u param(pv,ps)
u param(pv,ps);  vector ∗pv;
vsmult(&ve,ps,&vt1);  scalar ∗ps;
vadd(&vt1,pv);  {
}  extern logical

lsequal();

dunpack(pv0,ps,pv)

vector ∗pv0,∗pv;  extern reps
scalar ∗ps;  rmaxabs();

{

vector vt1;  extern scalar

sdnu;

wzero(pv0,pv);

w param(pv,ps);  int i;

vsmult(&ve,ps,&vt1);

vadd(&vt1,pv);  reps
}  rt1,rt2,rrad,rdnu;
wzero(pv0,pv)  scalar
vector ∗pv0;  srad,seps,st1,st2,st3;

vector ∗pv;

{  vector
int i;  vt1,vt2,vdnu;
reps  sequal(&szero,ps);

rl,ru;

/∗ calculate precision of center ∗/

scalar st1;  vequal(pv,&vt1);

r(&vt1,&vt2);

sequal(&szero,pv−>p);  rt1=rmaxabs(&vt1.g);
sequal(&szero,pv−>p+1);  rt2=rmaxabs(&vt2.g);
sprod(pv0−>p+2,&si2,pv−>p+2);  sscalar(−rt1,rt1,&st1);
sprod(pv0−>p+3,&si3,pv−>p+3);  sscalar(−rt2,rt2,&st2);
sprod(pv0−>p+4,&si4,pv−>p+4);  sadd(vt1.p+1,&st1);

sadd(vt2.p+1,&st2);

for(i=5;i<=n;i++)  sequal(&st2,&seps);
sequal(pv0−>p+i,pv−>p+i);  sminus(&st1,&seps);
ru=rmax2(rmax3(si2.u,si3.u,si4.u),rone);  /∗ estimate derivative on neighborhood ∗/
rl=rmin2(rmin3(si2.l,si3.l,si4.l),rone);  rrad=5.e−3;
sscalar(rl,ru,&st1);  sscalar(−rrad,rrad,&srad);
sprod(&pv0−>g,&st1,&pv−>g);  for(i=1;i<=2;i++){
sequal(&pv0−>h,&pv−>h);  sequal(&srad,&st1);
}  sadd(ps,&st1);

vequal(pv,&vt1);

pack(ps,pv,pv0)  vsmult(&ve,&st1,&vt2);
scalar ∗ps;  vadd(&vt2,&vt1);
vector ∗pv,∗pv0;  r(&vt1,&vt2);
{  dr(&ve,&vdnu);
int i;  rt2=rmaxabs(&vdnu.g);

sscalar(−rt2,rt2,&sdnu);

reps  sadd(vdnu.p+1,&sdnu);
rl,ru;  sequal(&sdnu,&st2);

sminus(&sone,&st2);

scalar st1;  squot(&seps,&st2,&st2);

sneg(&st2,&st2);

vector  sinter(&srad,&st2,&st3);
vt1,vt2;  if(!lsequal(&st3,&st2)){

printf(”error in u param, no contraction. step %d\n”,i);

vequal(pv,&vt1);  sequal(&st1,ps);
vsmult(&ve,ps,&vt2);  return;
vminus(&vt2,&vt1);  }

sequal(&st2,&srad);  sprod(&pb−>h,&s,&pbres−>h);
}  pbres−>e=pb−>e;
sadd(&srad,ps);  }

}

bmult(pb1,pb2,pseps,pbres)

w param(pv,ps)  borel function ∗pb1,∗pb2,∗pbres;
vector ∗pv;  scalar ∗pseps;
scalar ∗ps;  {
{  extern logical
extern reps  lsequal();

rmaxabs();

extern reps

extern scalar  rmaxabs();

sdnu;

int

reps rt1;  i,k;
scalar st1;  scalar

st1,st2,s1,s2,

vector vt1;  sb1,sb2,sb1eps,sb2eps,

sa1g,sa1h,sa2g,sa2h,

dr(pv,&vt1);  sagres,sbres;
rt1=rmaxabs(&vt1.g);  borel function
sscalar(−rt1,rt1,&st1);  bt1;

sadd(vt1.p+1,&st1);

sequal(&sone,ps);  sequal(&pb1−>g,&sa1g);
sminus(&sdnu,ps);  sequal(&pb1−>h,&sa1h);
squot(&st1,ps,ps);  sequal(&pb2−>g,&sa2g);
}  sequal(&pb2−>h,&sa2h);

srconst(pb1−>e,&sb1); srconst(pb2−>e,&sb2);

/∗ 5.5. ROUTINES INVOLVING BOREL FUNCTIONS ∗/

/∗ =================================== ∗/  /∗ compute polynomial part ∗/

bsetzo(&bt1);

/∗ bsetzo : set borel function equal zero  ∗/  for(k=0;k<=n;k++){
/∗ bequal : copy borel function  ∗/  sequal(&szero,&s1);
/∗ bsmult : scalar times borel function  ∗/  for(i=0;i<=k;i++){
/∗ bmult  : product of two borel functions  ∗/  sprod(pb1−>p+i,pb2−>p+k−i,&s2);
/∗ sabound : bound used in bmult  ∗/  sadd(&s2,&s1);
/∗ borel  : the borel transform  ∗/  }
/∗ iborel  : the inverse borel transform  ∗/  sequal(&s1,bt1.p+k);

}

bsetzo(pb)  /∗ determine constant in exponential ∗/
borel function ∗pb;  sequal(pseps,&sb1eps);
{  sadd(&sb1,&sb1eps);
int i;  sequal(pseps,&sb2eps);

sadd(&sb2,&sb2eps);

for(i=0;i<=n;i++)  sequal(pseps,&sbres);
sequal(&szero,pb−>p+i);  sadd(&sb1,&sbres);
sequal(&szero,&pb−>g);  sadd(&sb2,&sbres);

sequal(&szero,&pb−>h);

pb−>e=rzero;  /∗ bound error terms ∗/
}  sequal(&szero,&sagres);
bequal(pb,pbres)  for(k=n1;k<=n+n;k++){
borel function ∗pb,∗pbres;  sequal(&szero,&st1);
{  for(i=k−n;i<=n;i++){
int i;  sprod(pb1−>p+i,pb2−>p+k−i,&st2);

sadd(&st2,&st1);

for(i=0;i<=n;i++)  }
sequal(pb−>p+i,pbres−>p+i);  sabs(&st1,&st1);
sequal(&pb−>g,&pbres−>g);  sabound(&sbres,k−1,&st2);
sequal(&pb−>h,&pbres−>h);  sprod(&st1,&st2,&st1);
pbres−>e=pb−>e;  sadd(&st1,&sagres);
}  }
bsmult(pb,ps,pbres)  for(i=0;i<=n;i++){
borel function ∗pb,∗pbres;  sabs(pb1−>p+i,&st1);
scalar ∗ps;  sprod(&st1,&sa2g,&st1);
{  sabound(&sb1eps,i,&st2);
int i;  sprod(&st1,&st2,&st1);

sadd(&st1,&sagres);

scalar s;  sabs(pb1−>p+i,&st1);

sprod(&st1,&sa2h,&st1);

for(i=0;i<=n;i++)  sabound(&sb1eps,i+n,&st2);
sprod(ps,pb−>p+i,pbres−>p+i);  sprod(&st1,&st2,&st1);

sadd(&st1,&sagres);

sabs(ps,&s);  sabs(pb2−>p+i,&st1);
sprod(&pb−>g,&s,&pbres−>g);  sprod(&st1,&sa1g,&st1);

sabound(&sb2eps,i,&st2);

sprod(&st1,&st2,&st1);  int i;

sadd(&st1,&sagres);

sabs(pb2−>p+i,&st1);  reps
sprod(&st1,&sa1h,&st1);  rfive=5.e+00;

sabound(&sb2eps,i+n,&st2);

sprod(&st1,&st2,&st1);  scalar
sadd(&st1,&sagres);  st1,st2,st3,sag,sb,
}  srhohat,sfact,stdiff;
sabound(pseps,1,&st1);  sequal(&sone,&st1);
sprod(&sa1g,&sa2g,&st2);  sequal(&shalf,&st2);
sprod(&st1,&st2,&st1);  for(i=0;i<=n;i++){
sadd(&st1,&sagres);  sprod(pb−>p+i,&st1,pv−>p+i);

sprod(&st1,&st2,&st1);

if(!lsequal(&sa1h,&szero) &&  sadd(&sone,&st2);
!lsequal(&sa2h,&szero)){  }
sabound(pseps,n1+n,&st1);  sequal(&st1,&sfact);

sprod(&sa1h,&sa2h,&st2);

sprod(&st1,&st2,&st1);  /∗ extract exponential growth ∗/
sadd(&st1,&sagres);  /∗ and determine domains ∗/
}  srconst(pb−>e,&sb);

sprod(&sb,&srarg,&st1);

/∗ return result ∗/

sscalar(rzero,rmaxabs(&sagres),&bt1.g);  if(st1.u>rone)
srconst(rzero,&bt1.h);  printf(”domains in iborel incompatible %30.16\n”,st1.u);

bt1.e=rmaxabs(&sbres);

bequal(&bt1,pbres);  sprod(&s3half,&st1,&st1);
}  sadd(&sone,&st1);

squot(&st1,&sb,&st1);

sabound(pseps,k,ps)  srconst(rtwo/rfive,&st2);
scalar ∗pseps,∗ps;  sprod(&st1,&st2,&st1);
int k;  srconst(st1.l,&srhohat);

{

scalar  /∗ adapt polynomial part to this domain ∗/
st1,st2,sexp1;  sequal(&sone,&st1);
for(i=0;i<=n;i++){
if(k==0)  sprod(pv−>p+i,&st1,pv−>p+i);
sequal(&sone,ps);  sprod(&st1,&srarg,&st1);
else{  }

srconst(2.718281828e+00,&sexp1);

sprod(pseps,&sexp1,&st1);  /∗ calculate error terms on this domain ∗/
siconst(k,&st2);  srconst(pb−>g.u,&sag);
squot(&st2,&st1,&st1);  sprod(&sag,&shalf,&st1);
spower(&st1,k,&st1);  sprod(&sb,&srhohat,&st2);
srconst(st1.u,ps);  sequal(&sone,&stdiff);
}  sminus(&st2,&stdiff);
}  ssqrt(&stdiff,&stdiff);

spower(&stdiff,3,&st2);

borel(pv,pb)  squot(&srarg,&st2,&st2);
vector ∗pv;  sprod(&st1,&st2,&st1);
borel function ∗pb;  squot(&srarg,&srhohat,&st2);
{  sequal(&sone,&st3);
extern reps  sminus(&st2,&st3);
rmaxabs();  squot(&st1,&st3,&st1);

sscalar(rzero,rmaxabs(&st1),&pv−>g);

int i;  sequal(&szero,&pv−>h);

}

scalar st1,st2;

/∗ 5.6. OPERATIONS ON VECTORS ∗/

sequal(&sone,&st1);  /∗ ======================= ∗/

sequal(&shalf,&st2);

for(i=0;i<=n;i++){  /∗ vsetzo : set vector equal zero  ∗/
squot(pv−>p+i,&st1,pb−>p+i);  /∗ snorm : norm of vector  ∗/
sprod(&st1,&st2,&st1);  /∗ rnorm : bound partial norm of vector  ∗/
sadd(&sone,&st2);  /∗ vequal : copies vector  ∗/
}  /∗ vadd  : sum of vectors  ∗/
squot(&pv−>g,&shalf,&st2);  /∗ vminus : difference of vectors  ∗/
sscalar(rzero,rmaxabs(&st2),&pb−>g);  /∗ vsmult : product of vector and scalar  ∗/
squot(&pv−>h,&st1,&st2);  /∗ vmult : product of two vectors  ∗/
sscalar(rzero,rmaxabs(&st2),&pb−>h);  /∗ vargcmp : substitute varg in vector  ∗/
pb−>e=rone;  /∗ vsetli  : create vector v(z)=c+lz  ∗/
}  /∗ vsetin : create vector v(z)=1/(1+sz)  ∗/
/∗ vsetms : set vector v(z)=(1+sz)∗∗(−1/2)  ∗/

iborel(pb,pv) borel function ∗pb;

vector ∗pv;  vsetzo(pv)
{  vector ∗pv;
extern reps  {
rmaxabs();  int i;

for(i=0;i<=n;i++)
for(i=0;i<=n;i++)  sminus(pv−>p+i,pvres−>p+i);

sequal(&szero,pv−>p+i);

sequal(&szero,&pv−>g);  sequal(&pv−>g,&st1);
sequal(&szero,&pv−>h);  sminus(&pvres−>g,&st1);
}  sequal(&pv−>g,&st2);

sadd(&pvres−>g,&st2);

snorm(pv,psres)  sscalar(rminabs(&st1),rmaxabs(&st2),&pvres−>g);

vector ∗pv;

scalar ∗psres;  sequal(&pv−>h,&st1);
{  sminus(&pvres−>h,&st1);
int i;  sequal(&pv−>h,&st2);

sadd(&pvres−>h,&st2);

scalar  sscalar(rminabs(&st1),rmaxabs(&st2),&pvres−>h);
s,st1;  }
sequal(&szero,&s);  vsmult(pv,ps,pvres)
for(i=0;i<=n;i++){  vector ∗pv,∗pvres;
sabs(pv−>p+i,&st1);  scalar ∗ps;
sadd(&st1,&s);  {
}  int i;

sadd(&pv−>g,&s);

sadd(&pv−>h,&s);  scalar s;

sequal(&s,psres);

}  for(i=0;i<=n;i++)

sprod(pv−>p+i,ps,pvres−>p+i);

reps rnorm(imin,pv)

int imin;  sabs(ps,&s);
vector ∗pv;  sprod(&pv−>g,&s,&pvres−>g);
{  sprod(&pv−>h,&s,&pvres−>h);
extern reps  }

rmaxabs();

vmult(pva,pvb,pvres)

int i;  vector ∗pva,∗pvb,∗pvres;

{

scalar  int
s,st1;  i,k;
sequal(&szero,&s);  scalar
for(i=imin;i<=n;i++){  s1,s2,s3,st1,st2,
sabs(pv−>p+i,&st1);  seg,seh,sag,sah,sbg,sbh,
sadd(&st1,&s);  sx[degree plus one];

}

sadd(&pv−>g,&s);  sequal(&pva−>g,&sag);
sadd(&pv−>h,&s);  sequal(&pva−>h,&sah);
return(rmaxabs(&s));  sequal(&pvb−>g,&sbg);
}  sequal(&pvb−>h,&sbh);
vequal(pv,pvres)  /∗ calculate polynomial part ∗/
vector ∗pv,∗pvres;  for(i=0;i<=n;i++)
{  sequal(&szero,sx+i);

int i;

for(k=0;k<=n;k++){

for(i=0;i<=n;i++)  sequal(&szero,&s1);
sequal(pv−>p+i,pvres−>p+i);  for(i=0;i<=k;i++){
sequal(&pv−>g,&pvres−>g);  sprod(pva−>p+i,pvb−>p+k−i,&s2);
sequal(&pv−>h,&pvres−>h);  sadd(&s2,&s1);
}  }

sequal(&s1,sx+k);

vadd(pv,pvres)  }

vector ∗pv,∗pvres;

{  /∗ estimate higher order error ∗/
int i;  sequal(&szero,&seh);

for(k=n1;k<=n+n;k++){

for(i=0;i<=n;i++)  sequal(&szero,&s1);
sadd(pv−>p+i,pvres−>p+i);  for(i=k−n;i<=n;i++){
sadd(&pv−>g,&pvres−>g);  sprod(pva−>p+i,pvb−>p+k−i,&s2);
sadd(&pv−>h,&pvres−>h);  sadd(&s2,&s1);
}  }

sabs(&s1,&s1);

vminus(pv,pvres)  sadd(&s1,&seh);
vector ∗pv,∗pvres;  }
{  sequal(&szero,&s1);
extern reps  sequal(&szero,&s2);
rminabs(),rmaxabs();  for(i=1;i<=n;i++){

sabs(pva−>p+n−i,&s3);

int i;  sadd(&s3,&s1);

sabs(pvb−>p+n−i,&s3);

scalar  sadd(&s3,&s2);
st1,st2;  }

sabs(pva−>p+n,&st1);

sadd(&s1,&st1);  sscalar(rzero,rt1,&st1);
sprod(&sbh,&st1,&st1);  sadd(&st1,&vt1.g);

sabs(pvb−>p+n,&st2);

sadd(&s2,&st2);  vequal(&vt1,pvres);

sprod(&sah,&st2,&st2);

sadd(&st2,&st1);  }

sprod(&sah,&sbh,&st2);

sadd(&st2,&st1);  vsetli(pv,psc,psl)
sprod(&sah,&sbg,&st2);  vector ∗pv;
sadd(&st2,&st1);  scalar ∗psc,∗psl;
sprod(&sag,&sbh,&st2);  {
sadd(&st2,&st1);  vsetzo(pv);
sabs(pvb−>p+n,&st2);  sequal(psc,pv−>p);
sprod(&st2,&sag,&st2);  sequal(psl,pv−>p+1);
sadd(&st2,&st1);  }

sabs(pva−>p+n,&st2);

sprod(&st2,&sbg,&st2);  vsetin(pv,ps)
sadd(&st2,&st1);  vector ∗pv;
sadd(&st1,&seh);  scalar ∗ps;
{
/∗ estimate general error ∗/  extern reps
sprod(&s2,&sag,&st1);  rmaxabs();

sprod(&s1,&sbg,&st2);

sadd(&st2,&st1);  int i;

sprod(&sag,&sbg,&st2);

sadd(&st2,&st1);  reps
sequal(&st1,&seg);  rt1=rmaxabs(ps);
/∗ return result ∗/  scalar
for(i=0;i<=n;i++)  st1,st2,st3;

sequal(sx+i,pvres−>p+i);

sequal(&seg,&pvres−>g);  if(rt1>=rone)
sequal(&seh,&pvres−>h);  printf(”error in vsetin, radius= %30.16e\n”,rt1);
}  sequal(&sone,pv−>p);

for(i=2;i<=n1;i++){

vargcmp(pv,pvres)  sprod(pv−>p+i−2,ps,&st1);
vector ∗pv,∗pvres;  sneg(&st1,pv−>p+i−1);
{  }
extern reps  sequal(&szero,&pv−>g);
rmaxabs();  sabs(ps,&st1);

spower(&st1,n1,&st2);

int i;  sequal(&sone,&st3);

sminus(&st1,&st3);

reps rt1;  squot(&st2,&st3,&st1);

sscalar(rzero,st1.u,&pv−>h);

scalar  }

st1,st2;

vsetms(pv,ps)

vector  vector ∗pv;
vt1,vt2;  scalar ∗ps;

{

extern reps

/∗ calculate polymomial part ∗/  rmaxabs();

vsetzo(&vt1);

for(i=1;i<=n;i++){  int i;

vequal(&vpower[n1−i],&vt2);

vsmult(&vt2,pv−>p+n1−i,&vt2);  reps
vadd(&vt2,&vt1);  rt1=rmaxabs(ps);

}

sadd(pv−>p,vt1.p);  scalar

st1,st2,si,sf;

/∗ add error terms ∗/

srconst(r0save,&st1);  if(rt1>=rone)
sabs(&st1,&st1);  printf(”error in vsetms, radius= %30.16e\n”,rt1);

spower(&st1,n1,&st1);

srconst(r0save,&st2);  sequal(&sone,pv−>p);
sabs(&st2,&st2);  sneg(&shalf,&sf);
sprod(&st2,&pv−>g,&st2);  sprod(&sf,ps,&sf);
sprod(&st1,&pv−>h,&st1);  for(i=2;i<=n1;i++){
sadd(&st2,&st1);  sequal(&sf,pv−>p+i−1);
rt1=rmaxabs(&st1);  siconst(i,&si);
sscalar(−rt1,rt1,&st1);  sequal(&shalf,&st1);
sadd(&st1,vt1.p);  sminus(&si,&st1);
srconst(r3save,&st1);  squot(&st1,&si,&st1);
srconst(r1save,&st2);  sprod(&sf,&st1,&sf);
sprod(&st1,&st2,&st1);  sprod(ps,&sf,&sf);
sprod(&st1,&pv−>h,&st1);  }
srconst(r2save,&st2);  sequal(&szero,&pv−>g);
sprod(&st2,&pv−>g,&st2);  sequal(&sone,&st1);
sadd(&st2,&st1);  sabs(ps,&st2);
rt1=rmaxabs(&st1);  sminus(&st2,&st1);

sabs(&sf,&st2);

squot(&st2,&st1,&st1);  sneg(ps,psres)
sscalar(rzero,st1.u,&pv−>h);  scalar ∗ps,∗psres;
}  {

reps

r=−ps−>u;

/∗ 5.7. ROUTINES ACTING ON SCALARS ∗/

/∗ ============================ ∗/  psres−>u=−ps−>l;

psres−>l=r;

/∗ lsequal : compare two scalars ∗/ }
/∗ rmaxabs : largest distance from zero  ∗/
/∗ rminabs : smallest distance from zero  ∗/  sinv(ps,psres)
/∗ siconst : convert integer to scalar  ∗/  scalar ∗ps,∗psres;
/∗ srconst : convert rep to scalar  ∗/  {
/∗ sscalar : interval with given boundary  ∗/  extern reps
/∗ sabs  : absolute value of scalar  ∗/  rup();
/∗ sneg  : zero minus scalar  ∗/
/∗ sinv  : one divided by scalar  ∗/  if(ps−>l <= rzero && rzero <= ps−>u){
/∗ sequal : copies scalar  ∗/  printf(”error in sinv %30.16e %30.16e\n”,ps−>l,ps−>u);
/∗ sinter : intersection of two scalars  ∗/  sequal(ps,psres);
/∗ sadd  : sum of two scalars  ∗/  }
/∗ sminus : difference of two scalars  ∗/  else
/∗ sprod  : product of two scalars  ∗/  sscalar(rdquot(rone,ps−>u),
/∗ squot  : quotient of two scalars  ∗/  ruquot(rone,ps−>l),psres);
/∗ spower : scalar to the power k  ∗/  }
/∗ ssqrt  : square root of scalar  ∗/

sequal(ps,psres) scalar ∗ps,∗psres;

logical lsequal(ps1,ps2)  {
scalar ∗ps1,∗ps2;  psres−>l=ps−>l;
{  psres−>u=ps−>u;
return((ps1−>l==ps2−>l &&  }
ps1−>u==ps2−>u)? 1:0);
}  sinter(ps1,ps2,psres)

scalar ∗ps1,∗ps2,∗psres;

reps rmaxabs(ps)  {
scalar ∗ps;  reps
{  rl=rmax2(ps1−>l,ps2−>l),
return(rmax2(−ps−>l,ps−>u));  ru=rmin2(ps1−>u,ps2−>u);

}

if(rl>ru){

reps rminabs(ps)  printf(”error: no intersection. s1,s2= \n”);
scalar ∗ps;  printf(”%30.16e %30.16e %30.16e %30.16e \n”
{  ,ps1−>l,ps1−>u,ps2−>l,ps2−>u);
return(rmax3(rzero,ps−>l,−ps−>u));  sequal(ps1,psres);
}  }

else

siconst(i,ps)  sscalar(rl,ru,psres);
int i;  }

scalar ∗ps;

{  sadd(ps1,psres)
ps−>l=(reps)i;  scalar ∗ps1,∗psres;
ps−>u=(reps)i;  {
}  extern reps

rup();

srconst(r,ps)

reps r;  reps
scalar ∗ps;  r1l=ps1−>l,
{  r1u=ps1−>u,
ps−>l=r;  r2l=psres−>l,
ps−>u=r;  r2u=psres−>u;

}

if(r1l==rzero){

sscalar(rl,ru,ps)  if(r1u==rzero)
reps rl,ru;  ;
scalar ∗ps;  else if(r2u==rzero)
{  sscalar(r2l,r1u,psres);
ps−>l=rl;  else
ps−>u=ru;  sscalar(r2l,rusum(r1u,r2u),psres);
}  }

else if(r2l==rzero){

sabs(ps,psres)  if(r2u==rzero)
scalar ∗ps,∗psres;  sequal(ps1,psres);
{  else if(r1u==rzero)
reps  sscalar(r1l,r2u,psres);
rl=ps−>l,  else
ru=ps−>u;  sscalar(r1l,rusum(r1u,r2u),psres);

}

psres−>l=rmax3(rzero,rl,−ru);  else{
psres−>u=rmax2(−rl,ru);  if(r1u==rzero)
}  sscalar(rdsum(r1l,r2l),r2u,psres);

else if(r2u==rzero)  m,mh;

sscalar(rdsum(r1l,r2l),r1u,psres);

else  scalar
sscalar(rdsum(r1l,r2l),  s1,st1;

rusum(r1u,r2u),psres);

}  if(k<0)
}  printf(”error in spower, k not positive\n”);

else if(rmaxabs(ps)==rzero){

sminus(ps1,psres)  if(k==0)
scalar ∗ps1,∗psres;  sequal(&sone,psres);
{  else
scalar s;  sequal(&szero,psres);

}

sneg(ps1,&s);  else{
sadd(&s,psres);  m=k;
}  sequal(ps,&s1);

sequal(&sone,&st1);

sprod(ps1,ps2,psres)  while(m>0){
scalar ∗ps1,∗ps2,∗psres;  mh=m/2;
{  if(mh+mh!=m){
extern reps  sprod(&st1,&s1,&st1);
rup();  −−m;

}

reps  else{
rlu,rul,ruu,rll,  sprod(&s1,&s1,&s1);
r1l=ps1−>l,  m=mh;
r1u=ps1−>u,  }
r2l=ps2−>l,  }
r2u=ps2−>u;  sequal(&st1,psres);

}

if(r1u==rzero && r1l==r1u||  }

r2u==rzero && r2l==r2u)

sequal(&szero,psres);  ssqrt(ps,psres)
else if(r1l>=rzero){  scalar ∗ps,∗psres;
if(r2l>=rzero)  {
sscalar(rdcprod(r1l,r2l),rucprod(r1u,r2u),psres);  extern double
else if(r2u<rzero)  sqrt();

sscalar(rdcprod(r1u,r2l),rucprod(r1l,r2u),psres);

else  extern reps
sscalar(rdcprod(r1u,r2l),rucprod(r1u,r2u),psres);  rup();

}

else if(r1u<rzero){  reps
if(r2l>=rzero)  rl,ru;

sscalar(rdcprod(r1l,r2u),rucprod(r1u,r2l),psres);

else if(r2u<rzero)  if(ps−>l>=rzero){
sscalar(rdcprod(r1u,r2u),rucprod(r1l,r2l),psres);  rl=rmax2(rzero,sqrt(ps−>l));
else  while(rucprod(rl,rl)>ps−>l)
sscalar(rdcprod(r1l,r2u),rucprod(r1l,r2l),psres);  rl=rdown(rl);
}  ru=rmax2(rzero,sqrt(ps−>u));
else{  while(rdcprod(ru,ru)<ps−>u)
if(r2l==rzero)  ru=rup(ru);
sscalar(rdcprod(r1l,r2u),rucprod(r1u,r2u),psres);  sscalar(rl,ru,psres);
else if(r2u<rzero)  }
sscalar(rdcprod(r1u,r2l),rucprod(r1l,r2l),psres);  else{
else{  printf(”error in sqrt, argument negative %30.16d\n”,ps−>l);
rlu=rdcprod(r1l,r2u);  sequal(&szero,psres);
rul=rdcprod(r1u,r2l);  }
ruu=rucprod(r1u,r2u);  }

rll=rucprod(r1l,r2l); sscalar(rmin2(rlu,rul),rmax2(ruu,rll),psres);

}  /∗ 5.8. OPERATIONS INVOLVING ONLY REPS ∗/
}  /∗ =============================== ∗/

}

/∗ rup  : successor in the set of reps  ∗/
squot(ps1,ps2,psres)  /∗ rdown : predecessor in the set of reps  ∗/
scalar ∗ps1,∗ps2,∗psres;  /∗ rmax2 : maximum of two reps  ∗/
{  /∗ rmax3 : maximum of three reps  ∗/
scalar s;  /∗ rmin2 : minimum of two reps  ∗/
/∗ rmin3 : minimum of three reps  ∗/
sinv(ps2,&s); /∗ rusum : upwards rounded sum of reps ∗/
sprod(ps1,&s,psres);  /∗ rdsum : downwards rounded sum of reps ∗/
}  /∗ rddiff : downwards rounded diff. of reps ∗/
/∗ rucprod : upper bound on product of reps  ∗/
spower(ps,k,psres) /∗ rdcprod : lower bound on product of reps ∗/
scalar ∗ps,∗psres;  /∗ ruquot : upper bound on quotient of reps ∗/
int k; /∗ rdquot : lower bound on quotient of reps ∗/

{

extern reps  /∗ Remark.
rmaxabs();  The type reps (=double) is in 64−bit IEEE format.

The standard IEEE format ( see [24] ) is as follows.

int  Sign bit, 11−bit binary exponent with bias

of 1023, 52−bit mantissa ( with ”hidden” most

significant bit=1 ). The mantissa represents a  #define rucprod(x,y) (first(x,y):(second(x,y):rup((x)∗(y))))

number < 2 but not < 1. A zero exponent means a

value of zero, and the exponent 11111111111 means  #define rdcprod(x,y) (first(x,y):(second(x,y):rdown((x)∗(y))))

”not a number”. Bytes are in reversed order;

the lowest addressed byte is the least  #define ruquot(x,y) (first(x,y):rup((x)/(y)))

significant mantissa byte. ∗/

#define rdquot(x,y) (first(x,y):rdown((x)/(y)))

reps rup(r) reps r;

{  ! 5.9 COEFFICIENTS OF THE
union convert{ ! APPROXIMATE FIXED POINT
reps r;  ! =====================

unsigned long int i[2];

};  0.0000000000000000e+000

0.0000000000000000e+000

union convert  5.0910744110871690e−001
con;  9.5003399093995690e−001

3.3399471598550260e−001

unsigned long int  7.5228788724207240e−002
low,exp and high;  2.9166025360328340e−002

9.7809075474681730e−003

con.r=r;  2.8435242149942920e−003
low=con.i[0];  7.1603661837487370e−004
exp and high=con.i[1];  1.5546221099374870e−004

2.8821170499654220e−005

if(low==0){  4.4811455122526670e−006
if(exp and high==0)  5.6476122738011440e−007
exp and high=0X00100000;  5.3552629703583600e−008
else if(r>rzero)  3.0116983547275270e−009
low=1;  −5.6475304807471550e−011
else if(exp and high==0X80100000)  −3.0876243410073330e−011
exp and high=0;  −2.5344800033482360e−012
else{  4.5983279951728720e−014
low=0XFFFFFFFF;  2.5517826728299680e−014
−−exp and high;  1.2620701467321130e−015
}  −1.5954295084724120e−016
}  −1.9577090887265150e−017
else if(r>rzero){  5.6384521325252770e−019
if(low==0XFFFFFFFF){  1.6946227528291020e−019
if(exp and high==0X7FFFFFF){  −5.8050871684659000e−021
printf(”overflow error in up %30.16e\n”,r);  −1.5252782196365650e−021
return(r);  3.5918686798167610e−022
}  8.9115968244695160e−023
else{  −4.8724320918072530e−024
low=0;  −5.1256140504614330e−024
++exp and high;  −5.5171555999740080e−025
}  2.0355234891894370e−025
}  5.5448722017939900e−026
else  −6.2413537679217950e−027
++low;  −3.6414430696016150e−027
}  1.6372588275807340e−028
else  2.1526142159118610e−028
−−low;  −6.1726379503092100e−030

−1.2428705431372500e−029

con.i[0]=low;  5.9617241343875310e−031
con.i[1]=exp and high;  6.9991573813170870e−031
return(con.r);  −6.9906713111864140e−032
}  −3.6182464583042550e−032

6.9244538361489700e−033

/∗ Reps macros ∗/  1.4733779930696580e−033

−5.5893664296982620e−034

#define rdown(r) (−rup(−(r)))  −2.0798556973059780e−035

3.5409929658431840e−035

#define rmax2(x,y) (((x)>(y))? (x):(y))  −3.7338449518599320e−036

−1.5159188525456280e−036

#define rmax3(x,y,z) (((x)>(y))? rmax2(x,z):rmax2(y,z))  4.6350294425125140e−037

9.9086553208077940e−039

#define rmin2(x,y) (((x)<(y))? (x):(y))  −2.9079748380155910e−038

4.9369413872234490e−039

#define rmin3(x,y,z) (((x)<(y))? rmin2(x,z):rmin2(y,z)) #define rusum(x,y) rup((x)+(y)) #define rdsum(x,y) rdown((x)+(y)) #define rddiff(x,y) (−rup((y)−(x))) #define first(x,y) ((x)==rzero||(y)==rone)?(x) #define second(x,y) ((x)==rone||(y)==rzero)?(y)