!*****************************************************
!*     Euclidian division of two polynomials         *
!*             R(x) = P(x) / Q(x)                    *
!* ------------------------------------------------- *
!* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp *
!* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
!* [BIBLI 05].                                       *
!* ------------------------------------------------- *
!* SAMPLE RUN:                                       *
!*                                                   *
!* EUCLIDIAN DIVISION OF TWO POLYNOMIALS:            *
!*                                                   *
!* P(x) = '2x5 +2x3 -x2 +2x -1'                      *
!* Q(x) = 'x3 -1'                                    *
!*                                                   *
!*                                                   *
!* Quotient:                                         *
!*                                                   *
!*+     2X2 +     2                                  *
!*                                                   *
!* Remainder:                                        *
!*                                                   *
!*+ X2 +     2X +     1                              *
!*                                                   *
!* ------------------------------------------------- *
!* Functions used (of module Polynoms):              *
!*                                                   *
!*  AddNumber(), EnterPolynom(), DivNumber(),        *
!*  DisplayPolynom(), MultNumber() and SetNumber().  *
!*                                                   *
!*                       F90 version by J-P Moreau.  *
!*                           (www.jpmoreau.fr)       *
!*****************************************************
PROGRAM DIVPOL
Use Polynoms

type(ar_polynom)  P,Q,H,R
integer DivPolynom

  print *,''
  print *,' EUCLIDIAN DIVISION OF TWO POLYNOMIALS:'
  print *,''
  if (EnterPolynom(' P(X) = ', P).eq.0)  stop
  if (EnterPolynom(' Q(X) = ', Q).eq.0)  stop
  print *,''
  print *,''
  if (DivPolynom(P,Q,H,R).ne.0) then
    print *,' Quotient:'
    call DisplayPolynom(H)
    print *,''
    print *,' Remainder:'
    call DisplayPolynom(R)
  else
    print *,' Error in Euclidian division.'
  end if
  print *,''
  print *,''

END


  integer Function DivPolynom(P,Q,H,R)
    USE POLYNOMS
    type(ar_polynom) P,Q,H,R
    type(ar_number)  u, v
    !The Q polynomial must be <> zero
    if (Q%degree.eq.0.and.(Q%coeff(0)%value.eq.0.d0.or.Q%coeff(0)%p.eq.0)) then
      DivPolynom=0
	  return
	end if
    R=P; H%degree=P%degree - Q%degree
    if (H%degree<0) then
      H%degree=0; 
	  if (SetNumber(H%coeff(0),'0').eq.0) then
        DivPolynom=0
	    return
	  end if
    else
      do i=H%degree, 0, -1
        if (DivNumber(R%coeff(R%degree),Q%coeff(Q%degree),H%coeff(i)).eq.0) then
          DivPolynom=0
	      return
		end if
        do j=i, R%degree
          if (MultNumber(H%coeff(i),Q%coeff(j-i), u).eq.0) then
            DivPolynom=0
	        return
		  end if
          u%p=-u%p; u%value=-u%value
          if (AddNumber(R%coeff(j),u,v).eq.0) then
            DivPolynom=0
	        return
		  end if
		  R%coeff(j)=v
        end do
        if (R%degree > 0)  R%degree=R%degree-1
      end do
      do while (dabs(R%coeff(R%degree)%value) < AR_SMALL.and.R%degree>0)
	    R%degree=R%degree-1
      end do
    end if
    DivPolynom=1
  End !of function DivPolynom


!end of file divpol.f90