c.....given a funtion G(X)=b_1 +.b_2*X + ... + b_n* X^{n-1}
c.....we get g(x)=a_1+ a_2*x^2 +...+a_n*x^{2*(n-1)}
c.....where g(x)=G((x^2-c)/r)      
      subroutine rnscale(n,a,b,tran,dilat)
      integer n
      double precision a(n), b(n), tran, dilat
      double precision temp(n),aux
      call matmov(n,a,temp)
      aux=1.d0/dilat
      call matscale(n,temp,aux)
      aux=-tran
      call shift(n,temp,b,aux)
      return
      end
      
c...... given a funtion g(x)=a_1 + a_2 *x^2 +...+ a_n* x^{2*(n-1)}
c...... we get G((x^2-c)/r)=g(x)
      subroutine rgscale(n,a,b,tran,dilat)
      integer n
      double precision a(n), b(n), tran, dilat
      double precision aux
      aux=tran/dilat
      call matscale(n,a,dilat)
      call shift(n,a,b,aux)
      return
      end

      subroutine normerr(n,a ,tau)
      integer n
      double precision a(n), tau
      integer i
c     double precision aux
      tau = 0.d0
      do 10 i=1,n
         tau=tau + dabs(a(i))
10    continue
      return
      end

      
      subroutine matzer (n,a)            
      integer n 
      double precision a(n)
c     
c     z(t)=0                                     
c     
      integer i
      do 10 i=1,n 
         a(i)=0.0d0
10    continue
      return                                                           
      end                                                              

      subroutine matmov(n,a,b)                                         
      integer n
      double precision a(n),b(n)                                       
      integer i
      do 10 i=1,n                                                      
         b(i) = a(i)
10    continue
      return                                                           
      end                                                              

      subroutine matsca(n,a,z)                                         
      integer n 
      double precision a(n),z                                          
      integer i
c     
c     z*p(t)  where p(t)=a_1+   + a_n t^(n-1)
c     
      do 10 i=1,n                                                      
         a(i)=a(i)*z
10    continue
      return
      end                                                              

      subroutine matadd(n,a,b,c)                                       
      integer n 
      double precision a(n),b(n),c(n)                                  
      integer i
      do 10 i=1,n                                                      
         c(i)=a(i)+b(i) 
10    continue
      return
      end

      subroutine matsub(n,a,b,c)                                       
      integer n 
      double precision a(n),b(n),c(n)                                  
      integer i
      do 10 i=1,n                                                      
         c(i)=a(i)-b(i)
10    continue
      return
      end

      subroutine shift(n,a,b,c)                                       
      integer n 
      double precision a(n),b(n),c                                     
      integer j,k     
c     given polinomial p(t), we get polinomial p(t+c)
      call matmov(n,a,b)                                              
      do 11 j=1,n-1 
      do 12 k=n-1,j,-1                                                    
            b(k)=b(k)+c*b(k+1)                                               
12      continue                                                         
11    continue
      return
      end


      subroutine matscale(n, a, kappa)
      integer n 
      double precision a(n),kappa                                      
      integer i
      double precision w
c     q(t)=p(kappa* t)
      w=1.0d0                                                          
      do 10 i=2,n                                                      
         w=kappa*w                                                        
         a(i)=w*a(i)                                                      
10    continue
      return
      end                                                              

      subroutine prod(n,a,b,c)                                         
      integer n 
      double precision a(n),b(n),c(n)                                  
      integer i,ii,j
      double precision w                                               
      do 10 ii=1,n                                                     
         i=n+1-ii
         w=0.0d0
         do 11 j=1,i                                                      
            w=w+a(j)*b(i+1-j)                                                
 11      continue                                                         
         c(i)=w                                                           
 10   continue                                                         
      return
      end

      subroutine deriv(n,a,b)                                          
      integer n 
      double precision a(n),b(n)                                       
      integer i, nm1
      nm1=n-1                                                          
      do 10 i=1,nm1                                                    
         b(i)=dble(float(i))*a(i+1)                                       
10    continue                                                         
      b(n)=0.d0
      return                                                           
      end                                                              

      subroutine compos(n,x,y,z)                                     
      integer n                                                        
      double precision x(n),y(n),z(n)    
c     z(t) =  y( x(t))
      integer i,nm1,j
      call matzer(n,z)                                                
      z(1)=y(n)                                                        
      nm1=n-1                                                          
      do 10 j=1,nm1                                                    
         i=n-j                                                            
         call prod(n,x,z,z)                                               
         z(1) = z(1) + y(i)
 10   continue
      return                                                           
      end                                                              
      
      subroutine evalu(n,a,x,y)                                        
      integer n 
      double precision a(n),x,y                             
      integer k
c     given argument x, we compute y = p(x)
      y=0.0d0
      do 10 k=1,n                                                      
         y=(x*y)+a(n-k+1)                                                 
 10   continue
      return
      end   

      subroutine matprint(n, a)
      integer n
      double precision a(n)
      integer i
 
      do 10 i = 1,n 
         write(*,*) a(i)
 10   continue
      return
      end

      subroutine matvec(n, a, v, av)
      integer n
      double precision a(n,n), v(n), av(n)
      integer i,j
      double precision w
      do 10 i = 1, n
         w = 0.0d0 
         do 20 j = 1, n
            w = w + a(i,j)*v(j)
 20      continue
         av(i) = w
 10   continue
      return
      end
          
      subroutine matrmov(n,a,k,opmat)
      integer n,k
      integer i
      double precision a(n), opmat(n,n)
      do 10 i=1,n
         opmat(i,k)=a(i)
 10      continue
         return
         end


c   _______________________________________________________________
c     given pol. p(t)= a_1 + ... a_n*t^{n-1} and constants
c      constants r and h we get
c     q(t)= p(r* t + h)= b_1 +...+b_n*t^{n-1}      
      subroutine mataffine(n,a,b,h,r)
      integer n
      double precision a(n), h, r
      double precision b(n)

      call shift(n,a,b,h)
      call matscale(n,b,r)
      return
      end

c ____________________________________________________________
c     given a polynomial p(t)=a_1+...+a_n*t^{n-1} and a linear 
c     polinomial l(t)=b + m*t, this subroutine computes 
c     the product p(t)*(b+ m*t). p(t) enters as input and is
c     q(t)=c_1 +... c_n*t^{n-1} returns as output
      subroutine prodaffine(n,a,c,b,m)
      integer n
      integer i
      double precision a(n), b, m
      double precision c(n)
      double precision atemp(n)
      atemp(1)=b*a(1)
      do 10 i=2,n
         atemp(i)=b*a(i) +m*a(i-1)     
 10   continue
      call matmov(n,atemp,c)
      return
      end

c ------------------------------------------------------------------
c  given a polinomial f and integer m, the following soubroutine
c  computes f^m=f*...*f m times

      subroutine fm(n,m,f,ftom)
      integer n, m
      double precision f(n)
      double precision newf(n), newftemp(n)
      double precision ftom(n)
      integer mtemp, mres
c     initial values: f^m(x)=1 and newf(x)=f(x)
      call matzer(n,ftom)
      ftom(1)=1.d0
      call matmov(n,f,newf)
c     computes f^m(x)
      mtemp = m
 1    mres= mod(mtemp,2)
      if(mres.eq.1) then
         call prod(n,ftom,newf,newftemp)
         call matmov(n,newftemp,ftom)
         if(mtemp.eq.1) goto 2
         mtemp=mtemp-1
      else
         mtemp=mtemp/2
         call prod(n,newf,newf,newftemp)
         call matmov(n,newftemp,newf)
      endif
      goto 1
 2    return
      end