c ---------------------------------------------------------------------- c Menu-base cubic equation solver c c(1)*x**3 + c(2)*x**2 + c(3)*x + c(4) = 0 c c Same structure as MENUREA2.FOR c Variables used (cubic equation related) c c(4) ... coefficients c x(3) ... solutions c ireal ... number of real roots c nroot ... number of roots c Variables used (menu related) c ipick ... the field picked c irow(i) ... the row number of the ith field c icol(i) ... the column number of ith field c ifield =1..4 ... coefficients c ifield =5..6 ... number of real and complex solutions c ifield =7..9 ... answers c ascii,iscan ... used in PEAKKY & READKY to report the ASCII character and c scancodes c scan ... used to associate scancodes with key names c ---------------------------------------------------------------------- c Programming Note: c Link with asm.obj c ---------------------------------------------------------------------- c Instructor: Nam Sun Wang c ---------------------------------------------------------------------- integer irow(9),icol(9) real c(4) complex x(3) character scan(152)*2,ascii*1 character texti*2 ,texti1*1 character textr*13,textr1*12 character textc*27,textc1*26,blank*26 equivalence (texti1,texti) equivalence (textr1,textr) equivalence (textc1,textc) c The following definition makes easier the testing of scancodes ------- data scan/'es','11','22','33','44','55','66','77','88','99', 1- 10 1 '00','--','==','bs','tb','QQ','WW','EE','RR','TT', 11- 20 2 'YY','UU','II','OO','PP','[[',']]','cr','cl','AA', 21- 30 3 'SS','DD','FF','GG','HH','JJ','KK','LL',';;','''''', 31- 40 4 '``','ls','\\','ZZ','XX','CC','VV','BB','NN','MM', 41- 50 5 ',,','..','//','rs','ps','al',' ','ca','F1','F2', 51- 60 6 'F3','F4','F5','F6','F7','F8','F9','F0','nl','sl', 61 70 7 'HM','UP','PU','N-','LT','CT','RT','N+','EN','DN', 71 80 8 'PD','IN','DE','S1','S2','S3','S4','S5','S6','S7', 81 90 9 'S8','S9','S0','C1','C2','C3','C4','C5','C6','C7', 91 100 A 'C8','C9','C0','A1','A2','A3','A4','A5','A6','A7', 101-110 B 'A8','A9','A0','CP','CL','CR','CE','CD','CH','a1', 111-120 C 'a2','a3','a4','a5','a6','a7','a8','a9','a0','A-', 121-130 D 'A=','CU','##','C-','##','C+','##','CI','C.','AH', 131-140 E '##','AU','A-','AL','##','AR','A+','AE','##','AD', 141-150 F 'AI','A.'/ 151-152 c Add some dollar signs ------------------------------------------------ texti( 2: 2) = '$' textr(13:13) = '$' textc(27:27) = '$' blank=' ' c Set up variables ----------------------------------------------------- iattr = 112 ipick = 1 minpick = 1 maxpick = 4 c Display the menu ----------------------------------------------------- call menu c Rows and columns corresponding to each field ------------------------- irow(1) = 7 irow(2) = 9 irow(3) = 11 irow(4) = 13 irow(5) = 17 irow(6) = 17 irow(7) = 19 irow(8) = 21 irow(9) = 23 icol(1) = 30 icol(2) = 30 icol(3) = 30 icol(4) = 30 icol(5) = 25 icol(6) = 34 icol(7) = 30 icol(8) = 30 icol(9) = 30 c Set and display default values in fields 1 through maxpick ----------- do i=1, maxpick c(i)=1. write(textr1,602) c(i) call outpta(irow(i), icol(i), iattr, textr) enddo c Display answers in fields 5 through 9 -------------------------------- 100 call cubic(c(1), c(2), c(3), c(4), x, nroot) c Find the number of real roots ireal = 0 do i=1, nroot if(aimag(x(i)) .eq. 0.) ireal=ireal+1 enddo write(texti1,501) ireal call outpta(irow(5), icol(5), 7, texti) write(texti1,501) nroot-ireal call outpta(irow(6), icol(6), 7, texti) c Output solutions (actual ones followed by blank spaces for non-existent ones) do i=1, nroot write(textc1,650)x(i) call outpta(irow(i+6), icol(i+6), iattr, textc) enddo do i=nroot+1, 3 textc1=blank call outpta(irow(i+6), icol(i+6), iattr, textc) enddo c----------------------------------------------------------------------- c Wait for response ---------------------------------------------------- 200 call posit(irow(ipick),icol(ipick)) 50 call peekky(ascii,iscan,iyesno) if(iyesno .eq. 0)goto 50 c Proceed to check to see which key is pressed ------------------------- if(scan(iscan).eq.'F0')then call readky(ascii,iscan) c Exit to DOS call clrscr stop elseif(scan(iscan).eq.'UP')then c Go up one field ipick=ipick-1 if(ipick .eq. minpick-1)ipick=maxpick call readky(ascii,iscan) elseif(scan(iscan).eq.'DN' .or. scan(iscan).eq.'cr')then c Go down one field ipick=ipick+1 if(ipick .eq. maxpick+1)ipick=minpick call readky(ascii,iscan) elseif(ichar(ascii) .ne. 0 .and. scan(iscan).ne.'cr')then c Any other ascii character will be read into the variable call READX(irow(ipick),icol(ipick),iattr,c(ipick),-1.e30,1.e30) ipick=ipick+1 if(ipick .eq. maxpick+1)ipick=minpick goto 100 else c Any other key gets recycled back without any change in IPICK call readky(ascii,iscan) goto 200 endif goto 200 c Standard formats ----------------------------------------------------- 500 format(a) 501 format(i1) 602 format(1pe12.5) 650 format(1pe13.5,e13.5) end c ---------------------------------------------------------------------- subroutine menu cÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ Cubic Equation Solver By Nam Sun Wang ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ» cº º cº A*X**3 + B*X**2 + C*X + D = 0 º cº º cÌÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ͹ cº º cº COEFFICIENTS *** A: | º cº º cº B: | º cº º cº C: | º cº º cº D: | º cº º cÌÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ͹ cº º cº SOLUTIONS ****** ? Real & ? Complex Roots º cº º cº X1: 123456789012 123456789012 º cº º cº X2: | º cº º cº X3: | º cº º cÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍPress F10 to Exit to DOSÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ call clrscr call output( 1, 1,'ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ Cubic Equation Solver By N *am Sun Wang ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»$') call output( 2, 1,'º * º$') call output( 3, 1,'º A*X**3 + B*X**2 + C*X *+ D = 0 º$') call output( 4, 1,'º * º$') call output( 5, 1,'ÌÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ *ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ͹$') call output( 6, 1,'º * º$') call output( 7, 1,'º COEFFICIENTS *** A: * º$') call output( 8, 1,'º * º$') call output( 9, 1,'º B: * º$') call output(10, 1,'º * º$') call output(11, 1,'º C: * º$') call output(12, 1,'º * º$') call output(13, 1,'º D: * º$') call output(14, 1,'º * º$') call output(15, 1,'ÌÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ *ÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ͹$') call output(16, 1,'º * º$') call output(17, 1,'º SOLUTIONS ****** Real & Complex Root *s º$') call output(18, 1,'º * º$') call output(19, 1,'º X1: * º$') call output(20, 1,'º * º$') call output(21, 1,'º X2: * º$') call output(22, 1,'º * º$') call output(23, 1,'º X3: * º$') call output(24, 1,'º * º$') call output(25, 1,'ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍPress F10 to Exit t *o DOSÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ$') end c ---------------------------------------------------------------------- subroutine readx(irow,icol,iattr,xx,xmin,xmax) c ---------------------------------------------------------------------- c Read an integer ranged between xmin and xmax, including xmin and xmax. c Default/current values are displayed, and the new input is echoed back. c Automatically add the decimal point for a real number if it is missing. c The error message is displayed on the 25th row. c Public Variables ... c irow,icol ... row and column to display the message c iattr ... attribute (flashing, invert, bold, etc.) c xx ... integer variable whose value is to be read c xmin,xmax ... minimum and maxmum value that ix can take c If xmin and xmax are the same or xmin>xmax, accept all values. c Local Variables c fmt ... stores the format for output c textr, textr1 ... text representation of the integer variable ix c blank ... blank characters used for erase the trailing spaces c atext ... store input from the keyboard c nspace ... number of spaces in the field c ierror ... error flag c xxdef ... default value of ix c ilast ... number of non-blank spaces in the keyboard input c The field has 12 spaces. Change the number on the next line. c ---------------------------------------------------------------------- parameter (nspace=12,n1space=13,n=80) c character textr*n1space,textr1*nspace,blank*11,atext*n,fmt*40 character textr*13,textr1*12,blank*11,atext*80,fmt*40 equivalence (textr1,textr) write(fmt,650)nspace blank=' $' textr(n1space:n1space)='$' ierror = 0 c Save the current value to be used as default ------------------------- xxdef=xx c Echo the default value on screen ------------------------------------- write(textr1,fmt)xx 901 call outpta(irow,icol,iattr,textr) c Position the cursor to the desired (irow,icol) ----------------------- call posit(irow,icol) c Read the value into a text variable first ---------------------------- read(*,500)atext c Clear the error message, if it exists, after reading a new value ----- if(ierror .eq. 1)then call clrdwn(25,1) ierror=0 endif c Scan backward to find the first non-blank character ------------------ do 10 i=n,1,-1 if(atext(i:i).ne.char(0) .and. atext(i:i).ne.char(32))goto 11 10 continue return 11 ilast=i c Continue to scan backward to find the decimal point ------------------ do 20 i=ilast,1,-1 if(atext(i:i) .eq. '.')goto 21 20 continue c Append the decimal point if it has not already been entered ---------- atext(ilast+1:ilast+1)='.' c v3.20 and v3.20 compiler bug: cannot read ATEXT from the keyboard if the READ c statement is followed by another READ statement that read another variable c from an internal file ATEXT. For example, the following 2 lines give an c error message that the EOF is encountered: c read(*,500)atext c read(atext,*)y c The following line is a trick to overcome the above problem. 21 write(*,600)atext(n:n) c Assign the read value into the real variable XX. --------------------- c Use the default value if there is no entry. read(atext,'(e15.7)',err=901)xx c If xmin0, one real and two complex roots. c Step 3a: For D>0 and D=0, c Calculate u and v c u = cubic_root(-q/2 + sqrt(D)) c v = cubic_root(-q/2 - sqrt(D)) c Find the three transformed roots c y1 = u + v c y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2 c y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2 c Step 3a Alternately, for D<0, a trigonometric formulation is more convenient c y1 = 2 * sqrt(|p|/3) * cos(phi/3) c y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3) c y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3) c where phi = acos(-q/2/sqrt(|p|**3/27)) c pi = 3.141592654... c Step 4 Finally, find the three roots c x = y - b/a/3 c ---------------------------------------------------------------------- c Declare variables complex x(1) data pi/3.141592654/ c Step 0: If a is 0 use the quadratic formula. ------------------------- if(a .eq. 0.)then call quad(b, c, d, x, nroot) return endif c Cubic equation with 3 roots nroot = 3 c Step 1: Calculate p and q -------------------------------------------- p = c/a - b*b/a/a/3. q = (2.*b*b*b/a/a/a - 9.*b*c/a/a + 27.*d/a) / 27. c Step 2: Calculate DD (discriminant) ---------------------------------- DD = p*p*p/27. + q*q/4. c Step 3: Branch to different algorithms based on DD ------------------- if(DD .lt. 0.)then c Step 3b: c 3 real unequal roots -- use the trigonometric formulation phi = acos(-q/2./sqrt(abs(p*p*p)/27.)) temp1=2.*sqrt(abs(p)/3.) y1 = temp1*cos(phi/3.) y2 = -temp1*cos((phi+pi)/3.) y3 = -temp1*cos((phi-pi)/3.) else c Step 3a: c 1 real root & 2 conjugate complex roots OR 3 real roots (some are equal) temp1 = -q/2. + sqrt(DD) temp2 = -q/2. - sqrt(DD) u = abs(temp1)**(1./3.) v = abs(temp2)**(1./3.) if(temp1 .lt. 0.) u=-u if(temp2 .lt. 0.) v=-v y1 = u + v y2r = -(u+v)/2. y2i = (u-v)*sqrt(3.)/2. endif c Step 4: Final transformation ----------------------------------------- temp1 = b/a/3. y1 = y1-temp1 y2 = y2-temp1 y3 = y3-temp1 y2r=y2r-temp1 c Assign answers ------------------------------------------------------- if(DD .lt. 0.)then x(1) = cmplx( y1, 0.) x(2) = cmplx( y2, 0.) x(3) = cmplx( y3, 0.) elseif(DD .eq. 0.)then x(1) = cmplx( y1, 0.) x(2) = cmplx(y2r, 0.) x(3) = cmplx(y2r, 0.) else x(1) = cmplx( y1, 0.) x(2) = cmplx(y2r, y2i) x(3) = cmplx(y2r,-y2i) endif end