!***************************************************** !* GCD and SCM of two polynomials * !* * !* ------------------------------------------------- * !* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp * !* and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * !* [BIBLI 05]. * !* ------------------------------------------------- * !* SAMPLE RUN: * !* * !* GCD AND SCM OF TWO POLYNOMIALS: * !* * !* P(x) = 'x5 +5x4 +7x3 +5x2 +x -1' * !* Q(x) = 'x4 +4x3 -7x +2' * !* * !* * !* Greatest common divisor: * !* * !*+ X2 + 3 X - 1 * !* * !* * !* Smallest common multiple: * !* * !*+ X7 + 6 X6 + 10 X5 + 2 X4 - 8 X3 * !*- 10 X2 - 3 X + 2 * !* * !* ------------------------------------------------- * !* Functions used (of module Polynoms): * !* * !* AddNumber(), EnterPolynom(), DivNumber(), * !* DisplayPolynom(), MultNumber() and SetNumber(). * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !***************************************************** !Explanations: !The Greatest common divisor (GCD) and the Smallest Common !multiple (SCM) of two polynomials are defined but for a !coefficient. That is to say, if a polynomial R(x) is the !GCD of P(x) and Q(x), any polynomial µ*R(x) is also a GCD. !To find the GCD, after permutating P(x) and Q(x), if the !degree of P is lower than the degree of Q, we first perform !the Euclidian division P/Q: P=Q.H1+R1 (see program Divpol.f90). !if R1<>0, we then divide Q by R1: Q=R1.H2+R2. !if R2<>0, we continue until Rn=0. Rn-1 is then the GCD of !P(x) and Q(x). !To find the SCM of P(x) and Q(x), we use the relation: ! P(x).Q(x)=GCD(P,Q).SCM(P,Q). !--------------------------------------------------------------- PROGRAM GCDPOL Use Polynoms type(ar_polynom) P,Q,R, P1,Q1 integer GCDPolynom, SCMPolynom print *,'' print *,' GCD AND SCM OF TWO POLYNOMIALS:' print *,'' if (EnterPolynom(' P(X) = ', P).eq.0) stop if (EnterPolynom(' Q(X) = ', Q).eq.0) stop print *,'' print *,'' P1=P; Q1=Q if (GCDPolynom(P1,Q1,R).ne.0) then print *,' Greatest common divisor:' call DisplayPolynom(R) else print *,' Error in looking for GCD.' end if print *,'' print *,'' if (SCMPolynom(P,Q,R).ne.0) then print *,' Smallest common multiple:' call DisplayPolynom(R) else print *,' Error in looking for SCM.' 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 MultPolynom 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 ! GCD of two polynomials integer Function GCDPolynom(P,Q,R) USE POLYNOMS type(ar_polynom) P,Q,R, V type(ar_number) z integer*4 pg, pp, rr integer bb, DivPolynom if (P%degree.eq.0.and.P%coeff(P%degree)%value.eq.0) then GCDPolynom=0 return end if if (Q%degree.eq.0.and.Q%coeff(Q%degree)%value.eq.0) then GCDPolynom=0 return end if if (Q%degree>P%degree) then R=P; P=Q; Q=R end if R=Q do while (R%degree>0) if (DivPolynom(P,Q,V,R).eq.0) then GCDPolynom=0 return end if if (R%coeff(R%degree)%is_real.ne.0) then do while (R%degree>0.and.dabs(R%coeff(R%degree)%value)0.and.R%coeff(R%degree)%p.eq.0) R%degree=R%degree-1 end do endif if (R%degree.eq.0) then if (R%coeff(0)%is_real.ne.0) then if (dabs(R%coeff(0)%value)>AR_SMALL) Q%degree=0 else if (R%coeff(0)%p.ne.0) Q%degree=0 end if else P=Q; Q=R end if do i=0, Q%degree-1 if (DivNumber(Q%coeff(i),Q%coeff(Q%degree),Q%coeff(i)).eq.0) then GCDPolynom=0 return end if end do if (SetNumber(Q%coeff(Q%degree), '1').eq.0) then GCDPolynom=0 return end if end do R=Q if (Q%degree.eq.0) then bb=SetNumber(R%coeff(0), '1') else pg=0 do i=0, R%degree if (R%coeff(i)%is_real.eq.0) then if (pg.eq.0) then pg=R%coeff(i)%p else rr=GCD(pg,R%coeff(i)%p); pg=rr end if end if end do pp=0 do i=0, R%degree if (R%coeff(i)%is_real.eq.0) then if (pp.eq.0) then pp=R%coeff(i)%q else rr=GCD(pp,R%coeff(i)%q) pp=pp*R%coeff(i)%q/rr end if end if end do if (pg.ne.0) then z%is_real=0; z%p=pp z%q=pg; z%value=pp/pg do i=0, R%degree if (MultNumber(R%coeff(i), z, R%coeff(i)).eq.0) then GCDPolynom=0 return end if end do end if end if GCDPolynom=1 End !of function GCDPolynom ! SCM of two polynomials integer Function SCMPolynom(P,Q,R) USE POLYNOMS type(ar_polynom) P,Q,R, U,V,V1 integer GCDPolynom, DivPolynom V=P; V1=Q !to preserve P and Q after GCDPOlynom if (GCDPolynom(V,V1,U).eq.0) then SCMPolynom=0 return end if if (MultPolynom(P,Q,V).eq.0) then SCMPolynom=0 return end if if (DivPolynom(V,U,R,V).eq.0) then SCMPolynom=0 return end if SCMPolynom=1 End ! end of file gcdpol.f90