!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM Polynome
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Dieses Beispiel umfasst einige Unterprogramme zur
!  Manipulation von Polynomen. Jedes Polynom wird
!  dabei durch seine Koeffizienten p_n (n=0 bis npol)
!  dargestellt. Einige Anwendungen werden demonstriert.
!  Die hoechste Potenz wird durch den Wert von npol bestimmt.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

MODULE size
   integer, parameter::npol = 10
END MODULE size

PROGRAM Polynome

   use size
   implicit none


   real, dimension(0:npol)::p1,p2,p3
   real:: fac

   write(*,*) 'Zwei Arten, sin(x) zu bestimmen:'

   write(*,*) '1. Direkt, mit Psin(p)'
   call Psin(p1)
   call Pprint(p1)

   write(*,*) '2. Mit Substitution, P(x)=x, sin(P(x)).'
   p2=0
   p2(1) = 1
   call PvonP(p1, p2, p3)
   call Pprint(p3)

   write(*,*) 'Nun berechne ich sin(x-x^2+2x^3)'
!  Definition der Koeffizienten von p1: 0 + x -x^2 +2x^3
   p1=0
   p1(1) = 1
   p1(2) = -1
   p1(3) = 2
   call Psin(p2)
!  p2 ist noch immer =x
!  Berechne sin(p1):
   call PvonP(p2, p1, p3)
   call Pprint(p3)

   write(*,*) 'Nun berechne ich exp(1+2x):'
   p1=0
   p1(0) =1
   p1(1)=2
!  exp (x) -->p2 :
   call Pexp(p2)
!  p2(p1) --> p3 :}
!  Beachte: Der konstante Term p1(0) muss als e^p1(0)
!  gesondert behandelt werden. Es ist
!   exp(p1) = exp(p1(0)) * exp(p1-p1(0))
   fac = exp(p1(0))
   p1(0) = 0
   call PvonP(p2, p1, p3)
   call Smul(fac, p3, p1)
   call Pprint(p1)

   write(*,*) ' Versuchen Sie nun selbst, eine Polynomdivision mit'
   write(*,*) ' Hilfe von pmul und pinv durchzufuehren!'
   
END PROGRAM Polynome


SUBROUTINE Padd (a,p1,p2,p3)

   use size
   implicit none

   real, intent(in) ::a
   real, dimension(0:npol),intent(in)::p1,p2
   real, dimension(0:npol),intent(inout)::p3

   p3=p1+a*p2

END SUBROUTINE Padd 

SUBROUTINE Pmul (p1, p2, p3)
!  p1 * p2 --> p3

   use size
   implicit none

   real, dimension(0:npol),intent(in)::p1,p2
   real, dimension(0:npol),intent(inout)::p3
   integer::n1,n2,n3
 
   do n3 = 0,npol
     p3(n3) = 0
     do n2 = 0,n3
       n1 = n3 - n2;
       p3(n3) = p3(n3) + p1(n1) * p2(n2);
     end do
   end do
END SUBROUTINE Pmul 

SUBROUTINE Smul (a, p1,p2)

   use size
   implicit none

   real, intent(in) ::a
   real, dimension(0:npol),intent(in)::p1
   real, dimension(0:npol),intent(inout)::p2

   p2=a*p1
  
END SUBROUTINE Smul

SUBROUTINE Pprint (p)
!  Print p

   use size
   implicit none

   real, dimension(0:npol),intent(in)::p
   integer::i

   do i = 0,npol
     write(*,*) 'c(',i,') = ',p(i)
   end do
   write(*,*) 'Weiter mit CR...'
   read(*,*)
    
END SUBROUTINE Pprint

SUBROUTINE Pexp (p)
!  eps(x) -> p

   use size
   implicit none

   real, dimension(0:npol),intent(out)::p
   real::c
   integer:: n

   c=1
   p(0)=1
  
   do n=1,npol
     c=c/n
     p(n)=c
   end do

END SUBROUTINE Pexp

SUBROUTINE Pinv (p)
! 1/(1+x) --> p

   use size
   implicit none

   real, dimension(0:npol),intent(out)::p
   real::c
   integer:: n

   c = 1
  
   do n=1,npol
     p(n) = c
     c = -c
   end do
    
END SUBROUTINE Pinv

SUBROUTINE Psin (p)
!  sin(x) --> p 

   use size
   implicit none

   real, dimension(0:npol),intent(out)::p
   real::c
   integer:: n,i

   p=0
   c = 1
   if(npol > 0) then
     p(1)=c
     if(npol >2) then
       do i=1,int((npol-1)/2)
          n = 2*i + 1
          c = -c / (n*(n-1))
          p(n) = c
       end do
     end if
   end if
   
END SUBROUTINE Psin

SUBROUTINE PvonP (p1, p2, p3)
!  Da ein konstantes Glied in p2  (also p2(0)<>0)
!  zu jeder Ordnung in der Reihe p3 beitragen kann,
!  und dies vor allem bei 'unendlichen' Reihen p1 zu
!  Eindeutigkeitsproblemen fuehrt, wird nur
!  p1 ( p2 - p2(0) ) --> p3
!  berechnet. Auf entsprechende Verwendung ist dabei
!  zu achten!   

   use size
   implicit none

   real, dimension(0:npol),intent(in)::p1,p2
   real, dimension(0:npol),intent(inout)::p3

   real, dimension(0:npol)::pp, pa, pb
   integer::i
 
   pp=p2
   pa=0
   p3=0

!  p2-p2(0)--> pp 

   pp(0) = 0
   pa(0) = 1
   p3(0) = p1(0)
  
   do i = 1,npol

!    pp(=p2^(n-1))*p2 --> pb (=p2^n)}
     call Pmul(pa, pp, pb)

!    p1(n)*pb(=pp^n) wird zu p3 addiert}
     call Padd(p1(i),p3, pb, p3)

     pa = pb

   end do
   
END SUBROUTINE PvonP
