!***************************************************** !* Division of two polynomials by increasing powers * !* k * !* P(x) = Q(x) * H(x) + x * R(x) * !* ------------------------------------------------- * !* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp * !* and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * !* [BIBLI 05]. * !* ------------------------------------------------- * !* SAMPLE RUN: * !* * !* DIVISION OF TWO POLYNOMIALS BY INCREASING POWERS: * !* * !* P(x) = '2x2 +x' * !* Q(x) = 'x3 -x2 +1' * !* * !* Order: 4 * !* * !* * !* Quotient: * !* * !*+ X3 + 2X2 + X * !* * !* Remainder: * !* * !*- X2 - X + 1 * !* * !* ------------------------------------------------- * !* Functions used (of module Polynoms): * !* * !* AddNumber(), EnterPolynom(), DivNumber(), * !* DisplayPolynom(), MultNumber() and SetNumber(). * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !***************************************************** !Explanations: !Given two polynomials P(x) and Q(x) and an integer k, !it exists a couple of polynomials H(x) and R(x), such !as: Px) = Q(x).H(x) + x^k.R(x) !The degree of H(x) is < k. !Note: the constant coefficient of Q(x) must be <> 0. !----------------------------------------------------- PROGRAM DIVPOL1 Use Polynoms type(ar_polynom) P,Q,H,R integer DivPolynom1 print *,'' print *,' DIVISION OF TWO POLYNOMIALS BY INCREASING POWERS:' print *,'' if (EnterPolynom(' P(X) = ', P).eq.0) stop if (EnterPolynom(' Q(X) = ', Q).eq.0) stop print *,'' write(*,10,advance='no'); read *, k print *,'' print *,'' if (DivPolynom1(P,Q,k,H,R).ne.0) then print *,' Quotient:' call DisplayPolynom(H) print *,'' print *,' Remainder:' call DisplayPolynom(R) else print *,' Error in division by increasing powers.' end if print *,'' print *,'' 10 format(' Order: ') END !Division of two polynomials by increasing powers integer Function DivPolynom1(P,Q,k,H,R) USE POLYNOMS type(ar_polynom) P,Q,H,R type(ar_number) u !The Q polynomial constant term and k must be <> zero if (k<=0.or.Q%coeff(0)%value.eq.0.d0) then DivPolynom1=0 return end if mx=Q%degree+k-1 if (mxQ%degree) then do i=Q%degree+1, mx if (SetNumber(Q%coeff(i),'0').eq.0) then DivPolynom1=0 return end if end do end if !set polynomial R to zero R%degree=0 do i=0, AR_MAXPOL R%coeff(i)%value=0.d0 R%coeff(i)%p=0; R%coeff(i)%q=1 end do do i=0, P%degree R%coeff(i)=P%coeff(i) end do if (k>Q%degree) then do i=Q%degree+1, k if (SetNumber(Q%coeff(i),'0').eq.0) then DivPolynom1=0 return end if end do end if do i=0, k-1 if (DivNumber(R%coeff(i),Q%coeff(0),H%coeff(i)).eq.0) then DivPolynom1=0 return end if do j=1, mx if (MultNumber(H%coeff(i),Q%coeff(j-i), u).eq.0) then DivPolynom1=0 return end if u%p=-u%p; u%value=-u%value if (AddNumber(R%coeff(j), u, R%coeff(j)).eq.0) then DivPolynom1=0 return end if end do end do H%degree=k-1 do while (H%degree>0.and.H%coeff(H%degree)%value.eq.0.d0) H%degree = H%degree - 1 end do R%degree=mx-k do i=0, R%degree R%coeff(i)=R%coeff(i+k) end do do while (R%degree>0.and.R%coeff(R%degree)%value.eq.0.d0) R%degree = R%degree - 1 end do DivPolynom1=1 End !of function DivPolynom1 !end of file divpol1.f90