
      program main
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      parameter(imat=9)
      common/rconst/rzero,rone,rtwo
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      common/output/iout
      common/model /scov(0:1),sarg(0:1)
      common/norm  /sigma(0:1,1:imat,1:imat),sigmai(0:1,1:imat,1:imat),
     *		    sigmax(0:1),sigmxi(0:1)
      dimension rr(1:imat,1:imat)
      dimension p0(0:1,0:nt3,0:nu3),p0i(0:1,0:nt3,0:nu3)
      dimension pd0(0:1,0:nt3,0:nu3),pd0i(0:1,0:nt3,0:nu3)
      dimension p0save(0:1,0:nt3,0:nu3)
      dimension pt1(0:1,0:nt3,0:nu3),pt2(0:1,0:nt3,0:nu3)
      dimension sdelta(0:1),snu(0:1)
      dimension s(0:1),s1(0:1),s2(0:1),si(0:1),si1(0:1)
c
c  Initializations
c  ===============
c..open output file (this file will contain the results)
      open(12,file='RESULTS')
      iout=12
c..initialize constants
      rzero=0.D+00
      rone =1.D+00
      rtwo =2.D+00
      call siconst(0,szero)
      call siconst(1,sone)
      call siconst(2,stwo)
      call siconst(8,seight)
      call sinv2(stwo,shalf)
c..the following constants are used to bound derivatives
      call sexp(shalf,s)
      call ssqrt(s,s1)
      call sinv2(s1,semp25)
      call smult2(stwo,s1)
      call smult2(stwo,s1)
      call sexp(sone,s)
      call squot(s1,s,s4over)
c..read in (9x9)-matrix with change of basis
      open(11,file='SS')
      do 50 j=1,imat
      do 51 i=1,imat
      read(11,*)rr(i,j)
   51 continue
   50 continue
      close(11)
c..compute inverse matrix
      call matas(rr,sigma)
      call mati(sigma,sigmai)
c..compute constants for change of norm
      call siconst(0,sigmax)
      call siconst(0,sigmxi)
      do 60 j=1,imat
      call siconst(0,s)
      call siconst(0,si)
      do 70 i=1,imat
      call sabs2(sigma(0,i,j),s1)
      call sadd(s1,s)
      call sabs2(sigmai(0,i,j),si1)
      call sadd(si1,si)
   70 continue
      call sscalar(rzero,rmax2(s(1),sigmax(1)),sigmax)
      call sscalar(rzero,rmax2(si(1),sigmxi(1)),sigmxi)
   60 continue
      call sscalar(rzero,rmax2(rone,sigmax(1)),sigmax)
      call sscalar(rzero,rmax2(rone,sigmxi(1)),sigmxi)
c..define analyticity domains
      rtdom=3.D+00
      rudom=1.D-02
c..define model parameters (arg=(alpha/L^2)^2, 1/(2cov)=L^2)
      call siconst(32,s)
      call sinv2(s,sarg)
      call sinv2(seight,scov)
c..read input
      call pread(p0,rtdom,rudom)
      goto 111
c
c
c  Compute Precision of Approximate Fixed Point
c  ============================================
      call fr(p0,p0i,sdelta,1)
      call pequal(p0i,pt1)
      call psubtr(p0,pt1)
      reps=r1norm(pt1)
      write(iout,*)'epsilon=',reps,sdelta
c      call pack(pt1,pt2)
c      do 10 iu=0,nu3
c      do 11 it=0,nt3
c      write(12,*)pt2(0,it,iu),pt2(1,it,iu)
c   11 continue
c   10 continue
c..since the map is a contraction p0i is a better approximation
      call pwrite(p0i)
      close(12)
      stop
c
c
c  Estimate Norm of Tangent Map
c  ============================
c..save the original input
  111 continue
      call pequal(p0,p0save)
c..prepare p0 on neighborhood with radius beta
      rbeta=5.D-9
      call pack(p0,pt1)
      call sscalar(rzero,rbeta,pt1(0,nt+1,nu+1))
      call unpack(pt1,p0,rtdom,rudom,1)
c..estimate norm of tangent map
c..initialize norm of upper bound to zero
      rtheta=rzero
c..fix number of computed basis vectors
      iumax=2
      itmax=9
c..loop over tangent vectors
      goto 999
c      do 2 iu=0,iumax
c      do 3 it=1,itmax
      do 2 iu=2,2
      do 3 it=1,1
      if((iu.eq.1).and.(it.eq.1))goto 3
      call basvec(it,iu,pd0,rtdom,rudom)
      call fdr(p0,pd0,pd0i)
      write(iout,*)it,iu,r1norm(pd0i)
      rtheta=rmax2(rtheta,r1norm(pd0i))
    3 continue
    2 continue
999   continue
c..estimate norm of higher order
c..the following is the matrix which needs to be initialized
c------------|----------|---|---|
c  already   |	   Z	| Z | Z |
c  done      |	   E	| E | E |
c	     |	   R	| R | R |
c	     |	   O	| O | O |
c------------|----------|---|---|
c 1   ...   1|	   Z	| Z | Z |
c .	    .|		|   |	|
c .	    .|	   E	| E | E |
c	     |		|   |	|
c .	    .|	   R	| R | R |
c .	    .|		|   |	|
c 1   ...   1|	   O	| O | O |
c------------|----------|---|---|
c   ZERO     |1  ....  1| 0 | 1 |
c------------|----------|---|---|
c 1   ...   1|	 ZERO	| 0 | 0 |
c------------|----------|---|---|
c..initialize l1-sequence of higher order
      call pstzro(pt1)
c..initialize higher powers in t for low powers of u
      do 4 iu=0,iumax
      do 5 it=itmax+1,nt
      call sscalar(rzero,rone,pt1(0,it,iu))
    5 continue
      call sscalar(rzero,rone,pt1(0,nt+2,iu))
    4 continue
c..initialize  generals in t for higher powers of u
      do 6 iu=iumax+1,nu
      call sscalar(rzero,rone,pt1(0,nt+1,iu))
    6 continue
      call sscalar(rzero,rone,pt1(0,nt+1,nu+2))
c..map to space of analytic functions
      call unpack(pt1,pd0,rtdom,rudom,0)
c..apply tangent map
      call fdr(p0,pd0,pd0i)
      write(iout,*)'higher order ',r1norm(pd0i)
      rtheta=rmax2(rtheta,r1norm(pd0i))
c
c  Check if Contraction
c  ====================
      if(reps.lt.rdprod(rddiff(rone,rtheta),rbeta))then
        write(iout,*)'contraction on ball, fixed point exists.'
c
c  Compute Bounds on Critical Indices
c  ==================================
	call pack(p0save,pt1)
	call sscalar(rzero,ruquot(reps,rddiff(rone,rtheta)),
     *		     pt1(0,nt+1,nu+1))
	call unpack(pt1,p0,rtdom,rudom,1)
	call fr(p0,p0i,sdelta,1)
        write(iout,*)sdelta(0),'<= delta <=',sdelta(1)
	call slog(stwo,s1)
	call slog(sdelta,s2)
	call squot(s1,s2,snu)
        write(iout,*)snu(0),'<= nu <=',snu(1)
      else
        write(iout,*)'fixed point does not exist, estimation of ',
     *               'critical index not done.'
      endif
c
c..close output file
      close(12)
      end
c-----------------------------------------------------------------------
c     INPUT/OUTPUT OPERATIONS
c     =======================
c
c     twrite : write the coefficients of a function of t
c     pwrite : write the coefficients of a function of t and u
c     tread  : read the coefficients of a function of t
c     pread  : read the coefficients of a function of t and u
c-----------------------------------------------------------------------
      subroutine twrite(t)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
   99 format(D26.16)
      do 1 i=0,nt
      if(t(0,i).eq.t(1,i))then
	write(11,99)t(0,i)
      else
	write(11,99)(t(0,i)+t(1,i))/2
      endif
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine pwrite(p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3)
  100 format(I26)
      open(11,file='INPUT')
      write(11,100)nt
      write(11,100)nu
      do 1 i=0,nu
      call twrite(p(0,0,i))
    1 continue
      close(11)
      return
      end
