'****************************************************** '* 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 ia0 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)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