'******************************************************
'*         Basic subroutines for polynomials          *
'* -------------------------------------------------- *
'*  Ref.: "Mathematiques en Turbo-Pascal by M. Ducamp *
'*  and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
'*  [BIBLI 05].                                       *
'*                                                    *
'*                BASIC version by J-P Moreau, Paris. *
'*                         (www.jpmoreau.fr)          *
'* -------------------------------------------------- *
'* Version 1.1:   1. added safeguards against integer *
'*                overflow in AddNumber, MultNumber.  *
'*                2. correct bug in subroutine 4100   *
'*                in setting P(x) to zero.            *
'******************************************************
'
'(defined by calling program)
'MAXPOL=20  'Maximum degree for a polynomial
'SMALL=1e-20

' a number zz (real, integer or fractional) is defined by:
' 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, double, integer or fractional)

'predefine 3 polynomials P,Q,R (made by calling program)
'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)


'Copy number xx in zz
800
  isrzz=isrxx:zzv=xxv:ipzz=ipxx:iqzz=iqxx
return

'Copy number zz in xx
801
  isrxx=isrzz:xxv=zzv:ipxx=ipzz:iqxx=iqzz
return

'Copy number yy in zz
802
  isrzz=isryy:zzv=yyv:ipzz=ipyy:iqzz=iqyy
return

'Copy number zz in yy
803
  isryy=isrzz:yyv=zzv:ipyy=ipzz:iqyy=iqzz
return

'return greater common divisor of two integers
'ic = GCD(ia, ib)
900 ia=abs(ia)
ib=abs(ib)
if ia=0 or ib=0 then
  ic=1
  return
end if
if ia<ib then
  ir=ia: ia=ib: ib=ir
end if
ir=1
910 ir=ia MOD ib
ia=ib: ib=ir
if ir>0 then goto 910
ic=ia
return


'Convert a string ch$ into a number zz of type double, integer or fractional
'Examples of valid strings: 5.56874 or 15/187 or 255
'zz=SetNumber(ch$)
'integer ir, string temp$, integer isign
1000 isrzz=0: isfract=0: temp$="": isign=1
if LEFT$(ch$,1)="-" then isign=-1
if LEFT$(ch$,1)="-" or LEFT$(ch$,1)="+" then
'remove first character - or +
  ch$=RIGHT$(ch$,LEN(ch$)-1)
end if
for i=1 to LEN(ch$)
  if MID$(ch$,i,1)="." then isrzz=1    'detect a dot in string ch$
next i
if isrzz<>0 then
  zzv=VAL(ch$)
else
  for i=1 to LEN(ch$)
    if (MID$(ch$,i,1)="/") then
      ipos=i: isfract=1             'detect / in string ch$ at position i
      goto 1010
    end if
  next i
1010 if isfract<>0 then                'case of a fractional number
  temp$=LEFT$(ch$,ipos-1)              'take part before / & evaluate
  ipzz=VAL(temp$)                      'numerator of zz
  temp$=MID$(ch$,ipos+1,LEN(ch$)-ipos) 'take part after / & evaluate
  iqzz=VAL(temp$)
'simplify ipzz and iqzz ?
  ia=ipzz: ib=iqzz
  gosub 900                            'ic=GCD(ipzz,iqzz)
  ipzz=ipzz/ic: iqzz=iqzz/ic
  zzv=ipzz/iqzz
  else                                 'case of an integer
    zzv=VAL(ch$)
    ipzz=VAL(ch$): iqzz=1
  end if
end if
'case sign is negative
if isign=-1 then
  zzv=-zzv
  ipzz=-ipzz
end if
IERROR=0  'all ok
return 'SetNumber


'write to screen a number zz of type double, integer or fractional
2000 'Subroutine WriteNumber(zz)
if (isrzz<>0 and zzv<0) or ipzz<0 then
  print "-";
else
  print "+";
end if
if isrzz<>0 then   'real number
  print abs(zzv)
else               'integer or fractional number
  print abs(ipzz);
  if abs(iqzz)<>1 and ipzz<>0 then print "/"; abs(iqzz)
end if
return 'WriteNumber


'Add two numbers of type ar_number (double, integer or fractional)
'zz=AddNumber(xx, yy)
'integer negative
3000 isrzz=isrxx + isryy
zzv=xxv+yyv
if xxv=0 then
  gosub 802    'zz=yy
else
  if yyv=0 then
    gosub 800  'zz=xx
  else
    if abs(ipxx) > MAXINT or abs(iqxx) > MAXINT then isrzz=1
    if abs(ipyy) > MAXINT or abs(iqyy) > MAXINT then isrzz=1
    if isrzz=0 then
      ia=iqxx: ib=iqyy: gosub 900:        'ic=GCD(iqxx,iqyy)
      if 1#*ipxx*iqyy+1#*iqxx*ipyy>MAXINT or 1#*iqxx*iqyy>MAXINT then
        isrzz=1            'real number instead
      else
        ipzz=ipxx*iqyy+iqxx*ipyy: iqzz=iqxx*iqyy
        if 1#*ipzz*iqzz<0 then
          negative=1
        else
	  negative=0
        end if
        ipzz=abs(ipzz) / ic
        iqzz=abs(iqzz) / ic
        if negative<>0 then ipzz=-ipzz
      end if
    end if
  end if