c-----------------------------------------------------------------------
      subroutine tread(t,k)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
   99 format(D26.16)
      do 1 i=0,k
      read(11,99)r
      if(i.gt.nt)goto 1
      call srconst(r,t(0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine pread(p,rtdom,rudom)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3)
  100 format(I26)
c..define domains
      do 1 i=0,nu+2
      p(0,nt3,i)=rtdom
    1 continue
      p(0,0,nu3)=rudom
c..set coefficients zero
      call pstzro(p)
c..read coefficients
      open(11,file='INPUT')
      read(11,100)kt
      read(11,100)ku
      do 3 i=0,ku
      if(i.gt.nu)goto 3
      call tread(p(0,0,i),kt)
    3 continue
      close(11)
      return
      end
c-----------------------------------------------------------------------
c     DEFINITION OF NORMS AND BASIS  VECTORS
c     ======================================
c
c     tpack  : map: analytic functions -> l-1 spaces
c     tunpck : map: l-1 spaces -> analytic functions
c     pack   : map: analytic functions -> direct-sum of l-1 spaces
c     unpack : map: direct-sum of l-1 spaces -> analytic functions
c     basvec : make basisvector in space of analytic functions
c     st1nrm : compute the norm of an elment in direct-sum of l-1-spaces
c     r1norm : compute the norm of an element in l-1
c-----------------------------------------------------------------------
      subroutine tpack(t)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
      dimension s(0:1),s1(0:1),srhot(0:1)
c..maps the analytic function t into a l-1-sequence
c..extract domain
      call srconst(t(0,nt3),srhot)
c..pack coefficients and error terms
      do 1 it=1,nt+2
      call sfak(itgh(it),s)
      call spower(srhot,itgh(it),s1)
      call smult2(s1,s)
      call smult2(s,t(0,it))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine tunpck(t,rtdom)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
      dimension s(0:1),s1(0:1),srhot(0:1)
c..maps a l-1-sequence into a function,
c  analytic on a domain of radius rtdom.
c..make scalar from domain
      call srconst(rtdom,srhot)
c..unpack coefficients
      do 1 it=1,nt+2
      call sfak(itgh(it),s)
      call spower(srhot,itgh(it),s1)
      call smult2(s1,s)
      call sinv1(s)
      call smult2(s,t(0,it))
    1 continue
c..store domain
      t(0,nt3)=rtdom
      return
      end
c-----------------------------------------------------------------------
      subroutine pack(p,pi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      parameter(imat=9)
      common/norm  /sigma(0:1,1:imat,1:imat),sigmai(0:1,1:imat,1:imat),
     *		    sigmax(0:1),sigmxi(0:1)
      dimension p(0:1,0:nt3,0:nu3),pi(0:1,0:nt3,0:nu3)
      dimension v(0:1,1:imat)
      dimension s(0:1),s1(0:1),srhou(0:1)
c..extract domain
      call srconst(p(0,0,nu3),srhou)
c..transfer coefficients
      call pequal(p,pi)
c..take care of normalizations
      call siconst(0,pi)
      call siconst(0,pi(0,1,1))
c..pack functions
      do 1 iu=0,nu+2
      call tpack(pi(0,0,iu))
      call spower(srhou,iugh(iu),s)
      call tsmlt2(s,pi(0,0,iu))
    1 continue
c..apply coordinate transform on the first imat coefficients of F(t,0)
c..apply transformation to coefficients and store result in v
      do 3 i=1,imat
      call siconst(0,s)
      do 4 j=1,imat
      call smult3(sigmai(0,i,j),pi(0,j,0),s1)
      call sadd(s1,s)
    4 continue
      call sequal(s,v(0,i))
    3 continue
c..transfer result form v to pi
      do 5 i=1,imat
      call sequal(v(0,i),pi(0,i,0))
    5 continue
c..resize general error term for mu=0
      call smult2(sigmxi,pi(0,nt+1,0))
      return
      end
c-----------------------------------------------------------------------
      subroutine unpack(p,pi,rtdom,rudom,ispace)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      parameter(imat=9)
      common/norm  /sigma(0:1,1:imat,1:imat),sigmai(0:1,1:imat,1:imat),
     *		    sigmax(0:1),sigmxi(0:1)
      dimension p(0:1,0:nt3,0:nu3),pi(0:1,0:nt3,0:nu3)
      dimension v(0:1,1:imat)
      dimension s(0:1),s1(0:1),srhou(0:1)
c..make scalar with domain
      call srconst(rudom,srhou)
c..transfer coefficients
      call pequal(p,pi)
c..apply coordinate transform to first imat coefficients of F(t,0)
c..apply transformation to coefficients and store result in v
      do 3 i=1,imat
      call siconst(0,s)
      do 4 j=1,imat
      call smult3(sigma(0,i,j),pi(0,j,0),s1)
      call sadd(s1,s)
    4 continue
      call sequal(s,v(0,i))
    3 continue
c..transfer result form v to pi
      do 5 i=1,imat
      call sequal(v(0,i),pi(0,i,0))
    5 continue
c..resize general error term for mu=0
      call smult2(sigmax,pi(0,nt+1,0))
c..unpack functions
      do 1 iu=0,nu+2
      call tunpck(pi(0,0,iu),rtdom)
      call spower(srhou,iugh(iu),s)
      call sinv1(s)
      call tsmlt2(s,pi(0,0,iu))
    1 continue
c..domain
      pi(0,0,nu3)=rudom
c..fix normalizations for tangent space or basepoint
      call siconst(ispace,pi)
      call siconst(ispace,pi(0,1,1))
      return
      end
c-----------------------------------------------------------------------
      subroutine basvec(it,iu,pd0,rtdom,rudom)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension pd0(0:1,0:nt3,0:nu3)
      dimension pt1(0:1,0:nt3,0:nu3)
c..define coefficients
      call pstzro(pt1)
      call siconst(1,pt1(0,it,iu))
c..map to space of analytic functions
      call unpack(pt1,pd0,rtdom,rudom,0)
      return
      end
c-----------------------------------------------------------------------
      subroutine st1nrm(t,sres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),sres(0:1)
      dimension s(0:1),s1(0:1)
      call siconst(0,s)
      do 1 i=0,nt+2
      call sbound2(t(0,i),s1)
      call sadd(s1,s)
    1 continue
      call sbound2(s,sres)
      return
      end
c-----------------------------------------------------------------------
      function r1norm(p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3)
      dimension p0(0:1,0:nt3,0:nu3)
      dimension s(0:1),s1(0:1)
      call pack(p,p0)
      call siconst(0,s)
      do 1 i=0,nu+2
      call st1nrm(p0(0,0,i),s1)
      call sadd(s1,s)
    1 continue
      r1norm=s(1)
      return
      end
c-----------------------------------------------------------------------
c     MAPS IN FUNCTION SPACE
c     ======================
c
c     fr     : the extended renormalization group map
c     fdr    : the linearization of fr
c     fnstar : the map fn on the massless manifold
c     fdnstar: the linearization of fnstar
c     fs2    : shift in mu-argument (two arguments in subroutine)
c     fs3    : shift in mu-argument (three arguments in subroutine)
c     sustar : find critical mass
c     sdusta : the linearization of sustar
c     fstar  : localize the massless manifold
c     fdstar : the linearization of fstar
c     fn     : the standard RG operator with mass zero
c     fdn    : the linearization of fn
c     snomas : compute difference between a given mass and its image
c     sdnoma : the linearization of snomas
c-----------------------------------------------------------------------
      subroutine fr(p,pi,sdelta,irestr)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/output/iout
      dimension p(0:1,0:nt3,0:nu3),pi(0:1,0:nt3,0:nu3),sdelta(0:1)
      dimension s(0:1)
      print*,'in fr'
      call pcheck(p)
c..applies the extended renormalization group map R(.)
      call fnstar(p,pi)
      call scoefe(pi,1,1,sdelta)
      call sinv2(sdelta,s)
      call p2sdil(pi,s)
      call siconst(1,pi(0,1,1))
c..check that domain of image is larger
      if(pi(0,0,nu3).lt.p(0,0,nu3))
     *write(iout,*)'domain error in u, rhou=',pi(0,0,nu3)
      if(pi(0,nt3,0).lt.p(0,nt3,0))
     *write(iout,*)'domain error in t, rhot=',pi(0,nt3,0)
c..return if restriction of domain not requested
      if(irestr.eq.0)return
c..restrict domain
      do 1 i=0,nu+2
      pi(0,nt3,i)=p(0,nt3,0)
    1 continue
      pi(0,0,nu3)=p(0,0,nu3)
      print*,'in fr'
      call pcheck(pi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fdr(p,dp,dpi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      common/output/iout
      dimension p(0:1,0:nt3,0:nu3),
     *		dp(0:1,0:nt3,0:nu3),dpi(0:1,0:nt3,0:nu3)
      dimension pt1(0:1,0:nt3,0:nu3),pt2(0:1,0:nt3,0:nu3)
      dimension u(0:1,0:nu3)
      dimension sdelta(0:1),sddelt(0:1),s(0:1),s1(0:1),s2(0:1)
      print*,'in fdr'
      call pcheck(p)
      print*,'in fdr'
      call pcheck(dp)
c..applies the linearized map DR(.)
c..compute second term
      call fr(p,pt1,sdelta,0)
      call p2der(pt1,pt2)
      call sneg2(sone,s)
      call squot(s,sdelta,s1)
      call ustlin(u,szero,s1,pt2(0,0,nu3))
      call p2mult(pt2,u)
c..compute first term
      call fdnstar(p,dp,dpi)
      call sinv2(sdelta,s)
      call p2sdil(dpi,s)
c..compute ddelta
      call scoefe(dpi,1,1,s1)
      call scoefe(pt2,1,1,s2)
      call squot(s1,s2,s)
      call sneg2(s,sddelt)
      call psmadd(pt2,sddelt,dpi)
      call siconst(0,dpi(0,1,1))
c..check that domain of image is larger
      if(dpi(0,0,nu3).lt.dp(0,0,nu3))
     *write(iout,*)'domain error in u, rhou=',dpi(0,0,nu3)
      if(dpi(0,nt3,0).lt.dp(0,nt3,0))
     *write(iout,*)'domain error in t, rhot=',dpi(0,nt3,0)
c..restrict domain
      do 1 i=0,nu+2
      dpi(0,nt3,i)=dp(0,nt3,0)
    1 continue
      dpi(0,0,nu3)=dp(0,0,nu3)
      print*,'in fdr'
      call pcheck(dpi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fnstar(p,pi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),pi(0:1,0:nt3,0:nu3)
      dimension pt1(0:1,0:nt3,0:nu3)
c..applies the map N(S*(.))
      print*,'in fnstar'
      call pcheck(p)
      call fstar(p,pt1)
      print*,'in fnstar'
      call pcheck(pt1)
      call fn(pt1,pi)
      print*,'in fnstar'
      call pcheck(pi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fdnstar(p,dp,dpi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),
     *		dp(0:1,0:nt3,0:nu3),dpi(0:1,0:nt3,0:nu3)
      dimension p1(0:1,0:nt3,0:nu3),dp1(0:1,0:nt3,0:nu3)
c..computes DN DS*(.)
      print*,'in fdnstar'
      call pcheck(p)
      print*,'in fdnstar'
      call pcheck(dp)
      call fdstar(p,dp,dp1)
      print*,'in fdnstar'
      call pcheck(dp1)
      call fstar(p,p1)
      print*,'in fdnstar'
      call pcheck(p1)
      call fdn(p1,dp1,dpi)
      print*,'in fdnstar'
      call pcheck(dpi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fs2(p,su0)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension p(0:1,0:nt3,0:nu3),su0(0:1)
      dimension u0(0:1,0:nu3)
c..computes S(.)
      print*,'in fs'
      call pcheck(p)
c..choose domain for composition as large as possible.
      rrho=rddiff(p(0,0,nu3),rmaxabs(su0))
      call ustlin(u0,su0,sone,rrho)
      call p2comp(p,u0)
      print*,'in fs'
      call pcheck(p)
      return
      end
c-----------------------------------------------------------------------
      subroutine fs3(p,su0,pi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),su0(0:1),pi(0:1,0:nt3,0:nu3)
c..computes S(.)
      call pequal(p,pi)
      call fs2(pi,su0)
      return
      end
c-----------------------------------------------------------------------
      subroutine sustar(p,star)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/rconst/rzero,rone,rtwo
      common/output/iout
      dimension p(0:1,0:nt3,0:nu3),star(0:1)
      dimension su0(0:1),snu0(0:1)
      dimension seps(0:1)
      dimension st1(0:1),st2(0:1)
      logical lsequal
c..compute mu0* by Newtons method
c..first guess of center
      rcent=0.D+00
c..find a good guess for center (not done for final run)
      rder=1.904068D+00
      do 1 j=1,1
      call srconst(rcent,st1)
      call snomas(p,st1,st2)
      rcent=rcent-(st2(0)+st2(1))/(2*rder)
c      write(iout,*)'sustar:rc,prc',rcent,st2
    1 continue
c..calculate precision of center
      call srconst(rcent,st1)
      call snomas(p,st1,seps)
c..enlarge interval for seps to contain zero
      call sscalar(rmin2(seps(0),rzero),rmax2(rzero,seps(1)),seps)
c..estimate derivative on neighborhood
      ru0=1.D-05
      call sscalar(-ru0,ru0,su0)
      do 2 j=1,2
      call srconst(rcent,st1)
      call sadd(su0,st1)
      call sdnoma(p,st1,st2)
      call sneg2(st2,st1)
      call squot(seps,st1,snu0)
      call sinter(su0,snu0,st1)
      if(.not.lsequal(snu0,st1))then
        write(iout,*)'error in sustar, no contraction. step',j,
     *               'snu0=',snu0
	goto 3
      endif
      call sequal(snu0,su0)
    2 continue
    3 continue
      call srconst(rcent,star)
      call sadd(su0,star)
      return
      end
c-----------------------------------------------------------------------
      subroutine sdusta(p,dp,sdstar)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension p(0:1,0:nt3,0:nu3),dp(0:1,0:nt3,0:nu3),sdstar(0:1)
      dimension pt1(0:1,0:nt3,0:nu3),pt2(0:1,0:nt3,0:nu3),
     *		pt3(0:1,0:nt3,0:nu3)
      dimension su0(0:1)
      dimension st1(0:1),st2(0:1),st3(0:1)
c..compute variation of mu0*
      print*,'in sdusta'
      call pcheck(p)
      print*,'in sdusta'
      call pcheck(dp)
      call sustar(p,su0)
      call fs3(p,su0,pt1)
      call fs3(dp,su0,pt2)
      call fdn(pt1,pt2,pt3)
      call scoefe(pt3,1,0,st1)
      call scoefa(pt3,2,0,st1)
      call scoefm(pt2,1,0,st1)
      call scoefm(pt2,2,0,st1)
      call smult2(shalf,st1)
      call sdnoma(p,su0,st2)
      call squot(st1,st2,st3)
      call sneg2(st3,sdstar)
      return
      end
c-----------------------------------------------------------------------
      subroutine fstar(p,pi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),pi(0:1,0:nt3,0:nu3)
      dimension su0(0:1)
c..compute S*(.)
      print*,'in fstar'
      call pcheck(p)
      call sustar(p,su0)
      call fs3(p,su0,pi)
      print*,'in fstar'
      call pcheck(pi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fdstar(p,dp,dpi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),
     *		dp(0:1,0:nt3,0:nu3),dpi(0:1,0:nt3,0:nu3)
      dimension su0(0:1),sdu0(0:1)
      dimension pt1(0:1,0:nt3,0:nu3)
c..compute DS*(.)
      print*,'in fdstar'
      call pcheck(p)
      print*,'in fdstar'
      call pcheck(dp)
      call sdusta(p,dp,sdu0)
      call sustar(p,su0)
      call fs3(dp,su0,dpi)
      print*,'in fdstar'
      call pcheck(pi)
      call p2der(p,pt1)
      print*,'in fdstar'
      call pcheck(pt1)
      call fs2(pt1,su0)
      print*,'in fdstar'
      call pcheck(pt1)
      call psmadd(pt1,sdu0,dpi)
      print*,'in fdstar'
      call pcheck(dpi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fn(p,pi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/model /scov(0:1),sarg(0:1)
      dimension p(0:1,0:nt3,0:nu3),pi(0:1,0:nt3,0:nu3)
      dimension pt2(0:1,0:nt3,0:nu3),pt4(0:1,0:nt3,0:nu3),
     *		pt8(0:1,0:nt3,0:nu3)
      dimension uc(0:1,0:nu3)
c..compute N(.)
      print*,'in fn'
      call pcheck(p)
      call pmult3(p,p,pt2)
      call pmult3(pt2,pt2,pt4)
      call pmult3(pt4,pt4,pt8)
      print*,'in fn'
      call pcheck(pt8)
      call p1int(pt8,scov)
      print*,'in fn'
      call pcheck(pt8)
      call p2eql(pt8,uc)
      call utkinv(uc)
      call p2mult(pt8,uc)
      print*,'in fn'
      call pcheck(pt8)
      call p2zro(pt8)
      call siconst(1,pt8)
      call p1sdil(pt8,sarg,pi)
      print*,'in fn'
      call pcheck(pi)
      return
      end
c-----------------------------------------------------------------------
      subroutine fdn(p,dp,dpi)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      common/model /scov(0:1),sarg(0:1)
      dimension p(0:1,0:nt3,0:nu3),
     *		dp(0:1,0:nt3,0:nu3),dpi(0:1,0:nt3,0:nu3)
      dimension pt2(0:1,0:nt3,0:nu3),pt4(0:1,0:nt3,0:nu3),
     *		pt6(0:1,0:nt3,0:nu3),pt7(0:1,0:nt3,0:nu3),
     *		pt8(0:1,0:nt3,0:nu3)
      dimension uc(0:1,0:nu3),udc(0:1,0:nu3)
      dimension u(0:1,0:nu3)
c..compute DN(.)
      print*,'in fdn'
      call pcheck(p)
      print*,'in fdn'
      call pcheck(dp)
c..prepare powers of p
      call pmult3(p,p,pt2)
      call pmult3(pt2,pt2,pt4)
      call pmult3(pt4,pt2,pt6)
      call pmult3(p,pt6,pt7)
      call pmult3(pt4,pt4,pt8)
c..pt7 is used to store the contribution of dp
      call psmult(seight,pt7)
      call pmult2(dp,pt7)
c..integrate
      call p1int(pt7,scov)
      call p1int(pt8,scov)
c..normalize f^8
      call p2eql(pt8,uc)
      call utkinv(uc)
      call p2mult(pt8,uc)
      call p2zro(pt8)
      call siconst(1,pt8)
c..compute udc
      call p2eql(pt8,udc)
      call utkinv(udc)
      call p2eql(pt7,u)
      call umult(u,udc)
c..add things up
      call p2mult(pt8,udc)
      call psubtr(pt8,pt7)
      call p2mult(pt7,uc)
      call p2zro(pt7)
      call p1sdil(pt7,sarg,dpi)
      print*,'in fdn'
      call pcheck(dpi)
      return
      end
c-----------------------------------------------------------------------
      subroutine snomas(p,su0,smass)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension p(0:1,0:nt3,0:nu3),su0(0:1),smass(0:1)
      dimension pt1(0:1,0:nt3,0:nu3),pt2(0:1,0:nt3,0:nu3)
c..computes the function m(.)
      call fs3(p,su0,pt1)
      call fn(pt1,pt2)
      call scoefe(pt2,1,0,smass)
      call scoefa(pt2,2,0,smass)
      call scoefm(pt1,1,0,smass)
      call scoefm(pt1,2,0,smass)
      call smult2(shalf,smass)
      return
      end
c-----------------------------------------------------------------------
      subroutine sdnoma(p,su0,sdmass)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension p(0:1,0:nt3,0:nu3),su0(0:1),sdmass(0:1)
      dimension pt1(0:1,0:nt3,0:nu3),pt2(0:1,0:nt3,0:nu3),
     *		pt3(0:1,0:nt3,0:nu3)
c..needed to compute dmu0*
      call fs3(p,su0,pt1)
      call p2der(p,pt2)
      call fs2(pt2,su0)
      call fdn(pt1,pt2,pt3)
      call scoefe(pt3,1,0,sdmass)
      call scoefa(pt3,2,0,sdmass)
      call scoefm(pt2,1,0,sdmass)
      call scoefm(pt2,2,0,sdmass)
      call smult2(shalf,sdmass)
      return
      end
c-----------------------------------------------------------------------
c     OPERATIONS ON FUNCTIONS OF t and u
c     ==================================
c
c     scoefe : extract bounds for a coefficient
c     spack  : scale coeficients and error terms with domains
c     sunpac : invers of spack
c     scoefa : add a coefficient
c     scoefm : subtract a coefficient
c     pstzro : set function equal to zero (does not change domain)
c     p1zro  : set the function f(t,0) equal to zero (modulo general term)
c     p2zro  : set the function f(0,u) equal to zero (modulo general term)
c     pequal : set two functions equal
c     p1eql  : copy the function f(t,0) to a function of t only
c     p2eql  : copy the function f(0,u) to a function of u only
c     padd   : add two functions
c     p2der  : take the derivative with respect to the second argument
c     psubtr : subtract two functions
c     psmult : multiply a function with a scalar
c     psmadd : multiply a function with a number and add it
c     p1sdil : dilate a function in its first argument
c     p2sdil : dilate the second argument of a function
c     pmult2 : multiply two functions (two subroutine arguments)
c     pmult3 : multiply two functions (three subroutine arguments)
c     p2mult : multiply with a function of the second argument only
c     p2comp : compose with a map from second argument to second argument
c     p1int  : convolute the first argument with a gaussian
c     pq     : reorganize the storage from (t,u) to (u,t)
c     qp     : reorganize the storage from (u,t) to (t,u)
c-----------------------------------------------------------------------
      subroutine scoefe(p,i,j,s)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),s(0:1)
      dimension s1(0:1),srhot(0:1),srhou(0:1)
c..on the unit disk a coefficient is bounded by the l-1 norm.
c..extract domain
      call srconst(p(0,0,nu3),srhou)
      call srconst(p(0,nt3,0),srhot)
c..extract a bound on a coefficient
      if(j.eq.0)then
	if(i.eq.0)then
	  call sequal(p,s)
	else
	  call spack(p,i,j,srhot,srhou,s)
	  call spack(p,nt+1,j,srhot,srhou,s1)
	  call sscalar(-s1(1),s1(1),s1)
	  call sadd(s1,s)
	  call sunpack(s,i,j,srhot,srhou)
	endif
      else
	if(i.eq.0)then
	  call spack(p,i,j,srhot,srhou,s)
	  call spack(p,i,nu+1,srhot,srhou,s1)
	  call sscalar(-s1(1),s1(1),s1)
	  call sadd(s1,s)
	  call sunpack(s,i,j,srhot,srhou)
	else
	  call spack(p,i,j,srhot,srhou,s)
	  call spack(p,nt+1,j,srhot,srhou,s1)
	  call sscalar(-s1(1),s1(1),s1)
	  call sadd(s1,s)
	  call spack(p,i,nu+1,srhot,srhou,s1)
	  call sscalar(-s1(1),s1(1),s1)
	  call sadd(s1,s)
	  call spack(p,nt+1,nu+1,srhot,srhou,s1)
	  call sscalar(-s1(1),s1(1),s1)
	  call sadd(s1,s)
	  call sunpack(s,i,j,srhot,srhou)
	endif
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine spack(p,i,j,srhot,srhou,sres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),
     *		srhot(0:1),srhou(0:1),sres(0:1)
      dimension s1(0:1),s2(0:1)
c..scale coefficient or error term to its size on unit disk.
      it=itgh(i)
      call sfak(it,s1)
      call spower(srhot,it,s2)
      call smult2(s2,s1)
      call spower(srhou,iugh(j),s2)
      call smult2(s2,s1)
      call smult3(p(0,i,j),s1,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine sunpack(s,i,j,srhot,srhou)
      implicit double precision(a-h,o-z)
      dimension s(0:1),srhot(0:1),srhou(0:1)
      dimension s1(0:1),s2(0:1)
c..inverse of spack.
      call sfak(i,s1)
      call spower(srhot,i,s2)
      call smult2(s2,s1)
      call spower(srhou,j,s2)
      call smult2(s2,s1)
      call squot(s,s1,s2)
      call sequal(s2,s)
      return
      end
c-----------------------------------------------------------------------
      subroutine scoefa(p,i,j,s)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),s(0:1)
      dimension s1(0:1)
c..extract a bound on a coefficient and add it
      call scoefe(p,i,j,s1)
      call sadd(s1,s)
      return
      end
c-----------------------------------------------------------------------
      subroutine scoefm(p,i,j,s)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),s(0:1)
      dimension s1(0:1)
c..extract a bound on a coefficient and subtract it
      call scoefe(p,i,j,s1)
      call ssubtr(s1,s)
      return
      end
c-----------------------------------------------------------------------
      subroutine pstzro(p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/rconst/rzero,rone,rtwo
      dimension p(0:1,0:nt3,0:nu3)
c..sets only coefficients and error terms equal to zero;
c  does nothing with domain.
      do 1 i=0,nu+2
      call tstzro(p(0,0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine p1zro(p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3)
c. sets only coefficients and error terms of F(t,0) equal to zero.
      call tstzro(p)
      return
      end
c-----------------------------------------------------------------------
      subroutine p2zro(p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3)
c. sets only coefficients and error terms of F(0,u) equal to zero.
      do 1 i=1,nu+2
      call siconst(0,p(0,0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine pequal(p,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),pres(0:1,0:nt3,0:nu3)
c..sets two functions of two variables equal.
c..transfer coefficients and error terms
      do 1 i=0,nu+2
      call tequal(p(0,0,i),pres(0,0,i))
    1 continue
c..transfer domain
      pres(0,0,nu3)=p(0,0,nu3)
      return
      end
c-----------------------------------------------------------------------
      subroutine p1eql(p,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),tres(0:1,0:nt3)
c..extracts the function F(t,0)
      call tequal(p,tres)
      return
      end
c-----------------------------------------------------------------------
      subroutine p2eql(p,ures)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),ures(0:1,0:nu3)
c..extracts the function F(0,u).
c..extract coefficients and error terms
      do 1 i=0,nu+2
      call sequal(p(0,0,i),ures(0,i))
    1 continue
c..extract domain
      ures(0,nu3)=p(0,0,nu3)
      return
      end
c-----------------------------------------------------------------------
      subroutine padd(p,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),pres(0:1,0:nt3,0:nu3)
c..add two functions of two variables
c..treat coefficients and error terms
      do 1 i=0,nu+2
      call tadd(p(0,0,i),pres(0,0,i))
    1 continue
c..adapt domains
      pres(0,0,nu3)=rmin2(p(0,0,nu3),pres(0,0,nu3))
      return
      end
c-----------------------------------------------------------------------
      subroutine p2der(p,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),pres(0:1,0:nt3,0:nu3)
      dimension q(0:1,0:nu3,0:nt3)
c..take the partial derivative with respect to the second variable.
c..change representation
      call pq(p,q)
c..treat coefficients and error terms
      do 1 i=0,nt+2
      call uder(q(0,0,i))
    1 continue
c..clean error terms
c..general
      call sabs1(q(0,0,nt+1))
      call sabs1(q(0,nu,nt+1))
c..higher order
      call sabs1(q(0,0,nt+2))
      call sabs1(q(0,nu,nt+2))
c..change representation; domains are ported from p to q to pres.
      call qp(q,pres)
      return
      end
c-----------------------------------------------------------------------
      subroutine psubtr(p,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),pres(0:1,0:nt3,0:nu3)
c..subtract two functions of two variables.
c..treat coefficients
      do 1 i=0,nu
      call tsubtr(p(0,0,i),pres(0,0,i))
    1 continue
c..treat error terms
      call tadd(p(0,0,nu+1),pres(0,0,nu+1))
      call tadd(p(0,0,nu+2),pres(0,0,nu+2))
c..adapt domains
      pres(0,0,nu3)=rmin2(p(0,0,nu3),pres(0,0,nu3))
      return
      end
c-----------------------------------------------------------------------
      subroutine psmult(s,p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension s(0:1),p(0:1,0:nt3,0:nu3)
      dimension s1(0:1)
c..multiply a function of two variables with a scalar.
c..treat coefficients
      do 1 i=0,nu
      call tsmlt2(s,p(0,0,i))
    1 continue
c..treat error terms
      call sabs2(s,s1)
      call tsmlt2(s1,p(0,0,nu+1))
      call tsmlt2(s1,p(0,0,nu+2))
      return
      end
c-----------------------------------------------------------------------
      subroutine psmadd(p,s,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),s(0:1),pres(0:1,0:nt3,0:nu3)
      dimension s1(0:1)
c..multiply a function of two variables with a scalar;
c  add the result to an other such function.
c..treat coefficients
      do 1 i=0,nu
      call tsmadd(p(0,0,i),s,pres(0,0,i))
    1 continue
c..treat error terms
      call sabs2(s,s1)
      call tsmadd(p(0,0,nu+1),s1,pres(0,0,nu+1))
      call tsmadd(p(0,0,nu+2),s1,pres(0,0,nu+2))
c..check domains
      pres(0,0,nu3)=rmin2(p(0,0,nu3),pres(0,0,nu3))
      return
      end
c-----------------------------------------------------------------------
      subroutine p1sdil(p,s,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),s(0:1),pres(0:1,0:nt3,0:nu3)
      dimension s1(0:1)
c..dilate the first argument of a function of two variables
c..treat coefficients
      do 1 i=0,nu
      call tsdil(p(0,0,i),s,pres(0,0,i))
    1 continue
c..treat error terms
      call sabs2(s,s1)
      call tsdil(p(0,0,nu+1),s1,pres(0,0,nu+1))
      call tsdil(p(0,0,nu+2),s1,pres(0,0,nu+2))
c..transfer domain of second variable
      pres(0,0,nu3)=p(0,0,nu3)
      return
      end
c-----------------------------------------------------------------------
      subroutine p2sdil(p,s)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),s(0:1)
      dimension q(0:1,0:nu3,0:nt3)
      dimension s1(0:1)
c..dilate the second argument of a function of two variables
c..change representation
      call pq(p,q)
c..treat coefficients
      do 1 i=0,nt
      call usdil(q(0,0,i),s)
    1 continue
c..treat error terms
      call sabs2(s,s1)
      call usdil(q(0,0,nt+1),s1)
      call usdil(q(0,0,nt+2),s1)
c..change representation; domains are portet from p to q to p
      call qp(q,p)
      return
      end
c-----------------------------------------------------------------------
      subroutine pmult2(p1,p2)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p1(0:1,0:nt3,0:nu3),p2(0:1,0:nt3,0:nu3)
      dimension p3(0:1,0:nt3,0:nu3),
     *		t1(0:1,0:nt3),t2(0:1,0:nt3),
     *		t1n(0:1,0:nt3),t2n(0:1,0:nt3),
     *		thigh(0:1,0:nt3),tgnrl(0:1,0:nt3)
      dimension s(0:1),srhou(0:1)
c..multiply two functions of two variables.
c..check domains
      p3(0,0,nu3)=rmin2(p1(0,0,nu3),p2(0,0,nu3))
      call srconst(p3(0,0,nu3),srhou)
c..treat coefficients
      do 1 j1=0,nu
      call tmult(p1(0,0,0),p2(0,0,j1),p3(0,0,j1))
      do 2 j=1,j1
      call tmadd(p1(0,0,j),p2(0,0,j1-j),p3(0,0,j1))
    2 continue
    1 continue
c..initialize thigh (first coefficients, then domain)
      call tstzro(thigh)
      thigh(0,nt3)=rmin2(p1(0,nt3,0),p2(0,nt3,0))
c..coefficients times coefficients to higher order
      do 3 j1=nu+1,2*nu
      do 4 j=j1-nu,nu
      call tmult(p1(0,0,j),p2(0,0,j1-j),t1)
      call spower(srhou,j1-nu-1,s)
      call tsmadd(t1,s,thigh)
    4 continue
    3 continue
      call tbound1(thigh)
c..coefficients times general or higher order to general or higher order
      call p1eql(p1,t1)
      call p1eql(p2,t2)
      do 5 i=1,nu-1
      call spower(srhou,i,s)
      call tsmadd(p1(0,0,i),s,t1)
      call tsmadd(p2(0,0,i),s,t2)
    5 continue
      call tequal(t1,t1n)
      call tequal(t2,t2n)
      call spower(srhou,nu,s)
      call tsmadd(p1(0,0,nu),s,t1n)
      call tsmadd(p2(0,0,nu),s,t2n)
      call tbound1(t1)
      call tbound1(t2)
      call tmult(t1,p2(0,0,nu+1),tgnrl)
      call tmadd(p1(0,0,nu+1),t2,tgnrl)
      call tbound1(t1n)
      call tbound1(t2n)
      call tmadd(t1n,p2(0,0,nu+2),thigh)
      call tmadd(p1(0,0,nu+2),t2n,thigh)
c..last coefficient times general to higher order
      call tbound2(p1(0,0,nu),t1)
      call tbound2(p2(0,0,nu),t2)
      call tmadd(t1,p2(0,0,nu+1),thigh)
      call tmadd(p1(0,0,nu+1),t2,thigh)
c..general times general to general
      call tmult(p1(0,0,nu+1),p2(0,0,nu+1),t1)
      call tsmadd(t1,srhou,tgnrl)
c..higher order times higher order to higer order
      call tmult(p1(0,0,nu+2),p2(0,0,nu+2),t1)
      call spower(srhou,nu+1,s)
      call tsmadd(t1,s,thigh)
c..store error terms
      call tequal(tgnrl,p3(0,0,nu+1))
      call tequal(thigh,p3(0,0,nu+2))
c..transfer result
      call pequal(p3,p2)
      return
      end
c-----------------------------------------------------------------------
      subroutine pmult3(p1,p2,pres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p1(0:1,0:nt3,0:nu3),p2(0:1,0:nt3,0:nu3),
     *		pres(0:1,0:nt3,0:nu3)
      call pequal(p2,pres)
      call pmult2(p1,pres)
      return
      end
c-----------------------------------------------------------------------
      subroutine p2mult(pres,u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension pres(0:1,0:nt3,0:nu3),u(0:1,0:nu3)
      dimension q(0:1,0:nu3,0:nt3),u1(0:1,0:nu3)
c..multiply a function of two variables with a function of the
c  second variable.
c..change representation
      call pq(pres,q)
c..treat coefficients
      do 1 i=0,nt
      call umult(u,q(0,0,i))
    1 continue
c..treat error terms
      call ubound2(u,u1)
      call umult(u1,q(0,0,nt+1))
      call umult(u1,q(0,0,nt+2))
c..domain is transfered from p1 to q to pres
      call qp(q,pres)
      return
      end
c-----------------------------------------------------------------------
      subroutine p2comp(pres,u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension pres(0:1,0:nt3,0:nu3),u(0:1,0:nu3)
      dimension q(0:1,0:nu3,0:nt3),u1(0:1,0:nu3)
c..compose a function of two variables in the second argument
c  with a function of the second variable.
      call pq(pres,q)
c..treat coefficients
      do 1 i=0,nt
      call ucomp(q(0,0,i),u)
    1 continue
c..treat error terms
      call ubound2(u,u1)
      call ucomp(q(0,0,nt+1),u1)
      call ucomp(q(0,0,nt+2),u1)
c..clean error terms
      call ubound1(q(0,0,nt+1))
      call ubound1(q(0,0,nt+2))
c..domains are transferred from p1 to q to pres
      call qp(q,pres)
      return
      end
c-----------------------------------------------------------------------
      subroutine p1int(pres,scov)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension pres(0:1,0:nt3,0:nu3),scov(0:1)
c..compute the convolution with a gaussian
c  with respect to the first argument.
c..treat coefficients and error terms
      do 1 i=0,nu+2
      call tint(pres(0,0,i),scov)
    1 continue
c..clean error terms
      call tbound1(pres(0,0,nu+1))
      call tbound1(pres(0,0,nu+2))
      return
      end
c-----------------------------------------------------------------------
      subroutine pq(p,q)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3),q(0:1,0:nu3,0:nt3)
c..change the representation from (t,u) to (u,t)
      do 1 j=0,nt+2
      do 2 i=0,nu+2
      call sequal(p(0,j,i),q(0,i,j))
    2 continue
      q(0,nu3,j)=p(0,0,nu3)
    1 continue
      q(0,0,nt3)=p(0,nt3,0)
      return
      end
c-----------------------------------------------------------------------
      subroutine qp(q,p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension q(0:1,0:nu3,0:nt3),p(0:1,0:nt3,0:nu3)
c..change the representation from (u,t) to (t,u)
      do 1 i=0,nu+2
      do 2 j=0,nt+2
      call sequal(q(0,i,j),p(0,j,i))
    2 continue
      p(0,nt3,i)=q(0,0,nt3)
    1 continue
      p(0,0,nu3)=q(0,nu3,0)
      return
      end
c-----------------------------------------------------------------------
c     OPERATIONS ON FUNCTIONS OF t
c     ============================
c
c     stnorm : compute norm of function of t
c     tstzro : set function zero
c     tequal : set two functions equal
c     tbound1: bounds on absolute value of coefficients (one argument)
c     tbound2: bounds on absolute value of coefficients (two arguments)
c     tadd   : add two functions
c     tsubtr : subtract two functions
c     tsmlt2 : multiply a function with a scalar (two arguments)
c     tsmlt3 : multiply a function with a scalar (three arguments)
c     tsmadd : multiply a function with a scalar and add it
c     tsdil  : dilate the argument
c     tmult  : multiply two functions of t
c     tmadd  : multiply two functions and add the result to a third one
c     tint   : convlute a function with a gaussian
c-----------------------------------------------------------------------
      subroutine stnorm(t,sres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),sres(0:1)
      dimension s(0:1),s1(0:1),s2(0:1),srhot(0:1)
c..computes upper and lower bound on the norm of t
c..extract growth rate
      call srconst(t(0,nt3),srhot)
c..initialize temporary variable
      call siconst(0,s)
c..add up coefficients and error term
      do 1 i=0,nt+2
      call spower(srhot,itgh(i),s1)
      call sfak(itgh(i),s2)
      call smult2(s2,s1)
      call sbound2(t(0,i),s2)
      call smadd(s1,s2,s)
    1 continue
c..transfer result
      call sbound2(s,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine tstzro(t)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
c..sets coefficients and error terms of t equal to zero;
c  domain needs not be defined.
      do 1 i=0,nt+2
      call siconst(0,t(0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine tequal(t,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),tres(0:1,0:nt3)
c..sets two functions of t equal.
c..transfer coefficients and error terms
      do 1 i=0,nt+2
      call sequal(t(0,i),tres(0,i))
    1 continue
c..transfer domain
      tres(0,nt3)=t(0,nt3)
      return
      end
c-----------------------------------------------------------------------
      subroutine tbound1(t)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
c..replace the coefficients of t by bounds on their absolute values.
      do 1 i=0,nt+2
      call sbound1(t(0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine tbound2(t,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),tres(0:1,0:nt3)
c..computes bounds on the absolute values of the coefficients of t.
      call tequal(t,tres)
      call tbound1(tres)
      return
      end
c-----------------------------------------------------------------------
      subroutine tadd(t,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),tres(0:1,0:nt3)
c..add two functions of t
c..add coefficients and error terms
      do 1 i=0,nt+2
      call sadd(t(0,i),tres(0,i))
    1 continue
c..adapt domains
      tres(0,nt3)=rmin2(t(0,nt3),tres(0,nt3))
      return
      end
c-----------------------------------------------------------------------
      subroutine tsubtr(t,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),tres(0:1,0:nt3)
c..computes the difference of two functions
c..subtract coefficients
      do 1 i=0,nt
      call ssubtr(t(0,i),tres(0,i))
    1 continue
c..add errror terms
      call sadd(t(0,nt+1),tres(0,nt+1))
      call sadd(t(0,nt+2),tres(0,nt+2))
c..adapt domains
      tres(0,nt3)=rmin2(t(0,nt3),tres(0,nt3))
      return
      end
c-----------------------------------------------------------------------
      subroutine tsmlt2(s,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension s(0:1),tres(0:1,0:nt3)
      dimension s1(0:1)
c..multiplies t with the scalar s
c..treat coefficients
      do 1 i=0,nt
      call smult2(s,tres(0,i))
    1 continue
c..treat error terms
      call sbound2(s,s1)
      call smult2(s1,tres(0,nt+1))
      call smult2(s1,tres(0,nt+2))
      return
      end
c-----------------------------------------------------------------------
      subroutine tsmlt3(t,s,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),s(0:1),tres(0:1,0:nt3)
      call tequal(t,tres)
      call tsmlt2(s,tres)
      return
      end
c-----------------------------------------------------------------------
      subroutine tsmadd(t,s,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),s(0:1),tres(0:1,0:nt3)
      dimension s1(0:1)
c..multiply t with s and add the result to tres.
c..treat coefficients
      do 1 i=0,nt
      call smadd(s,t(0,i),tres(0,i))
    1 continue
c..treat error terms
      call sbound2(s,s1)
      call smadd(s1,t(0,nt+1),tres(0,nt+1))
      call smadd(s1,t(0,nt+2),tres(0,nt+2))
c..adapt domains
      tres(0,nt3)=rmin2(t(0,nt3),tres(0,nt3))
      return
      end
c-----------------------------------------------------------------------
      subroutine tsdil(t,s,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3),s(0:1),tres(0:1,0:nt3)
      dimension s1(0:1)
c..dilate the function t with the scalar s.
c..initialize temporary variable
      call siconst(1,s1)
c..dilate coefficients
      do 1 i=0,nt
      call smult3(t(0,i),s1,tres(0,i))
      call smult2(s,s1)
    1 continue
c..treat error terms
      call sbound1(s1)
      call smult3(t(0,nt+2),s1,tres(0,nt+2))
      call sbound2(s,s1)
      call smult3(t(0,nt+1),s1,tres(0,nt+1))
c..adjust domain
      tres(0,nt3)=rdquot(t(0,nt3),s1(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine tmult(t1,t2,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension t1(0:1,0:nt3),t2(0:1,0:nt3),
     *		tres(0:1,0:nt3)
      dimension t3(0:1,0:nt3),
     *		s(0:1),s1(0:1),s2(0:1),si(0:1),
     *		s1n(0:1),s2n(0:1),
     *		shigh(0:1),sgnrl(0:1),srhot(0:1)
      dimension s3(0:1)
c..multiply two functions of t
c..check domains
      t3(0,nt3)=rdinv(rusum(ruinv(t1(0,nt3)),ruinv(t2(0,nt3))))
      call srconst(t3(0,nt3),srhot)
c..treat coefficients
      do 1 i1=0,nt
      call siconst(0,t3(0,i1))
      do 2 i=0,i1
      call smadd(t1(0,i),t2(0,i1-i),t3(0,i1))
    2 continue
    1 continue
c..coefficients times coefficients to higher order
      call siconst(0,s3)
      do 3 i1=nt+1,2*nt
      call siconst(0,s)
      do 4 i=i1-nt,nt
      call smadd(t1(0,i),t2(0,i1-i),s)
    4 continue
      call spower(srhot,i1-nt-1,s1)
      call sfak(i1,s2)
      call smult2(s2,s1)
      call smadd(s1,s,s3)
    3 continue
      call sfak(nt+1,s1)
      call squot(s3,s1,shigh)
      call sbound1(shigh)
c..coefficients times error terms to error terms
      call sequal(t1,s1)
      call sequal(t2,s2)
      do 5 i=1,nt-1
      call spower(srhot,i,s)
      call sfak(i,si)
      call smult2(si,s)
      call smadd(t1(0,i),s,s1)
      call smadd(t2(0,i),s,s2)
    5 continue
      call sequal(s1,s1n)
      call sequal(s2,s2n)
      call spower(srhot,nt,s)
      call sfak(nt,si)
      call smult2(si,s)
      call smadd(t1(0,nt),s,s1n)
      call smadd(t2(0,nt),s,s2n)
      call sbound1(s1)
      call sbound1(s2)
      call smult3(s1,t2(0,nt+1),sgnrl)
      call smadd(t1(0,nt+1),s2,sgnrl)
      call sbound1(s1n)
      call sbound1(s2n)
      call smadd(s1n,t2(0,nt+2),shigh)
      call smadd(t1(0,nt+2),s2n,shigh)
c..last coefficient  times general to higher order
      call siconst(nt+1,si)
      call squot(t1(0,nt),si,s1)
      call squot(t2(0,nt),si,s2)
      call sbound1(s1)
      call sbound1(s2)
      call smadd(s1,t2(0,nt+1),shigh)
      call smadd(t1(0,nt+1),s2,shigh)
c..general times general to general
      call smult3(t1(0,nt+1),t2(0,nt+1),s1)
      call smadd(s1,srhot,sgnrl)
c..higher order times higher order to higer order
      call smult3(t1(0,nt+2),t2(0,nt+2),s1)
      call spower(srhot,nt+1,s)
      call sfak(nt+1,si)
      call smult2(si,s)
      call smadd(s1,s,shigh)
c..store error terms
      call sequal(sgnrl,t3(0,nt+1))
      call sequal(shigh,t3(0,nt+2))
c..transfer result
      call tequal(t3,tres)
      return
      end
c-----------------------------------------------------------------------
      subroutine tmadd(t1,t2,tres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t1(0:1,0:nt3),t2(0:1,0:nt3),
     *		tres(0:1,0:nt3)
      dimension t3(0:1,0:nt3)
c..multiply t1 with t2 and add the result to tres.
      call tmult(t1,t2,t3)
      call tadd(t3,tres)
      return
      end
c-----------------------------------------------------------------------
      subroutine tint(t,scov)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/rconst/rzero,rone,rtwo
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      common/output/iout
      dimension t(0:1,0:nt3),scov(0:1)
      dimension t1(0:1,0:nt3)
      dimension s(0:1),sfact(0:1)
      dimension st1(0:1),st2(0:1),si(0:1),srho(0:1)
      dimension s3(0:1)
c..compute the convolution of t with a Gaussian with covariance
c  scov, unless the domains do not admit this.
      rrho=t(0,nt3)
      ru2cov=ruprod(rtwo,scov(1))
c..transfer domain
      rrhohat=rddiff(rrho,ru2cov)
      t1(0,nt3)=rrhohat
c..check domain
      if(rrhohat.lt.rzero)then
        write(iout,*)'error in integration, function grows to fast',
     *	t(0,nt3)
      endif
c..treat coefficients
      call tstzro(t1)
      call squot(scov,stwo,s)
      do 1 i=0,nt
      i2=2*i+2
      i1=i2-1
      call siconst(1,sfact)
      call siconst(0,st1)
      do 2 j=0,nt-i
      call smult3(t(0,i+j),sfact,st2)
      call sadd(st2,st1)
      j2=2*j
      call smult2(s,sfact)
      call siconst((i2+j2)*(i1+j2),si)
      call smult2(si,sfact)
      call siconst(j+1,si)
      call squot(sfact,si,s3)
      call sequal(s3,sfact)
    2 continue
      call sequal(st1,t1(0,i))
    1 continue
c..treat error terms
      call srconst(rrho,srho)
c..general to general
      call sequal(t(0,nt+1),t1(0,nt+1))
c.. general to first coefficient
      call smult3(t(0,nt+1),srho,st1)
      call sscalar(-st1(1),st1(1),st2)
      call sadd(st2,t1)
c..higher order to general
      call spower(srho,nt,st2)
      call sfak(nt+1,sfact)
      call smult2(sfact,st2)
      call smult2(t(0,nt+2),st2)
      call sbound1(st2)
      call sadd(st2,t1(0,nt+1))
c..higher order to first coefficient
      call smult2(srho,st2)
      call sscalar(-st2(1),st2(1),st1)
      call sadd(st1,t1)
c..higher order is zero
      call siconst(0,t1(0,nt+2))
c..transfer the result
      call tequal(t1,t)
      return
      end
c-----------------------------------------------------------------------
c     OPERATIONS ON FUNCTIONS OF u
c     ============================
c
c     sunorm : compute norm of function of u
c     ustzro : set function zero
c     uequal : set two functions equal
c     ubound1: compute bounds on the absolute value of coefficients
c     ubound2: compute bounds on the absolute value of coefficients
c     usmult : multiply function with scalar
c     uder   : take derivative
c     usdil  : dilate the argument
c     umult  : multiply two functions
c     ucomp  : compose two functions
c     utkinv : compute one over a function
c     ustlin : set a function equal to a linear polynom
c-----------------------------------------------------------------------
      subroutine sunorm(u,sres)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3),sres(0:1)
      dimension s(0:1),s1(0:1),s2(0:1),srhou(0:1)
c..computes upper and lower bound on the norm of u
c..extract domain
      call srconst(u(0,nu3),srhou)
c..initialize temporary variables
      call siconst(1,s2)
      call siconst(0,s)
c..add up coefficients
      do 1 i=0,nu
      call sbound2(u(0,i),s1)
      call smadd(s1,s2,s)
      call smult2(srhou,s2)
    1 continue
c..add contribution of higher order term
      call smadd(u(0,nu+2),s2,s)
c..add contribution of general term
      call smadd(u(0,nu+1),srhou,s)
c..return value
      call sbound2(s,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine ustzro(u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3)
c..sets coefficients and error terms of u equal to zero;
c  domain needs not be defined.
      do 1 i=0,nu+2
      call siconst(0,u(0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine uequal(u,ures)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3),ures(0:1,0:nu3)
c..sets two functions equal.
c..transfer coefficients and error terms
      do 1 i=0,nu+2
      call sequal(u(0,i),ures(0,i))
    1 continue
c..transfer domain
      ures(0,nu3)=u(0,nu3)
      return
      end
c-----------------------------------------------------------------------
      subroutine ubound1(u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3)
c..replaces the coefficients of u by bounds on their absolute values.
      do 1 i=0,nu+2
      call sbound1(u(0,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine ubound2(u,ures)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3),ures(0:1,0:nu3)
c..computes bounds on the absolute values of the coefficients of u.
      call uequal(u,ures)
      call ubound1(ures)
      return
      end
c-----------------------------------------------------------------------
      subroutine usmult(u,s,ures)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3),s(0:1),ures(0:1,0:nu3)
      dimension s1(0:1)
c..multiplies u with the scalar s
      call uequal(u,ures)
c..treat coefficients
      do 1 i=0,nu
      call smult2(s,ures(0,i))
    1 continue
c..treat error terms
      call sbound2(s,s1)
      call smult2(s1,ures(0,nu+1))
      call smult2(s1,ures(0,nu+2))
      return
      end
c-----------------------------------------------------------------------
      subroutine uder(u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension u(0:1,0:nu3)
      dimension u1(0:1,0:nu3),
     *		s(0:1),si(0:1),srhou(0:1),
     *		s1(0:1),s2(0:1),sg(0:1),sl(0:1),sh(0:1)
c..takes the drivative of u
c..treat coefficients
      do 1 i=1,nu
      call siconst(i,si)
      call smult3(u(0,i),si,u1(0,i-1))
    1 continue
c..u(0,nu) is treated below
c..treat error terms (the domain is restricted by r=s(0)
      call srconst(u(0,nu3),srhou)
c..compute error terms
c..contribution from general to constant term
      call smult3(s4over,u(0,nu+1),s1)
      call sscalar(-s1(1),s1(1),s2)
      call sadd(s2,u1)
c..contribution from general to general
      call squot(s4over,srhou,s1)
      call squot(s1,semp25,sg)
      call smult3(sg,u(0,nu+1),s1)
      call sbound2(s1,u1(0,nu+1))
c..contribution from higher order to last coefficient
      if(nu+1.gt.4)then
	call siconst(nu+1,sl)
      else
	call spower(semp25,nu,s1)
	call squot(s4over,s1,sl)
      endif
      call smult3(u(0,nu+2),sl,s1)
      call sscalar(-s1(1),s1(1),s2)
      call sbound2(s2,u1(0,nu))
c..contribution from higher order to higer order
      if(nu+2.gt.4)then
	call siconst(nu+2,s)
      else
	call spower(semp25,nu+1,s1)
	call squot(s4over,s1,s)
      endif
      call squot(s,srhou,sh)
      call smult3(u(0,nu+2),sh,s1)
      call sbound2(s1,u1(0,nu+2))
c..adapt domain
      u1(0,nu3)=rdprod(u(0,nu3),semp25(1))
c..transfer result
      call uequal(u1,u)
      return
      end
c-----------------------------------------------------------------------
      subroutine usdil(u,s)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3),s(0:1)
      dimension s1(0:1)
c..dilate the function u with the scalar s.
c..initialize temporary variable
      call siconst(1,s1)
c..dilate coefficients
      do 1 i=0,nu
      call smult2(s1,u(0,i))
      call smult2(s,s1)
    1 continue
c..treat error terms
      call sbound1(s1)
      call smult2(s1,u(0,nu+2))
      call sbound2(s,s1)
      call smult2(s1,u(0,nu+1))
c..adjust domain
      u(0,nu3)=rdquot(u(0,nu3),s1(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine umult(u1,u2)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u1(0:1,0:nu3),u2(0:1,0:nu3)
      dimension u3(0:1,0:nu3),
     *		s(0:1),s1(0:1),s2(0:1),
     *		s1n(0:1),s2n(0:1),
     *		shigh(0:1),sgnrl(0:1),srhou(0:1)
c..multiply two functions of u
c..check domains
      u3(0,nu3)=rmin2(u1(0,nu3),u2(0,nu3))
      call srconst(u3(0,nu3),srhou)
c..treat coefficients
      do 1 i1=0,nu
      call siconst(0,u3(0,i1))
      do 2 i=0,i1
      call smadd(u1(0,i),u2(0,i1-i),u3(0,i1))
    2 continue
    1 continue
c..coefficients times coefficients to higher order
      call siconst(0,shigh)
      do 3 i1=nu+1,2*nu
      call spower(srhou,i1-nu-1,s)
      do 4 i=i1-nu,nu
      call smult3(u1(0,i),u2(0,i1-i),s1)
      call smadd(s1,s,shigh)
    4 continue
    3 continue
      call sbound1(shigh)
c..coefficients times error terms to error terms
      call sequal(u1,s1)
      call sequal(u2,s2)
      do 5 i=1,nu-1
      call spower(srhou,i,s)
      call smadd(u1(0,i),s,s1)
      call smadd(u2(0,i),s,s2)
    5 continue
      call sequal(s1,s1n)
      call sequal(s2,s2n)
      call spower(srhou,nu,s)
      call smadd(u1(0,nu),s,s1n)
      call smadd(u2(0,nu),s,s2n)
      call sbound1(s1)
      call sbound1(s2)
      call smult3(s1,u2(0,nu+1),sgnrl)
      call smadd(u1(0,nu+1),s2,sgnrl)
      call sbound1(s1n)
      call sbound1(s2n)
      call smadd(s1n,u2(0,nu+2),shigh)
      call smadd(u1(0,nu+2),s2n,shigh)
c..last coefficient times general to higher order
      call sbound2(u1(0,nu),s1)
      call sbound2(u2(0,nu),s2)
      call smadd(s1,u2(0,nu+1),shigh)
      call smadd(u1(0,nu+1),s2,shigh)
c..general times general to general
      call smult3(u1(0,nu+1),u2(0,nu+1),s1)
      call smadd(s1,srhou,sgnrl)
c..higher order times higher order to higer order
      call smult3(u1(0,nu+2),u2(0,nu+2),s1)
      call spower(srhou,nu+1,s)
      call smadd(s1,s,shigh)
c..store error terms
      call sequal(sgnrl,u3(0,nu+1))
      call sequal(shigh,u3(0,nu+2))
c..transfer result
      call uequal(u3,u2)
      return
      end
c-----------------------------------------------------------------------
      subroutine ucomp(ures,uarg)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/output/iout
      dimension ures(0:1,0:nu3),uarg(0:1,0:nu3)
      dimension u1(0:1,0:nu3)
      dimension s(0:1),s1(0:1),s2(0:1)
c..compute the composition ures(uarg(.)) if no domain violation occurs.
c..check that the range of the argument is smaller than
c  the domain of analyticity of the function into which the
c  substitution is performed.
      call sunorm(uarg,s)
      if(s(1).gt.ures(0,nu3))then
        write(iout,*)'error in ucomp',ures(0,nu3),s(1)
      endif
c..transfer domain (result is analytic on domain of argument)
      u1(0,nu3)=uarg(0,nu3)
c..initializae coefficients and error terms
      call ustzro(u1)
c..compose into polynomial part of ures
      call sequal(ures(0,nu),u1)
      do 1 i=1,nu
      call umult(uarg,u1)
      call sadd(ures(0,nu-i),u1)
    1 continue
c..treat error terms
c..contribution  of general to first coefficient
      call smult3(s,ures(0,nu+1),s1)
      call sscalar(-s1(1),s1(1),s2)
      call sadd(s2,u1)
c..contribution to general from general
      call smult2(s,s1)
      call sbound1(s1)
      call sadd(s1,u1(0,nu+1))
c..contribution to general from higher order
      call spower(s,nu,s1)
      call sbound1(s1)
      call smadd(s1,ures(0,nu+2),u1(0,nu+1))
c..transfer result
      call uequal(u1,ures)
      return
      end
c-----------------------------------------------------------------------
      subroutine utkinv(u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension u(0:1,0:nu3)
      dimension uinv(0:1,0:nu3)
      dimension ut1(0:1,0:nu3)
      dimension scof0i(0:1),s1(0:1),s2(0:1)
      print*,'in utkinv'
      call ucheck(u)
c..compute reciprocal of u1 by using that
c  1/u(.)=1/(u(0)+(u(.)-u(0)))=(1/u(0))*(1/(1+(u(.)-u(0))/u(0)));
c  for the last term one uses the composition of 1/(1+z) with
c  (u(.)-u(0))/u(0).
c..compute 1/(1+z)
c..fix coefficients of 1/(1+z)
      call siconst(1,uinv)
      do 1 i=1,nu
      call sneg2(uinv(0,i-1),uinv(0,i))
    1 continue
c..error terms for uinv are fixed later to their optimal values
c..compute argument for composition
      call sequal(u,scof0i)
      call sinv1(scof0i)
      call usmult(u,scof0i,ut1)
      call siconst(0,ut1)
c..choose domain for 1/(1+z)
      call sunorm(ut1,s1)
      uinv(0,nu3)=s1(1)
c..general error term for 1/(1+z) is zero
      call siconst(0,uinv(0,nu+1))
c..compute higher order error term for 1/(1+z)
      call spower(s1,nu+1,s2)
      call sneg1(s1)
      call sadd(sone,s1)
      call squot(s2,s1,uinv(0,nu+2))
c..compose (the composition is well defined by construction)
      call ucomp(uinv,ut1)
      call usmult(uinv,scof0i,u)
      print*,'in utkinv'
      call ucheck(u)
      return
      end
c-----------------------------------------------------------------------
      subroutine ustlin(u,s0,su,rrho)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3),s0(0:1),su(0:1)
c..makes a linear function on prescribed domain
c..set domain
      u(0,nu3)=rrho
c..set coefficients and error terms
      call ustzro(u)
      call sequal(s0,u(0,0))
      call sequal(su,u(0,1))
      return
      end
c-----------------------------------------------------------------------
c     5.8. OPERATIONS INVOLVING MATRICES
c     ==================================
c
c     amati  : invers of a matrix		   (computer arithmetic)
c     amatm  : product of two matricesc 	   (computer arithmetic)
c     amate  : set matrices equal		   (computer arithmetic)
c     amatn  : replace a matrix by its negative    (computer arithmetic)
c     amatd2 : add two to the diagonal of a matrix (computer arithmetic)
c
c     matas  : convert a numerical matrix to a scalar matrix
c     matsa  : convert a scalar matrix to a numerical matrix
c     mats   : make a matrix of intervalls
c
c     mati  : invers of a matrix
c     mata  : add two matrices
c     matm  : product of two matrices
c     mate  : set matrices equal
c     matn  : replce a matrix by its negative
c     matd1 : add two to the diagonal of a matrix
c
c     lin   : checks if one matrix is included in an other
c-----------------------------------------------------------------------
      subroutine amati(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      common/rconst/rzero,rone,rtwo
      dimension aa(1:imat,1:imat),bb(1:imat,1:imat)
      dimension cc(1:imat,1:imat),dd(1:imat,1:imat)
c..compute the numerical inverse of a matrix by the Pan/Reiff-Method.
c..compute l1-norm of rows.
      r1=0
      do 1 i=1,imat
      r=0
      do 2 j=1,imat
      r=r+rabs(aa(i,j))
    2 continue
      r1=rmax2(r,r1)
    1 continue
c..compute l1-norm of colums.
      r2=0
      do 3 j=1,imat
      r=0
      do 4 i=1,imat
      r=r+rabs(aa(i,j))
    4 continue
      r2=rmax2(r,r2)
    3 continue
c..compute normalized transposed matrix as the guess of order zero.
      do 5 i=1,imat
      do 6 j=1,imat
      dd(i,j)=aa(j,i)/(r1*r2)
    6 continue
    5 continue
c..iterate the contracting map
      do 7 i=1,100
      call amatm(dd,aa,cc)
      call amatn(cc)
      call amatd2(cc)
      call amatm(cc,dd,bb)
      call amate(bb,dd)
    7 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine amatm(aa,bb,cc)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(1:imat,1:imat),bb(1:imat,1:imat),
     *		cc(1:imat,1:imat)
c..multoply two numerical matrices.
      do 1 i=1,imat
      do 2 j=1,imat
      r=0
      do 3 k=1,imat
      r=r+aa(i,k)*bb(k,j)
    3 continue
      cc(i,j)=r
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine amate(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(1:imat,1:imat),bb(1:imat,1:imat)
c..set the numerical matrix bb equal to the numerical matrix aa.
      do 1 i=1,imat
      do 2 j=1,imat
      bb(i,j)=aa(i,j)
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine amatn(aa)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(1:imat,1:imat)
c..replace a matrix by its inverse
      do 1 i=1,imat
      do 2 j=1,imat
      aa(i,j)=-aa(i,j)
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine amatd2(aa)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(1:imat,1:imat)
c..add two in the diagonal of the numerical matrix aa.
      do 1 i=1,imat
      aa(i,i)=2+aa(i,i)
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine matas(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(1:imat,1:imat),bb(0:1,1:imat,1:imat)
c..make a matrix of scalars from a numerical matrix.
      do 1 j=1,imat
      do 2 i=1,imat
      call srconst(aa(i,j),bb(0,i,j))
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine matsa(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat),bb(1:imat,1:imat)
c..make a numerical matrix from a matrix of scalars.
      do 1 j=1,imat
      do 2 i=1,imat
      bb(i,j)=(aa(0,i,j)+aa(1,i,j))/2
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine mats(r,aa)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat)
c..make a matrix of scalars (-r,r).
      do 1 j=1,imat
      do 2 i=1,imat
      call sscalar(-r,r,aa(0,i,j))
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine mata(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat),bb(0:1,1:imat,1:imat)
c..add the matrix aa to the matrix bb.
      do 1 i=1,imat
      do 2 j=1,imat
      call sadd(aa(0,i,j),bb(0,i,j))
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine mati(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      common/rconst/rzero,rone,rtwo
      common/output/iout
      dimension aa(0:1,1:imat,1:imat),bb(0:1,1:imat,1:imat)
      dimension bb1(0:1,1:imat,1:imat),bb2(0:1,1:imat,1:imat)
      dimension rr(0:1,1:imat,1:imat)
      dimension eps(0:1,1:imat,1:imat)
      dimension cc(0:1,1:imat,1:imat),zz(0:1,1:imat,1:imat)
      dimension tt(0:1,1:imat,1:imat)
      dimension raa(1:imat,1:imat),rbb(1:imat,1:imat)
      logical lin
c..computes a bound on the inverse of a matrix given a "best possible"
c  numerical guess.
c..find a good guess for inverse
      call matsa(aa,raa)
      call amati(raa,rbb)
      call matas(rbb,rr)
c..prove existence of inverse
      call matm(rr,aa,cc)
      call matn(cc)
      call matd1(cc)
      call matm(aa,rr,tt)
      call matn(tt)
      call matd1(tt)
      call matm(rr,tt,zz)
      r=rup(rzero)
      call mats(r,eps)
      call mate(eps,bb1)
    1 continue
      call matm(cc,bb1,bb2)
      call mata(zz,bb2)
      if(.not.lin(bb2,bb1))then
        write(iout,*)'in mati: not yet a contraction.'
	call mata(eps,bb2)
	call mate(bb2,bb1)
	goto 1
      else
	call mata(rr,bb1)
      endif
      call mate(bb1,bb)
      return
      end
c-----------------------------------------------------------------------
      subroutine matm(aa,bb,cc)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat),bb(0:1,1:imat,1:imat),
     *		cc(0:1,1:imat,1:imat)
      dimension s(0:1),s1(0:1)
c..computes the product of two matrices.
      do 1 i=1,imat
      do 2 j=1,imat
      call siconst(0,s)
      do 3 k=1,imat
      call smult3(aa(0,i,k),bb(0,k,j),s1)
      call sadd(s1,s)
    3 continue
      call sequal(s,cc(0,i,j))
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine mate(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat),bb(0:1,1:imat,1:imat)
c..sets the matrix bb equal to the matrix aa.
      do 1 i=1,imat
      do 2 j=1,imat
      call sequal(aa(0,i,j),bb(0,i,j))
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine matn(aa)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat)
c..stores the negative of the matrix aa in bb.
      do 1 i=1,imat
      do 2 j=1,imat
      call sneg1(aa(0,i,j))
    2 continue
    1 continue
      return
      end
c-----------------------------------------------------------------------
      subroutine matd1(aa)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension aa(0:1,1:imat,1:imat)
c..adds one in the diagonal of the matrix aa.
      do 1 i=1,imat
      call sadd(sone,aa(0,i,i))
    1 continue
      return
      end
c-----------------------------------------------------------------------
      logical function lin(aa,bb)
      implicit double precision(a-h,o-z)
      parameter(imat=9)
      dimension aa(0:1,1:imat,1:imat),bb(0:1,1:imat,1:imat)
      logical l
c..checks if all the elements of the matrix aa are included in
c  the elements of the matrix bb.
      l=.true.
      do 1 j=1,imat
      do 2 i=1,imat
      l=l.and.((aa(0,i,j).gt.bb(0,i,j))
     *	       .and.(aa(1,i,j).lt.bb(1,i,j)))
    2 continue
    1 continue
      lin=l
      return
      end
c-----------------------------------------------------------------------
c
c     5.7. ROUTINES ACTING ON SCALARS
c     ===============================
c
c     lsequal :  compare two scalars
c     rmaxabs :  largest distance from zero
c     rminabs :  smallest distance from zero
c     siconst :  convert integer to scalar
c     srconst :  convert rep to scalar
c     sscalar :  interval with given boundary
c     sabs1   :  absolute value of scalar (1 argument)
c     sabs2   :  absolute value of scalar (2 arguments)
c     sbound1 :  bound on absolute value of scalar (1 argument)
c     sbound2 :  bound on absolute value of scalar (2 arguments)
c     sneg1   :  replace scalar by its negative
c     sneg2   :  compute the negative of a scalar
c     sinv1   :  replace a scalar by its inverse
c     sinv2   :  compute the inverse of a scalar
c     sequal  :  copies scalar
c     sinter  :  intersection of two scalars
c     sadd    :  sum of two scalars
c     ssubtr  :  difference of two scalars
c     smult3  :  product of two scalars (three arguments)
c     smult2  :  product of two scalars (two arguments)
c     smadd   :  multiplies to scalars and adds the result to a third one
c     squot   :  quotient of two scalars
c     sfak    :  Euler Gamma of k+1/2
c     spower  :  scalar to the power k
c     ssqrt   :  replace scalar by its square root
c     slog    :  natural log
c     sexp    :  exponential function
c-----------------------------------------------------------------------
      logical function lsequal(s1,s2)
      implicit double precision(a-h,o-z)
      dimension s1(0:1),s2(0:1)
c..answer is "true" if the two scalars are equal.
      lsequal=(s1(0).eq.s2(0)).and.(s1(1).eq.s2(1))
      return
      end
c-----------------------------------------------------------------------
      double precision function rmaxabs(s)
      implicit double precision(a-h,o-z)
      dimension s(0:1)
c..compute the maximum of the absolute values in s.
      rmaxabs=rmax2(-s(0),s(1))
      return
      end
c-----------------------------------------------------------------------
      double precision function rminabs(s)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      dimension s(0:1)
c..compute the minimum of the absolute values in s.
      rminabs=rmax3(rzero,s(0),-s(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine siconst(i,s)
      implicit double precision(a-h,o-z)
      dimension s(0:1)
c..produce a scalar from the integer i.
      s(0)=i
      s(1)=i
      return
      end
c-----------------------------------------------------------------------
      subroutine srconst(r,s)
      implicit double precision(a-h,o-z)
      dimension s(0:1)
c..produce a scalar of with zero from the representable number r.
      s(0)=r
      s(1)=r
      return
      end
c-----------------------------------------------------------------------
      subroutine sscalar(rl,ru,s)
      implicit double precision(a-h,o-z)
      dimension s(0:1)
c..make scalar from two reps
      s(0)=rl
      s(1)=ru
      return
      end
c-----------------------------------------------------------------------
      subroutine sabs2(s,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      dimension s(0:1),sres(0:1)
c..compute the absolute value of a scalar.
      r=s(0)
      sres(0)=rmax3(rzero,r,-s(1))
      sres(1)=rmax2(-r,s(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine sabs1(s)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      dimension s(0:1)
c..replaces a scalar by its absolute value
      r=s(0)
      s(0)=rmax3(rzero,r,-s(1))
      s(1)=rmax2(-r,s(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine sbound2(s,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      dimension s(0:1),sres(0:1)
c..computes a bound on the absolute value.
      sres(0)=rzero
      sres(1)=rmax2(-s(0),s(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine sbound1(s)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      dimension s(0:1)
c..computes a bound on the absolute value.
      r=s(0)
      s(0)=rzero
      s(1)=rmax2(-r,s(1))
      return
      end
c-----------------------------------------------------------------------
      subroutine sneg2(s,sres)
      implicit double precision(a-h,o-z)
      dimension s(0:1),sres(0:1)
c..compute the negative of a scalar.
      r=-s(1)
      sres(1)=-s(0)
      sres(0)=r
      return
      end
c-----------------------------------------------------------------------
      subroutine sneg1(s)
      implicit double precision(a-h,o-z)
      dimension s(0:1)
c..replace a scalar by its negative
      r=-s(1)
      s(1)=-s(0)
      s(0)=r
      return
      end
c-----------------------------------------------------------------------
      subroutine sinv1(s)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/output/iout
      dimension s(0:1)
c..replace a scalar by its inverse.
      if(s(0).le.rzero.and.rzero.le.s(1))then
        write(iout,*)'error in sinv',s(0),s(1)
      else
	rl=rdquot(rone,s(1))
	ru=ruquot(rone,s(0))
	call sscalar(rl,ru,s)
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine sinv2(s,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/output/iout
      dimension s(0:1),sres(0:1)
c..compute the inverse of a scalar.
      if(s(0).le.rzero.and.rzero.le.s(1))then
        write(iout,*)'error in sinv',s(0),s(1)
	call sequal(s,sres)
      else
	call sscalar(rdquot(rone,s(1)),ruquot(rone,s(0)),sres)
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine sequal(s,sres)
      implicit double precision(a-h,o-z)
      dimension s(0:1),sres(0:1)
c..set two scalars equal.
      sres(0)=s(0)
      sres(1)=s(1)
      return
      end
c-----------------------------------------------------------------------
      subroutine sinter(s1,s2,sres)
      implicit double precision(a-h,o-z)
      common/output/iout
      dimension s1(0:1),s2(0:1),sres(0:1)
c..compute the intersection of two scalars.
      rl=rmax2(s1(0),s2(0))
      ru=rmin2(s1(1),s2(1))
      if(rl.gt.ru)then
        write(iout,*)'error: no intersection, s1,s2='
	write(iout,*)s1(0),s1(1)
	write(iout,*)s2(0),s2(1)
	call sequal(s1,sres)
      else
	call sscalar(rl,ru,sres)
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine sadd(s1,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      dimension s1(0:1),sres(0:1)
c..computes bounds on the sum of two scalars.
c..extract bounds of arguments.
      r1l=s1(0)
      r1u=s1(1)
      r2l=sres(0)
      r2u=sres(1)
c..compute bounds for sum; treat special cases specially.
      if(r1l.eq.rzero)then
	if(r1u.eq.rzero)then
	  return
	elseif(r2u.eq.rzero)then
	  call sscalar(r2l,r1u,sres)
	else
	  call sscalar(r2l,rusum(r1u,r2u),sres)
	endif
      elseif(r2l.eq.rzero)then
	if(r2u.eq.rzero)then
	  call sequal(s1,sres)
	elseif(r1u.eq.rzero)then
	  call sscalar(r1l,r2u,sres)
	else
	  call sscalar(r1l,rusum(r1u,r2u),sres)
	endif
      else
	if(r1u.eq.rzero)then
	  call sscalar(rdsum(r1l,r2l),r2u,sres)
	elseif(r2u.eq.rzero)then
	  call sscalar(rdsum(r1l,r2l),r1u,sres)
	else
	  call sscalar(rdsum(r1l,r2l),rusum(r1u,r2u),sres)
	endif
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine ssubtr(s1,sres)
      implicit double precision(a-h,o-z)
      dimension s1(0:1),sres(0:1)
      dimension s(0:1)
c..computes bounds on the difference of two scalars.
      call sneg2(s1,s)
      call sadd(s,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine smult2(s,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension s(0:1),sres(0:1)
c..computes bounds on the product of two scalars
c..extract lower and upper bounds
      r1l=s(0)
      r1u=s(1)
      r2l=sres(0)
      r2u=sres(1)
c..treat special cases specially
      if(r1u.eq.rzero.and.r1l.eq.r1u.or.
     *	 r2u.eq.rzero.and.r2l.eq.r2u)then
	call siconst(0,sres)
      elseif(r1l.ge.rzero)then
	if(r2l.ge.rzero)then
	  call sscalar(rdprod(r1l,r2l),ruprod(r1u,r2u),sres)
	elseif(r2u.lt.rzero)then
	  call sscalar(rdprod(r1u,r2l),ruprod(r1l,r2u),sres)
	else
	  call sscalar(rdprod(r1u,r2l),ruprod(r1u,r2u),sres)
	endif
      elseif(r1u.lt.rzero)then
	if(r2l.ge.rzero)then
	  call sscalar(rdprod(r1l,r2u),ruprod(r1u,r2l),sres)
	elseif(r2u.lt.rzero)then
	  call sscalar(rdprod(r1u,r2u),ruprod(r1l,r2l),sres)
	else
	  call sscalar(rdprod(r1l,r2u),ruprod(r1l,r2l),sres)
	endif
      else
	if(r2l.eq.rzero)then
	  call sscalar(rdprod(r1l,r2u),ruprod(r1u,r2u),sres)
	elseif(r2u.lt.rzero)then
	  call sscalar(rdprod(r1u,r2l),ruprod(r1l,r2l),sres)
	else
	  rlu=rdprod(r1l,r2u)
	  rul=rdprod(r1u,r2l)
	  ruu=ruprod(r1u,r2u)
	  rll=ruprod(r1l,r2l)
	  call sscalar(rmin2(rlu,rul),rmax2(ruu,rll),sres)
	endif
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine smult3(s1,s2,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension s1(0:1),s2(0:1),sres(0:1)
c..replaces a scalar by its product with an other one.
      call sequal(s2,sres)
      call smult2(s1,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine smadd(s1,s2,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      dimension s1(0:1),s2(0:1),sres(0:1)
      dimension s(0:1)
c..computes bounds for the product of two scalars and
c  adds it to the result.
      call smult3(s1,s2,s)
      call sadd(s,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine squot(s1,s2,sres)
      implicit double precision(a-h,o-z)
      dimension s1(0:1),s2(0:1),sres(0:1)
c..computes bounds for the quotient of the scalers s1 and s2.
      call sinv2(s2,sres)
      call smult2(s1,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine sfak(k,sres)
      implicit double precision(a-h,o-z)
      dimension sres(0:1)
      dimension s(0:1),s1(0:1)
c..computes k faktorial
      call siconst(1,s)
      do 1 i=1,k
      call siconst(i,s1)
      call smult2(s1,s)
    1 continue
      call sequal(s,sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine spower(s,k,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/sconst/szero(0:1),sone(0:1),stwo(0:1),
     *		    shalf(0:1),seight(0:1),semp25(0:1),s4over(0:1)
      common/output/iout
      dimension s(0:1),sres(0:1)
      dimension s1(0:1),st1(0:1)
      dimension s2(0:1)
c..computes the k'th power of s.
c..k negative
      if(k.lt.0)then
        write(iout,*)'error in spower, k not positive'
c..s=0
      elseif(rmaxabs(s).eq.rzero)then
	if(k.eq.0)then
	  call siconst(1,sres)
	else
	  call siconst(0,sres)
	endif
      else
c..k>0, s<>0
	m=k
	call sequal(s,s1)
	call siconst(1,st1)
    1	if(m.gt.0)then
	  mh=m/2
	  if(mh+mh.ne.m)then
	    call smult2(s1,st1)
	    m=m-1
	  else
	    call smult3(s1,s1,s2)
	    call sequal(s2,s1)
	    m=mh
	  endif
	  goto 1
	endif
	call sequal(st1,sres)
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine ssqrt(s,sres)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
      common/output/iout
      dimension s(0:1),sres(0:1)
c..computes bound on the square root of s.
      if(s(0).ge.rzero)then
c..lower bound
	rl=rmax2(rzero,dsqrt(s(0)))
    1	if(ruprod(rl,rl).gt.s(0))then
	  rl=rdown(rl)
	  goto 1
	endif
c..upper bound
	ru=rmax2(rzero,dsqrt(s(1)))
    2	if(rdprod(ru,ru).lt.s(1))then
	  ru=rup(ru)
	  goto 2
	endif
	call sscalar(rl,ru,sres)
c..argument negative
      else
        write(iout,*)'error in ssqrt, argument negative',s(0),s(1)
	call siconst(0,sres)
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine slog(s,sres)
      implicit double precision(a-h,o-z)
      dimension s(0:1),sres(0:1)
c..computes bound on logarithm of s.
      call sscalar(rdown(dlog(s(0))),rup(dlog(s(1))),sres)
      return
      end
c-----------------------------------------------------------------------
      subroutine sexp(s,sres)
      implicit double precision(a-h,o-z)
      dimension s(0:1),sres(0:1)
c..computes bound on the exponential function of s.
      call sscalar(rdown(dexp(s(0))),rup(dexp(s(1))),sres)
      return
      end
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
      subroutine pcheck(p)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension p(0:1,0:nt3,0:nu3)
      print*,'mu=0.',p(1,nt+1,0),p(1,nt+2,0)
      print*,'mu=gn',p(1,nt+1,nu+1),p(1,nt+2,nu+1)
      print*,'mu=ho',p(1,nt+1,nu+2),p(1,nt+2,nu+2)
      return
      end
c-----------------------------------------------------------------------
      subroutine tcheck(t)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension t(0:1,0:nt3)
      print*,t(0,nt+1)
      return
      end
c-----------------------------------------------------------------------
      subroutine ucheck(u)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      dimension u(0:1,0:nu3)
      print*,u(0,0),u(1,0)
      print*,u(0,nu+1)
      return
      end
c-----------------------------------------------------------------------
c
c     5.8. OPERATIONS INVOLVING ONLY REPS
c     ===================================
c
c     rup     :  successor in the set of reps
c     rdown   :  predecessor in the set of reps
c     rabs    :  absolute value of a reps
c     rmax2   :  maximum of two reps
c     rmax3   :  maximum of three reps
c     rmin2   :  minimum of two reps
c     rmin3   :  minimum of three reps
c     rusum   :  upwards rounded sum of reps
c     rdsum   :  downwards rounded sum of reps
c     rddiff  :  downwards rounded diff. of reps
c     ruprod  :  upper bound on product of reps
c     rdprod  :  lower bound on product of reps
c     ruquot  :  upper bound on quotient of reps
c     rdquot  :  lower bound on quotient of reps
c     ruinv   :  upper bound on invers of a rep
c     rdinv   :  lower bound on invers of a rep
c-----------------------------------------------------------------------
c
c   Remark 1: The type reps is in 64-bit IEEE format.
c	      The standard IEEE format ( see [24] ) is as follows.
c	      Sign bit, 11-bit binary exponent with bias
c             of 1023, 52-bit mantissa ( with "hidden" most
c	      significant bit=1 ). The mantissa represents a
c	      number < 2 but not < 1. A zero exponent means a
c	      value of zero, and the exponent 11111111111 means
c             "not a number".
c
c   Remark 2: Bytes are addressed in order; the highest addressed byte
c	      is the least significant mantissa byte.
c	      The order in which the bytes are numbered
c	      does (unfortunately) not form part of the standard
c	      and varies form implementation to implementation.
c-----------------------------------------------------------------------
      double precision function rup(r)
      implicit double precision(a-h,o-z)
      parameter(i20=2**20,i30=2**30)
      parameter(i20m30=i20-i30,i20m31=i20m30-i30,i31m20=-i20m31)
      common/rconst/rzero,rone,rtwo
      common/output/iout
      dimension i(0:1)
      equivalence(r1,i(0))
c..computes nearest representable number to the right of r.
      r1=r
      iexp=i(0)
      ilow=i(1)
c..case where r>0
      if(iexp.gt.0)then
	ilow=ilow+1
	if(ilow.eq.0)then
	  iexp=iexp+1
	  if(iexp.eq.i31m20)then
            write(iout,*)'error in rup, exponent overflow',r
	    rup=r
	    return
	  endif
	endif
c..case where r<0
      elseif(r1.lt.rzero)then
	if(ilow.eq.0)then
	  if(iexp.eq.i20m31)then
	    rup=rzero
	    return
	  else
	    ilow=ilow-1
	    iexp=iexp-1
	  endif
	else
	  ilow=ilow-1
	endif
c..case r=0
      else
	iexp=i20
      endif
      i(0)=iexp
      i(1)=ilow
      rup=r1
      return
      end
c-----------------------------------------------------------------------
      double precision function rdown(rl)
      implicit double precision(a-h,o-z)
c..computes nearest representable number to the left of rl.
      rdown=-rup(-rl)
      return
      end
c-----------------------------------------------------------------------
      double precision function rabs(r)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
      common/rconst/rzero,rone,rtwo
c..computes the absolute value of a number
      r1=r
      if(r1.lt.rzero)r1=-r1
      rabs=r1
      return
      end
c-----------------------------------------------------------------------
      double precision function rmax2(r1,r2)
      implicit double precision(a-h,o-z)
c..computes the maximum of two numbers
      r=r1
      if(r2.gt.r1)r=r2
      rmax2=r
      return
      end
c-----------------------------------------------------------------------
      double precision function rmin2(r1,r2)
      implicit double precision(a-h,o-z)
computes the minimum of two numbers
      r=r1
      if(r2.lt.r1)r=r2
      rmin2=r
      return
      end
c-----------------------------------------------------------------------
      double precision function rmax3(r1,r2,r3)
      implicit double precision(a-h,o-z)
c..computes the maximum  of three numbers
      r=rmax2(r1,r3)
      if(r1.lt.r2)r=rmax2(r2,r3)
      rmax3=r
      return
      end
c-----------------------------------------------------------------------
      double precision function rmin3(r1,r2,r3)
      implicit double precision(a-h,o-z)
c..computes the minimum of three numbers
      r=rmin2(r1,r3)
      if(r1.gt.r2)r=rmin2(r2,r3)
      rmin3=r
      return
      end
c-----------------------------------------------------------------------
      double precision function rusum(ru1,ru2)
      implicit double precision(a-h,o-z)
c..computes an upper bound on the sum ru1+ru2
      rusum=rup(ru1+ru2)
      return
      end
c-----------------------------------------------------------------------
      double precision function rdsum(rl1,rl2)
      implicit double precision(a-h,o-z)
c..computes a lower bound on the sum rl1+rl2
      rdsum=rdown(rl1+rl2)
      return
      end
c-----------------------------------------------------------------------
      double precision function rddiff(r1,r2)
      implicit double precision(a-h,o-z)
c..computes a lower bound on the difference r1-r2.
      rddiff=rdown(r1-r2)
      return
      end
c-----------------------------------------------------------------------
      double precision function ruprod(r1,r2)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
c..computes an upper bound on the product r1*r2.
      if(r1.eq.rzero.or.r2.eq.rone)then
	ruprod=r1
      else
	if(r1.eq.rone.or.r2.eq.rzero)then
	  ruprod=r2
	else
	  ruprod=rup(r1*r2)
	endif
      endif
      return
      end
c-----------------------------------------------------------------------
      double precision function rdprod(r1,r2)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
c..computes a lower bound on the product r1*r2.
      if(r1.eq.rzero.or.r2.eq.rone)then
	rdprod=r1
      else
	if(r1.eq.rone.or.r2.eq.rzero)then
	  rdprod=r2
	else
	  rdprod=rdown(r1*r2)
	endif
      endif
      return
      end
c-----------------------------------------------------------------------
      double precision function ruquot(r1,r2)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
c..computes an upper bound on the quotient r1/r2.
      if(r1.eq.rzero.or.r2.eq.rone)then
	ruquot=r1
      else
	ruquot=rup(r1/r2)
      endif
      return
      end
c-----------------------------------------------------------------------
      double precision function rdquot(r1,r2)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
c..computes a lower bound on the quotient r1/r2.
      if(r1.eq.rzero.or.r2.eq.rone)then
	rdquot=r1
      else
	rdquot=rdown(r1/r2)
      endif
      return
      end
c-----------------------------------------------------------------------
      double precision function ruinv(r)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
c..computes an upper bound on the inverse of r.
      ruinv=rup(rone/r)
      return
      end
c-----------------------------------------------------------------------
      double precision function rdinv(r)
      implicit double precision(a-h,o-z)
      common/rconst/rzero,rone,rtwo
c..computes a lower bound on the inverse of r.
      rdinv=rdown(rone/r)
      return
      end
c-----------------------------------------------------------------------

c     5.8. OPERATIONS ON INDECES
c     ==========================
c
c     itgh     :  gives power of general/higher order of first variable
c     iugh     :  gives power of general/higher order of second variable
c-----------------------------------------------------------------------
      integer function itgh(i)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
c..function returning the order of coefficients,
c  general and higher order bounds for functions of t.
      ii=i
      if(i.eq.nt+1)ii=1
      if(i.eq.nt+2)ii=nt+1
      itgh=ii
      return
      end
c-----------------------------------------------------------------------
      integer function iugh(i)
      implicit double precision(a-h,o-z)
      parameter(nt=75,nt3=nt+3,nu=7,nu3=nu+3)
c..function returning the order of coefficients,
c  general and higher order bounds for functions of u.
      ii=i
      if(i.eq.nu+1)ii=1
      if(i.eq.nu+2)ii=nu+1
      iugh=ii
      return
      end
c-----------------------------------------------------------------------
