'*****************************************************
'*     Euclidian division of two polynomials         *
'*             R(x) = P(x) / Q(x)                    *
'* ------------------------------------------------- *
'* Ref.: "Mathematiques 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:                                         *
'*                                                   *
'* 2 X2 + 2                                          *
'*                                                   *
'* Remainder:                                        *
'*                                                   *
'*  X2 + 2 X + 1                                     *
'*                                                   *
'*                                                   *
'*                     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 " EUCLIDIAN DIVISION OF TWO POLYNOMIALS:"
  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)
  '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 DivPolynom(P,Q,H,R)

  if IERROR<>0 then
    print " Error in euclidian division"
  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


'DivPolynom(P,Q,H,R)
'NUMBER: u
500 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
  ihdeg=ipdeg - iqdeg
  if ihdeg<0 then
    ihdeg=0:
    isrch(0)=0:chv(0)=0#:ipch(0)=0:iqch(0)=1
  else
    for ii=ihdeg to 0 step -1
      'DivNumber(cr(irdeg),cq(iqdeg),ch(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
      isrch(ii)=isrzz:chv(ii)=zzv:ipch(ii)=ipzz:iqch(ii)=iqzz   'ch(ii)=zz
      for jj=ii to irdeg
        'MultNumber(ch(ii),cq[jj-ii), u)
        isrxx=isrch(ii):xxv=chv(ii):ipxx=ipch(ii):iqxx=iqch(ii) 'xx=ch(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
510 if abs(crv(irdeg))>=SMALL or irdeg<=0 then goto 520
      irdeg=irdeg-1
      goto 510
520 end if
  IERROR=0
return 'DivPolynom


$INCLUDE "numeric\polynoms\polynoms.bas"

'end of file divpol.bas