'***************************************************** '* GCD and SCM of two polynomials * '* * '* ------------------------------------------------- * '* Ref.: "Mathematiques 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: * '* * '* 2 * '* X + 3X - 1 * '* * '* * '* Smallest common multiple: * '* * '* 7 6 5 4 3 2 * '* X + 6X + 10X + 2X -8X - 10 X -3X +2 * '* * '* * '* BASIC 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.pas). ' 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). '----------------------------------------------------------------- defdbl a-h,o-z defint i-n MAXINT = 32767 'Maximum integer number MAXPOL = 20 'Maximum degree for a polynomial SMALL = 1e-20 'small real number ' number zz (real, integer or fractional) ' isrzz=1: zz is real (double precision) ' isrzz=0: zz is integer or fractional ' zzv: value of zz if real ' ipzz: integer value of zz if integer (or numerator value if fractional) ' iqzz: denominator value if fractional (=1 if integer) 'a polynomial P is defined by: 'integer ipdeg (degree of polynomial) 'number cp(MAXPOL) (coefficients of polynomial, real, integer or fractional) 'predefine 4 polynomials P,Q,R DIM isrcp(MAXPOL),cpv(MAXPOL),ipcp(MAXPOL),iqcp(MAXPOL) DIM isrcq(MAXPOL),cqv(MAXPOL),ipcq(MAXPOL),iqcq(MAXPOL) DIM isrcr(MAXPOL),crv(MAXPOL),ipcr(MAXPOL),iqcr(MAXPOL) 'auxiliary polynomials P1,Q1,U,V DIM isrcp1(MAXPOL),cpv1(MAXPOL),ipcp1(MAXPOL),iqcp1(MAXPOL) DIM isrcq1(MAXPOL),cqv1(MAXPOL),ipcq1(MAXPOL),iqcq1(MAXPOL) DIM isrcu(MAXPOL),cuv(MAXPOL),ipcu(MAXPOL),iqcu(MAXPOL) DIM isrcv(MAXPOL),cvv(MAXPOL),ipcv(MAXPOL),iqcv(MAXPOL) cls print print " GCD AND SCM OF TWO POLYNOMIALS:" print tx$=" P(x) = ": gosub 4000 'Enter P(x) 'save P(x) in P1(x) ipdeg1=ipdeg for i=0 to MAXPOL isrcp1(i)=isrcp(i):cpv1(i)=cpv(i):ipcp1(i)=ipcp(i):iqcp1(i)=iqcp(i) next i tx$=" Q(x) = ": gosub 4000 'Enter P(x) 'save P(x) in Q(x) iqdeg=ipdeg for i=0 to MAXPOL isrcq(i)=isrcp(i):cqv(i)=cpv(i):ipcq(i)=ipcp(i):iqcq(i)=iqcp(i) next i 'save Q(x) in Q1(x) iqdeg1=iqdeg for i=0 to MAXPOL isrcq1(i)=isrcq(i):cqv1(i)=cqv(i):ipcq1(i)=ipcq(i):iqcq1(i)=iqcq(i) next i 'restore P(x) from P1(x) ipdeg=ipdeg1 for i=0 to MAXPOL isrcp(i)=isrcp1(i):cpv(i)=cpv1(i):ipcp(i)=ipcp1(i):iqcp(i)=iqcp1(i) next i print:print gosub 700 'call GCDPolynom(P,Q,R) 'Note: P(x) and Q(x) are modified during the process if IERROR<>0 then print " Error in looking for GCD." else 'copy R(x) in P(x) for printing ipdeg=irdeg for i=0 to ipdeg isrcp(i)=isrcr(i):cpv(i)=crv(i):ipcp(i)=ipcr(i):iqcp(i)=iqcr(i) next i print " Greatest common divisor:" gosub 5000 'Display P(x) end if 'restore P(x) from P1(x) ipdeg=ipdeg1 for i=0 to MAXPOL isrcp(i)=isrcp1(i):cpv(i)=cpv1(i):ipcp(i)=ipcp1(i):iqcp(i)=iqcp1(i) next i 'restore Q(x) from Q1(x) iqdeg=iqdeg1 for i=0 to MAXPOL isrcq(i)=isrcq1(i):cqv(i)=cqv1(i):ipcq(i)=ipcq1(i):iqcq(i)=iqcq1(i) next i print gosub 750 'call SCMPolynom(P,Q,R) if IERROR<>0 then print " Error in looking for SCM." else 'copy R(x) in P(x) for printing ipdeg=irdeg for i=0 to ipdeg isrcp(i)=isrcr(i):cpv(i)=crv(i):ipcp(i)=ipcr(i):iqcp(i)=iqcr(i) next i print " Smallest common multiple:" gosub 5000 'Display P(x) end if END 'of main program 'Function MultPolynom P(X) * Q(X) = R(X) 500 IERROR=0 'set R polynomial to zero irdeg=0 for i=0 to MAXPOL isrcr(i)=0 crv(i)=0# ipcr(i)=0 iqcr(i)=1 next i 'verify that P and Q are not void if ipdeg=0 and cpv(0)=0 then IERROR=1 return end if if iqdeg=0 and cqv(0)=0 then IERROR=1 return end if irdeg=ipdeg+iqdeg if irdeg > MAXPOL then 'R degree is too big IERROR=1 return end if for n=0 to irdeg 'set R coeff(n) to zero isrcr(n)=0:crv(n)=0#:ipcr(n)=0:iqcr(n)=1 for i=0 to ipdeg j=n-i if j>=0 and j<=iqdeg then isrxx=isrcp(i):xxv=cpv(i):ipxx=ipcp(i):iqxx=iqcp(i) 'xx=P coeff(i) isryy=isrcq(j):yyv=cqv(j):ipyy=ipcq(j):iqyy=iqcq(j) 'yy=Q coeff(j) gosub 3100 'zz=xx*yy gosub 803 'yy=zz isrxx=isrcr(n):xxv=crv(n):ipxx=ipcr(n):iqxx=iqcr(n) 'xx=R coeff(n) gosub 3000 'zz=xx+yy isrcr(n)=isrzz:crv(n)=zzv:ipcr(n)=ipzz:iqcr(n)=iqzz 'R coeff(n)=zz end if next i next n return 'DivPolynom(P,Q,V,R) Euclidian division: P(x)=Q(x)*V(x)+R(x) 'NUMBER: u 600 IERROR=1 'The Q polynomial must be <> zero if iqdeg=0 and cqv(0)=0 then return 'Divide by zero error 'put P(x) in R(x) irdeg=ipdeg for i=0 to irdeg isrcr(i)=isrcp(i):crv(i)=cpv(i) ipcr(i)=ipcp(i):iqcr(i)=iqcp(i) next i ivdeg=ipdeg - iqdeg if ivdeg<0 then ivdeg=0: isrcv(0)=0:cvv(0)=0#:ipcv(0)=0:iqcv(0)=1 else for ii=ivdeg to 0 step -1 'DivNumber(cr(irdeg),cq(iqdeg),cv(ii)) isrxx=isrcr(irdeg):xxv=crv(irdeg):ipxx=ipcr(irdeg):iqxx=iqcr(irdeg)'xx=cr(irdeg) isryy=isrcq(iqdeg):yyv=cqv(iqdeg):ipyy=ipcq(iqdeg):iqyy=iqcq(iqdeg)'yy=cq(iqdeg) gosub 3200 'zz=xx/yy if IERROR<>0 then return 'divide by zero error isrcv(ii)=isrzz:cvv(ii)=zzv:ipcv(ii)=ipzz:iqcv(ii)=iqzz 'cv(ii)=zz for jj=ii to irdeg 'MultNumber(cv(ii),cq[jj-ii), u) isrxx=isrcv(ii):xxv=cvv(ii):ipxx=ipcv(ii):iqxx=iqcv(ii) 'xx=cv(ii) k=jj-ii isryy=isrcq(k):yyv=cqv(k):ipyy=ipcq(k):iqyy=iqcq(k) 'yy=cq(jj-ii) gosub 3100 'zz=xx*yy isru=isrzz:uv=zzv:ipu=ipzz:iqu=iqzz 'u=zz ipu=-ipu: uv=-uv 'u=-u 'AddNumber(cr(jj),u, cr(jj)) isrxx=isrcr(jj):xxv=crv(jj):ipxx=ipcr(jj):iqxx=iqcr(jj) 'xx=cr(jj) isryy=isru:yyv=uv:ipyy=ipu:iqyy=iqu 'yy=u gosub 3000 'zz=xx+yy isrcr(jj)=isrzz:crv(jj)=zzv:ipcr(jj)=ipzz:iqcr(jj)=iqzz 'cr(jj)=zz next jj if irdeg > 0 then irdeg=irdeg-1 next ii 610 if abs(crv(irdeg))>=SMALL or irdeg<=0 then goto 620 irdeg=irdeg-1 goto 610 620 end if IERROR=0 return 'DivPolynom 'GCD of two polynomials 'GCDPolynom(P,Q,R) 'NUMBER: z 'integer ipg,ipp 700 IERROR=1 if ipdeg=0 and abs(cpv(ipdeg))ipdeg then irdeg=ipdeg 'R=P for i=0 to irdeg isrcr(i)=isrcp(i):crv(i)=cpv(i):ipcr(i)=ipcp(i):iqcr(i)=iqcp(i) next i ipdeg=iqdeg 'P=Q for i=0 to irdeg isrcp(i)=isrcq(i):cpv(i)=cqv(i):ipcp(i)=ipcq(i):iqcp(i)=iqcq(i) next i iqdeg=irdeg 'Q=R for i=0 to irdeg isrcq(i)=isrcr(i):cqv(i)=crv(i):ipcq(i)=ipcr(i):iqcq(i)=iqcr(i) next i end if irdeg=iqdeg 'R=Q for i=0 to irdeg isrcr(i)=isrcq(i):crv(i)=cqv(i):ipcr(i)=ipcq(i):iqcr(i)=iqcq(i) next i 710 if irdeg<=0 then goto 720 gosub 600 'DivPolynom(P,Q,V,R) if IERROR<>0 then return 715 if irdeg>0 and abs(crv(irdeg))SMALL then iqdeg=0 else ipdeg=iqdeg 'P=Q for i=0 to ipdeg isrcp(i)=isrcq(i):cpv(i)=cqv(i):ipcp(i)=ipcq(i):iqcp(i)=iqcq(i) next i iqdeg=irdeg 'Q=R for i=0 to irdeg isrcq(i)=isrcr(i):cqv(i)=crv(i):ipcq(i)=ipcr(i):iqcq(i)=iqcr(i) next i end if for i1=0 to iqdeg-1 'DivNumber(Q.coeff[i1],Q.coeff[Q.degree],Q.coeff[i1]) isrxx=isrcq(i1):xxv=cqv(i1):ipxx=ipcq(i1):iqxx=iqcq(i1) 'xx=cq(i1) i2=iqdeg isryy=isrcq(i2):yyv=cqv(i2):ipyy=ipcq(i2):iqyy=iqcq(i2) 'yy=cq(iqdeg) gosub 3200 'zz=xx/yy if IERROR<>0 then return 'divide by zero error isrcq(i1)=isrzz:cqv(i1)=zzv:ipcq(i1)=ipzz:iqcq(i1)=iqzz 'cq(i1)=zz next i1 isrcq(i2)=0:cqv(i2)=1#:ipcq(i2)=1:iqcq(i2)=1 'cq(iqdeg)=1 goto 710 720 irdeg=iqdeg 'R=Q for i=0 to irdeg isrcr(i)=isrcq(i):crv(i)=cqv(i):ipcr(i)=ipcq(i):iqcr(i)=iqcq(i) next i if iqdeg=0 then isrcr(0)=0:crv(0)=1#:ipcr(0)=1:iqcr(0)=1 'cr(0)=1 else ipg=0 for i=0 to irdeg if isrcr(i)=0 then 'R coeff(i) not real if ipg=0 then ipg=ipcr(i) else ia=ipg:ib=ipcr(i):gosub 900 'ic=GCD(ia,ib) ipg=ic end if end if next i ipp=0 for i=0 to irdeg if isrcr(i)=0 then 'R coeff(i) not real if ipp=0 then ipp=iqcr(i) else ia=ipp:ib=iqcr(i):gosub 900 'ic=GCD(ia,ib) ipp=ipp*iqcr(i)/ic end if end if next i if ipg<>0 then isrz=0: ipz=ipp iqz=ipg: zv=ipp/ipg for i=0 to irdeg 'MultNumber(R coeff(i), z, R coeff(i)) isrxx=isrcr(i):xxv=crv(i):ipxx=ipcr(i):iqxx=iqcr(i) 'xx=cr(i) isryy=isrz:yyv=zv:ipyy=ipz:iqyy=iqz 'yy=z gosub 3100 'zz=xx*yy isrcr(i)=isrzz:crv(i)=zzv:ipcr(i)=ipzz:iqcr(i)=iqzz 'cr(i)=zz next i end if end if IERROR=0 return 'SCM of two polynomials 'SCMPolynom(P,Q,R) 'POLYNOMIALS: U,V 750 IERROR=1 'GCDPolynom(P,Q,U) gosub 700 'GCDPolynom(P,Q,R) 'Note: P(x) and Q(x) are modified during the process if IERROR<>0 then return 'Error in GCDPolynom iudeg=irdeg 'U=R for i=0 to iudeg isrcu(i)=isrcr(i):cuv(i)=crv(i):ipcu(i)=ipcr(i):iqcu(i)=iqcr(i) next i 'restore P(x) from P1(x) ipdeg=ipdeg1 for i=0 to MAXPOL isrcp(i)=isrcp1(i):cpv(i)=cpv1(i):ipcp(i)=ipcp1(i):iqcp(i)=iqcp1(i) next i 'restore Q(x) from Q1(x) iqdeg=iqdeg1 for i=0 to MAXPOL isrcq(i)=isrcq1(i):cqv(i)=cqv1(i):ipcq(i)=ipcq1(i):iqcq(i)=iqcq1(i) next i 'MultPolynom(P,Q,V) gosub 500 'MultPolynom(P,Q,R) ivdeg=irdeg 'V=R for i=0 to ivdeg isrcv(i)=isrcr(i):cvv(i)=crv(i):ipcv(i)=ipcr(i):iqcv(i)=iqcr(i) next i 'DivPolynom(V,U,R,V) ipdeg=ivdeg 'P=V for i=0 to ipdeg isrcp(i)=isrcv(i):cpv(i)=cvv(i):ipcp(i)=ipcv(i):iqcp(i)=iqcv(i) next i iqdeg=iudeg 'Q=U for i=0 to iqdeg isrcq(i)=isrcu(i):cqv(i)=cuv(i):ipcq(i)=ipcu(i):iqcq(i)=iqcu(i) next i gosub 600 'DivPolynom(P,Q,V,R) if IERROR<>0 then return 'Error in DivPolynom irdeg=ivdeg 'R=V for i=0 to irdeg isrcr(i)=isrcv(i):crv(i)=cvv(i):ipcr(i)=ipcv(i):iqcr(i)=iqcv(i) next i IERROR=0 return $INCLUDE "Numeric\Polynoms\Polynoms.bas" 'end of file gcdpol.bas