end if
IERROR=0
return 'AddNumber


'Multiply two numbers of type ar_number (double, integer or fractional)
'Function zz=MultNumber(xx,yy)
'integer negative
3100 isrzz = isrxx + isryy
zzv=xxv * yyv
if isrzz<>0 and zzv=0 then
  ipzz=0: iqzz=1
else
  if abs(ipxx) > MAXINT or abs(iqxx) > MAXINT then isrzz=1
  if abs(ipyy) > MAXINT or abs(iqyy) > MAXINT then isrzz=1
  if isrzz=0  then           'integer or fractional number
    if 1#*ipxx*ipyy>MAXINT or 1#*iqxx*iqyy>MAXINT then
      isrzz=1                             'real number instead
    else
      ipzz=ipxx*ipyy: iqzz=iqxx*iqyy
      ia=ipzz: ib=iqzz: gosub 900         'ic=GCD(ipzz,iqzz)
      if 1#*ipzz*iqzz<0 then
        negative=1
      else
        negative=0
      end if
      ipzz=abs(ipzz) / ic
      iqzz=abs(iqzz) / ic
      if negative<>0  then ipzz=-ipzz
    end if
  end if
end if
IERROR=0
return 'MultNumber

'Divide two numbers of type NUMBER (real, integer or fractional)
'DivNumber: zz=xx/yy
3200 IERROR=2
  if isryy<>0 and abs(yyv)<SMALL then return     'Divide by zero error
  if isryy=0 and abs(ipyy)=0 then return         'IERROR=2
  isrzz=isrxx+isryy
  if (abs(ipxx)>MAXINT) or (abs(iqxx)>MAXINT) then isrzz=1
  if (abs(ipyy)>MAXINT) or (abs(iqyy)>MAXINT) then isrzz=1
  if isrzz=0 then
    if 1#*ipxx*iqyy>MAXINT or 1#*iqxx*ipyy>MAXINT then
      isrzz=1
    else
      ipzz=ipxx*iqyy: iqzz=iqxx*ipyy
      'ic=GCD(z.p,z.q)
      ia=ipzz:ib=iqzz:gosub 900  'ic=GCD(ia,ib)
      if 1#*ipzz*iqzz<0 then
        negative=1
      else
        negative=0
      end if
      ipzz=abs(ipzz)/ic: iqzz=abs(iqzz)/ic
      if negative=1 then ipzz=-ipzz
    end if
  end if
  zzv=xxv/yyv
  IERROR=0
return 'DivNumber

'Subroutine called by EnterPolynom
'Analysis of elementary monom, valid examples: +5.25x3 or +5x3 or 5/3x3
'Subroutine ExtractMonom(ch$, p)
'p is a polynomial
'type number: coeff,signe0
'temporary string chn$
'integer ideg,ixpos,isppos,isigne
4100 chn$=""
'eliminate initial blank spaces in ch$
4101 if LEFT$(ch$,1)=" " then
    ch$=RIGHT$(ch$,Len(ch$)-1)
  else
    goto 4102
  end if
  goto 4101
4102 if LEN(ch$)<1 then
  IERROR=1
  return
end if
if LEFT$(ch$,1)="-" then
  isigne=-1
else
  isigne=1
end if
if LEFT$(ch$,1)="+" or LEFT$(ch$,1)="-" then
  ch$=RIGHT$(ch$,LEN(ch$)-1)
end if
if LEN(ch$)<1 then
  IERROR=1
  return
end if
ixpos=-1
'seek x characters
for i=1 to LEN(ch$)
  if MID$(ch$,i,1)="x" then
    ixpos=i: goto 4110
  end if
next i
4110 if ixpos<0 then          'no more x found = constant term
  gosub 1000                                                 'zz=SetNumber(ch$)
  isrcp(0)=isrzz: cpv(0)=zzv: ipcp(0)=ipzz: iqcp(0)=iqzz     'cp(0)=zz
  if isigne=-1 then
    ipcp(0)=-ipcp(0)
    cpv(0)=-cpv(0)
  end if
  ch$=""
  return
else  'term with x found
  if ixpos=1 then
    chn$="1"
  else
    chn$=MID$(ch$,1,ixpos-1)
  end if
  'evaluate chn$
  ch1$=ch$ 'save ch$ in ch1$
  ch$=chn$: gosub 1000  'zz=SetNumber(ch$)
  isrcoeff=isrzz: coeffv=zzv: ipcoeff=ipzz: iqcoeff=iqzz  'coeff=zz
  ch$=ch1$  'restore ch$
  'seek a blank space
  isppos=-1
  for i=1 to LEN(ch$)
    if MID$(ch$,i,1)=" " then
      isppos=i
      goto 4120
    end if
  next i 		
