c-----------------------------------------------------------
c Chapter 21: Composition of Power Series(p187)
c-----------------------------------------------------------
c   Name of subroutine: POWSER
c
c   Algorithm:Compose power series.
c
c   input: n=10   
c   complier: f77 powser.f
c-----------------------------------------------------------

      parameter(n=10)
      double precision  alpha, a(n),b(n),c(n),d(n)
      integer i,option
      alpha = 7.
      option = 1
      write(*,30)
30    format('Input the coefficient of C(1)...C(N)')
      read(*,40),(c(j),j=1,n)
40    format(10(d5.2))
      call powser(a,b,c,n,alpha,option,d)
      write(*,50)
50    format('The output is ')
      write(*,60),(a(i),i=1,n)
60    format (10(d9.2))
      stop
      end

c-----Subroutine begins here--------------------------------
      subroutine powser(a,b,c,n,alpha,option,d)
      double precision a(n),b(n),c(n),d(n),alpha,alp,r,s,t,v,dfloat
      integer q, option
      ind=1
      q=0
      if(option-3)10,30,40
10    m1=0
      s=1.
      if(option.eq.2) ind=0
      alp=1.
      if(option.eq.1) alp=alpha
      n1=n
15    do 11 j=1,n1
      v=0.
      if(j.eq.1) go to 11
      j1=j-1
      do 12 i=1,j1
12    v=v+a(i)*c(j-i+q)*(alp*(j-i)-ind*i)
11    a(j)=(alp*c(j)+v/dfloat(j))*s
      if(option-2) 43,43,36
30    do 31 i=1,n
31    d(i)=b(1)*c(i)
      do 33 q=1,n
      if(c(q) .ne. 0) go to 34
33    continue
      go to 38
34    s=1./c(q)
      m=1
35    m=m+1
      m1=m*q
      if(m1 .gt. n) go to 38
      if(b(m) .eq. 0.) go to 35
      alp = m
      r = b(m)*c(q) **m
      d(m1) = d(m1) + r
      n1 = n-m1
      if (n1) 38,38,15
36    m1 = m*q 
      do 37 i=1,n1
37    d(i+m1) = d(i+m1)+a(i)*r
      go to 35
38    do 39 i = 1,n
39    a(i) = d(i)
      return
40    t = 1.
      do 41 i = 1,n
      t = t/c(1)
      b(i) = a(i) * t
41    d(i) = c(i) * t
      if (n .eq. 1) return
      do 42 m = 2,n
      s = -d(m)
      m0 = m-1
      do 42 i = m,n
      do 42 l = i,n
      b(l) = b(l) + s*b(l-m0)
42    d(l) = d(l) + s*d(l-m0)
43    return
      end