'******************************************************* '* Division of two polynomials by increasing powers * '* k * '* P(x) = Q(x).H(x) + x.R(x) * '* --------------------------------------------------- * '* Ref.: "Mathematiques 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 + 2 X2 + X * '* * '* Remainder: * '* * '* - X2 - X + 1 * '* * '* * '* DIVISION OF TWO POLYNOMIALS BY INCREASING POWERS: * '* * '* P(x) = 3/5x3 +2x -4 * '* Q(x) = x2 -5/4 * '* Order: 9 * '* * '* * '* Quotient: * '* * '* 4096/3125 X8 - 704/625 X7 + 1024/625 X6 -176/125 X5 * '* + 256/125 X4 - 44/25 X3 + 64/25 X2 - 8/5 X + 16/5 * '* * '* Remainder: * '* * '* - 4096/3125 X + 704/625 * '* * '* * '* BASIC Version By J-P Moreau. * '* (www.jpmoreau.fr) * '******************************************************* 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,H 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) DIM isrch(MAXPOL),chv(MAXPOL),ipch(MAXPOL),iqch(MAXPOL) cls print print " DIVISION OF TWO POLYNOMIALS BY INCREASING POWERS:" print tx$=" P(x) = ": gosub 4000 'Enter P(x) 'save 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 tx$=" Q(x) = ": gosub 4000 'Enter P(x) print input " Order: ", iorder 'Note: name is k in title 'put P(x) in Q(x) iqdeg=ipdeg for i=0 to iqdeg isrcq(i)=isrcp(i) cqv(i)=cpv(i) ipcq(i)=ipcp(i) iqcq(i)=iqcp(i) next i 'restore P(x) from R(x) 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 gosub 500 'call DivPolynom1(P,Q,iorder,H,R) if IERROR<>0 then print " Error in division by increasing powers." else 'put H(x) in P(x) for printing ipdeg=ihdeg for i=0 to ipdeg isrcp(i)=isrch(i) cpv(i)=chv(i) ipcp(i)=ipch(i) iqcp(i)=iqch(i) next i print print " Quotient:":gosub 5000 'display P(x) print 'put 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 print " Remainder:":gosub 5000 'display P(x) end if END 'of main program 'Division of two polynomials by increasing powers 'DivPolynom1(P,Q,iorder,H,R) 'NUMBER: u 500 IERROR=1 'The Q polynomial constant term and iorder must be <> zero if iorder<=0 then return if isrcq(0)<>0 and abs(cqv(0))iqdegree then for i=iqdeg+1 to mx isrcq(i)=0:cqv(i)=0#:ipcq(i)=0:iqcq(i)=1 'cq(i)=0 next i end if 'set R(x) to 0 irdeg=0 for i=0 to MAXPOL isrcr(i)=0:crv(i)=0#:ipcr(i)=0:iqcr(i)=1 next i for i=0 to ipdeg 'cr(i)=cp(i) isrcr(i)=isrcp(i):crv(i)=cpv(i):ipcr(i)=ipcp(i):iqcr(i)=iqcp(i) next i if iorder>iqdeg then for i=iqdeg+1 to iorder isrcq(i)=0:cqv(i)=0#:ipcq(i)=0:iqcq(i)=1 'cq(i)=0 next i end if for ii=0 to iorder-1 'DivNumber(R coeff(ii),Q coeff(0),H coeff(ii)) isrxx=isrcr(ii):xxv=crv(ii):ipxx=ipcr(ii):iqxx=iqcr(ii) 'xx=cr(ii) isryy=isrcq(0):yyv=cqv(0):ipyy=ipcq(0):iqyy=iqcq(0) 'yy=cq(0) gosub 3200 'zz=xx/yy if IERROR<>0 then print "Error in DivNumber: ierror=";IERROR return end if isrch(ii)=isrzz:chv(ii)=zzv:ipch(ii)=ipzz:iqch(ii)=iqzz 'ch(ii)=zz for jj=1 to mx 'MultNumber(H.coeff[ii],Q.coeff[jj-ii], u) isrxx=isrch(ii):xxv=chv(ii):ipxx=ipch(ii):iqxx=iqch(ii) 'xx=ch(ii) kk=jj-ii if kk<0 then 'BASIC does not accept negative index isryy=0:yyv=0#:ipyy=0:iqyy=1 else isryy=isrcq(kk):yyv=cqv(kk):ipyy=ipcq(kk):iqyy=iqcq(kk) 'yy=cq(jj-ii) end if gosub 3100 'zz=xx*yy isru=isrzz:uv=zzv:ipu=ipzz:iqu=iqzz 'u=zz ipu=-ipu: uv=-uv 'u=-u 'AddNumber(R.coeff[jj], u, R.coeff[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 next ii ihdeg=iorder-1 510 if ihdeg<=0 or abs(chv(ihdeg)) > SMALL then goto 520 ihdeg=ihdeg-1 goto 510 520 irdeg=mx-iorder for i=0 to irdeg 'R coeff(i)=R coeff(i+iorder) isrcr(i)=isrcr(i+iorder):crv(i)=crv(i+iorder) ipcr(i)=ipcr(i+iorder):iqcr(i)=iqcr(i+iorder) next i 530 if irdeg<=0 or abs(crv(irdeg)) > SMALL then goto 540 irdeg=irdeg-1 goto 530 540 IERROR=0 return $INCLUDE "numeric\polynoms\polynoms.bas" 'end of file divpol1.bas