4120 if isppos<0  then isppos=LEN(ch$)+1
  if ixpos=isppos-1 then
    ideg=1
  else
    chn$=MID$(ch$,ixpos+1,isppos-ixpos)
    ideg=VAL(chn$)
    if ideg > MAXPOL then
      IERROR=1
      return
    end if
  end if
  if isppos > LEN(ch$) then
    ch$=" "
  else
    ch$=RIGHT$(ch$,LEN(ch$)-isppos)
  end if
end if
if isigne=-1 then
  ipcoeff = -ipcoeff
  coeffv = -coeffv
end if
'AddNumber(coeff, cp(ideg),cp(ideg))
isrxx=isrcoeff: xxv=coeffv: ipxx=ipcoeff: iqxx=iqcoeff             'xx=coeff
isryy=isrcp(ideg): yyv=cpv(ideg): ipyy=ipcp(ideg): iqyy=iqcp(ideg) 'yy=cp(ideg)
gosub 3000 'zz=xx+yy
isrcp(ideg)=isrzz: cpv(ideg)=zzv: ipcp(ideg)=ipzz: iqcp(ideg)=iqzz 'cp(ideg)=zz
if ideg > ipdeg  then ipdeg=ideg
IERROR=0
return 'ExtractMonom


'Convert a valid string into a polynomial P
'Example of valid string: ch$="x5 +3/5x4 -12x2 +8x -1/4"
'tx$ is a prompting message, example: tx$="Enter polynomial P(x):"
'Note that a blank space is required between monomes.
'integer Function EnterPolynom(tx$, p)
4000 'prompt string: tx$, polynomial: p
'string defining polynomial: ch$
ch$=""
'put polynomial P(x) to zero
ipdeg=0
for i=0 to MAXPOL
  isrcp(i)=0:cpv(i)=0#:ipcp(i)=0:iqcp(i)=1
next i
for i=1 to 3   '3 tries authorized
  print tx$;
  input "",ch$
'suppress initial blanks
4010 if LEFT$(ch$,1)=" " then
  ch$=RIGHT$(ch$,Len(ch$)-1)
  else
    goto 4020
  end if
  goto 4010
4020 gosub 4100        'call ExtractMonom(ch$,p)
  if Len(ch$)<1 then   'string is completely analysed
    IERROR=0           'with no error
    return
  else
    goto 4020          'string is not finished, continue
  end if	
next i
'something went wrong
IERROR=1
return 'EnterPolynom


'This subroutine displays to screen the symbolic
'representation of a polynomial P(x)
'Subroutine DisplayPolynom(p)
5000 print
  if ipdeg=0 and cpv=0 then
    print "0"
  else
    for i=ipdeg to 0 step -1
      if abs(cpv(i)) > SMALL or ipcp(i)<>0 then  'print if coeff <>0
	if (isrcp(i)<>0 and cpv(i)<0) or ipcp(i)<0 then
	  print "-";
        else
	  print "+";
        end if
        if isrcp(i)<>0 then 'case of a real number
	  print using "####.####"; abs(cpv(i));
        else                'case of an integer or fractional number
	  lp=abs(ipcp(i))
          lq=abs(iqcp(i))
	  if lp<>1 or lq<>1 or i=0  then print lp;
	  if lq<>1 and lp<>0 then print  "/"; lq;
        end if
        'print X+exponent
        if i<>0 then
          if i<>1 then
            print " X"; i; " ";
          else
            print " X ";
          end if
        end if
      end if
    next i
  end if
  print
return 'DisplayPolynom

'return the value y of a polynomial P(x) for argument x=zz
'x and y are of type number (real, integer or fractional)
'Function y = EvalPolynom(p,zz)
6000 if ipdeg > MAXPOL then
    IERROR=1
    return
  end if
  'save initial zz in z0
  isrz0=isrzz:z0v=zzv:ipz0=ipzz:iqz0=iqzz
  'set y to zero
  isry=0: yv=0: ipy=0: iqy=1
  for i=0 to ipdeg
    'y=y*z0
    isrxx=isrz0:xxv=z0v:ipxx=ipz0:iqxx=iqz0  'xx=z0
    isryy=isry:yyv=yv:ipyy=ipy:iqyy=iqy      'yy=y
    gosub 3100                               'zz=xx*yy
    isry=isrzz:yv=zzv:ipy=ipzz:iqy=iqzz      'y=zz
    'y=y+cp(ipdeg-i)
    j=ipdeg-i
    isrxx=isrcp(j):xxv=cpv(j):ipxx=ipcp(j):iqxx=iqcp(j)  'xx=cp(ipdeg-i)
    isryy=isry:yyv=yv:ipyy=ipy:iqyy=iqy      'yy=y
    gosub 3000                               'zz=xx+yy
    isry=isrzz:yv=zzv:ipy=ipzz:iqy=iqzz      'y=zz
  next i
  IERROR=0
return 'EvalPolynom

'--------------------------------------------------------
'NOTE: this file must be included by any program using
'      polynomials at the end of it, using the line:
'
'      $INCLUDE polynoms.bas
'
'      (this file is supposed to be in the same directory
'       than the main program).
'--------------------------------------------------------
'end of file Polynoms.bas