!*****************************************************
!*  This program calculates R(x) = P(Q(x)),          *
!*  P(x), Q(x) and R(x) being polynomials.           *
!* ------------------------------------------------- *
!* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp *
!* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
!* [BIBLI 05].                                       *
!* ------------------------------------------------- *
!* SAMPLE RUN:                                       *
!*                                                   *
!* SUBSTITUTION OF TWO POLYNOMIALS:                  *
!*                                                   *
!* P(x) = '-5x2 +3'                                  *
!* Q(x) = '2x3 -x +5'                                *
!*                                                   *
!*-    20 X6 +    20 X4 -    100 X3 -     5 X2       *
!*+    50 X -   122                                  *
!*                                                   *
!*                                                   *
!* SUBSTITUTION OF TWO POLYNOMIALS:                  *
!*                                                   *
!* P(x) = x4                                         *
!* Q(x) = '2x2 -3'                                   *
!*                                                   *
!*+    16 X8 -    96 X6 +    216 X4 -   216 X2       *
!*+    81                                            *
!*                                                   *
!* ------------------------------------------------- *
!* Functions used (of module Polynoms):              *
!*                                                   *
!*  AddNumber(), EnterPolynom(), DisplayPolynom(),   *
!*  MultNumber() and SetNumber().                    *
!*                                                   *
!*                       F90 version By J-P Moreau.  *
!*                           (www.jpmoreau.fr)       *
!*****************************************************
PROGRAM SUBSTPOL
Use Polynoms

type(ar_polynom) P,Q,R
integer res, SubstPolynom

  print *,''
  print *,' SUBSTITUTION OF TWO POLYNOMIALS:'
  print *,''
  if (EnterPolynom('  P(X) = ', P).eq.0)  stop ' Error in Enter Polynomial.'
  if (EnterPolynom('  Q(X) = ', Q).eq.0)  stop ' Error in Enter Polynomial.'
  print *,''
  res=SubstPolynom(P,Q,R)
  if (res.ne.0) then
    call DisplayPolynom(R)
  else
    print *,' Error in polynomial substitution.'
  end if
  print *,''
  print *,''
END

  ! P(X) * Q(X) = R(X)
  integer Function MultPolynom(P,Q,R1)
  USE POLYNOMS
    type(ar_polynom) P,Q,R,R1
    type(ar_number)  u
    R%degree=0                !set R polynomial to zero
	do i=0, AR_MAXPOL
      R%coeff(i)%value=0
	  R%coeff(i)%p=0
      R%coeff(i)%q=1
	end do
    !verify that P and Q are not void
    if (P%degree.eq.0.and.P%coeff(0)%value.eq.0) then 
	  MultPolynom=0
	  return
    end if
    if (Q%degree.eq.0.and.Q%coeff(0)%value.eq.0) then 
	  MultPolynom=0
	  return
    end if
    R%degree=P%degree+Q%degree
	
    if (R%degree>AR_MAXPOL) then     !R degree is too big
	  MultPolynom=0
	  return
	end if
    do n=0, R%degree
      if (SetNumber(R%coeff(n),'0').eq.0) then
	    MultPolynom=0
	    return
	  end if
      do i=0, P%degree
        j=n-i
        if (j>=0.and.j<=Q%degree) then
          if (MultNumber(P%coeff(i),Q%coeff(j),u).eq.0) then
	        MultPolynom=0
	        return
		  end if
          if (AddNumber(R%coeff(n),u,R%coeff(n)).eq.0) then
	        MultPolynom=0
	        return
		  end if
        end if
      end do
    end do
    MultPolynom=1
	R1=R
  End !of function MultPol
 
  ! a P(X) + b Q(X) = R(X)
  integer Function CombiPolynom(P,Q,a,b,R1)
    USE POLYNOMS
    type(ar_polynom) P,Q,R1, R 
	type(ar_number)  a, b, u, v
	integer degree
    if (Q%degree > P%degree) then 
	  degree=Q%degree
    else 
	  degree=P%degree
    end if 
    if (degree > AR_MAXPOL) then       !degree of R too big
	  CombiPolynom=0
	  return
    end if
    do i=0, degree
      if (SetNumber(u,'0').eq.0) then
	    CombiPolynom=0
	    return
      end if
      if (SetNumber(v,'0').eq.0) then
	    CombiPolynom=0
	    return
      end if
      if (i<=P%degree) then
        if (MultNumber(a,P%coeff(i), u).eq.0) then
	      CombiPolynom=0
	      return
        end if
      end if
      if (i<=Q%degree) then
        if (MultNumber(b,Q%coeff(i), v).eq.0) then
	      CombiPolynom=0
	      return
        end if
      end if 
      if (AddNumber(u,v,R%coeff(i)).eq.0) then
	    CombiPolynom=0
	    return
      end if
    end do
    do while (degree>0.and.dabs(R%coeff(degree)%value)<AR_SMALL)
      degree=degree-1
    end do
	R1%degree=degree
	do i=0, degree
	  R1%coeff(i)=R%coeff(i)
    end do
    CombiPolynom=1
  End !of function CombiPolynom

  ! R(X) = P(Q(X))
  integer Function SubstPolynom(P,Q,R)
    USE POLYNOMS 
	type(ar_polynom) P,Q,R, V
    type(ar_number)  a
	integer CombiPolynom
    !set R polynomial to zero
    R%degree=0
	do i=0, AR_MAXPOL
      R%coeff(i)%value=0
	  R%coeff(i)%p=0
      R%coeff(i)%q=1
	end do
    R%degree=P%degree * Q%degree
    if (R%degree>AR_MAXPOL) then           !R degree is too big
	  SubstPolynom=0
	  return
    end if
    R%coeff(0)=P%coeff(0); V%degree=0
    if (SetNumber(V%coeff(0),'1').eq.0) then
	  SubstPolynom=0
	  return
	end if
    if (SetNumber(a,'1').eq.0) then
	  SubstPolynom=0
	  return
	end if
    do i=1, P%degree
      if (MultPolynom(Q,V,V).eq.0) then
	    SubstPolynom=0
	    return
	  end if
      if (CombiPolynom(R,V,a,P%coeff(i),R).eq.0) then
	    SubstPolynom=0
	    return
	  end if
    end do
    SubstPolynom=1
  End !of function SubsPolynom


!end of file substpol.f90