Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmod69.bas
There are 2 other files named cmod69.bas in the archive. Click here to see a list.
00010Q9=0
00020REM CMOD69 CMOD69 CMOD69 CMOD69
00025REM ANOVA CONTROL
00035 DIM M(10),N(10),S(10),A(10),C(10) ,P(6,6),Q(6),U(1),G(1,10)
00040DIML(10,10),I(10),V(10),T(10,10),B(10),F(10)
00045 C$="111213212223313233"
00050GOSUB 5160
00055 FILES RFILE1,RFILE2,RFILE3,RF4,,,RF7,RF8
00075RESTORE#1
00076 INPUT# 1,I1,I2,I3
00080SCRATCH#1
00081 PRINT # 1,69,I2,I3
00085I9=I2
00090I2=I3
00095S7=0
00100X=0
00105S9=0
00110IF I1=89 THEN 135
00115IF I1=75 THEN 135
00120RESTORE#8
00121 INPUT# 8 ,G
00125IF G=0 THEN 210
00130I2=0
00135Q9=1
00140RESTORE#7
00141INPUT#7,M
00142INPUT#7,S1
00143INPUT#7,G
00144INPUT#7,R9
00145IF M=0 THEN 210
00146MATN=ZER(M)
00147MATB=ZER(M)
00148MATL=ZER(R9,M)
00149MATM=ZER(R9)
00150MATA=ZER(R9)
00151FORI=1TOM
00152INPUT#7,N(I)
00153NEXTI
00154FORI=1TOM
00155INPUT#7,B(I)
00156NEXTI
00157FOR I=1TOR9
00158FORJ=1TOM
00159INPUT#7,L(I,J)
00160NEXTJ
00161NEXTI
00162FORI=1TOR9
00163INPUT#7,M(I)
00164NEXTI
00165FORI=1TOR9
00166INPUT#7,A(I)
00167NEXTI
00180RESTORE#8
00181 INPUT#8,G0
00182 INPUT#8,M
00183 INPUT#8,S2
00185S7=1
00190MAT E=ZER(M,M)
00191 MAT T=ZER(M,M)
00195MAT F=ZER(M)
00200FORI=1TOM
00201FORJ=1TOM
00202INPUT#8,E(I,J)
00203NEXTJ
00204NEXTI
00205FORI=1TOM
00206INPUT#8,F(I)
00207NEXTI
00208GOTO1151
00210PRINT L$
00215PRINT " POSTERIOR DISTRIBUTIONS FOR ANOVA"
00220RESTORE#2
00221 INPUT# 2,M
00225IF M <> 0 THEN 265
00230PRINT " "
00235PRINT "BEFORE YOU CAN DO THE POSTERIOR ANALYSIS YOU MUST ENTER"
00240PRINT "THE SUMMARY STATISTICS FOR YOUR DATA. IF YOU HAVE THESE"
00245PRINT "YOU SHOULD USE MODULE 1 OF THIS MODEL TO ENTER THEM."
00250PRINT "TO CONTINUE TYPE '1'.";
00255GOSUB 9000
00260CHAIN "RSTRT"
00265N0=0
00270S0=0
00275MAT E=ZER(M,M)
00276 MAT T=ZER(M,M)
00280MAT N=ZER(M)
00285MAT S=ZER(M)
00290MAT M=ZER(M)
00295MAT B=ZER(M)
00300MAT F=ZER(M)
00301FORI=1TOM
00302INPUT#2,N(I)
00303NEXTI
00304FORI=1TOM
00305INPUT#2,M(I)
00306NEXTI
00307FORI=1TOM
00308INPUT#2,S(I)
00309NEXTI
00310MAT B=M
00315MAT F=M
00320FOR I=1 TO M
00325S0=S0+N(I)*S(I)*S(I)
00330N0=N0+N(I)
00335NEXT I
00340G=N0-M
00345FOR I=1 TO M
00350S(I)=SQR(S0/(G-2)/N(I))
00355E(I,I)=1/N(I)
00360NEXT I
00365S1=S0
00370S2=S1
00375G0=G
00380PRINT " "
00385GOTO 395
00390PRINT L$
00395PRINT "JOINT DISTRIBUTION OF THE PARAMETERS IS MULTIVARIATE T."
00400PRINT "MARGINAL DISTRIBUTIONS ARE UNIVARIATE T. DF = ";G
00405PRINT "---------------------------------------------------"
00410PRINT " MEAN STANDARD DEVIATION"
00415FOR I=1 TO M
00420IF S7=1 THEN 460
00425IF I2=16 THEN 445
00430: CELL## ######.## ######.###
00435PRINT USING 430,I,M(I),S(I)
00440GOTO 470
00445:## CELL 'C ######.## ######.###
00450 PRINT USING 445,I,MID$(C$,(I-1)*2+1,(I-1)*2+2-((I-1)*2+1)+1),M(I),S(I)
00455GOTO 470
00460:PARAMETER## ######.## ######.###
00465PRINT USING 460,I,M(I),S(I)
00470NEXT I
00475PRINT "--------------------------------------------------"
00480PRINT " OPTIONS"
00485PRINT "1. EXAMINE MARGINAL"
00490IF M>1 THEN 505
00495PRINT "2. RETURN TO ORIGINAL BASIS"
00500GOTO 555
00505PRINT "2. SINGLE LINEAR COMBINATION"
00510PRINT "3. PAIR-WISE COMPARISONS"
00515PRINT "4. EXAMINE CONTOURS"
00520IF Q9=1 THEN 535
00525PRINT "5. MULTIPLE LINEAR COMBINATIONS"
00530GOTO 540
00535PRINT "5. RETURN TO ORIGINAL BASIS"
00540IF Q9=1 THEN 555
00545IF I3 <> 16 THEN 555
00550PRINT "6. ANOVA EFFECTS"
00555PRINT "INPUT THE NUMBER OF THE OPTION YOU WANT OR '0' TO EXIT.";
00560GOSUB 9000
00565IF M>1 THEN 595
00570IF O1 <> 1 THEN 580
00575CHAIN "CMOD86"
00580IF O1=2 THEN 1440
00585PRINT "MUST BE 1 OR 2. RESPECIFY."
00590GOTO 560
00595IF O1 <> 6 THEN 630
00600IF Q9=1 THEN 700
00605IF I3 <> 16 THEN 640
00610SCRATCH#3
00611 PRINT # 3,1
00615IF I9=1 THEN 625
00620CHAIN "CMOD89"
00625CHAIN "CMOD75"
00630IF O1 <> 0 THEN 640
00635CHAIN "RSTRT"
00640IF O1 <> 1 THEN 650
00645CHAIN "CMOD86"
00650IF O1=4 THEN 1180
00655IF O1 <> 2 THEN 665
00660CHAIN "CMOD87"
00665IF O1 <> 3 THEN 675
00670CHAIN "CMOD88"
00675IF Q9 <> 1 THEN 690
00680IF O1 <> 5 THEN 700
00685GOTO 1440
00690IF O1=5 THEN 720
00695IF Q9=1 THEN 710
00700PRINT "MUST BE 1,2,3,4 OR 5. RESPECIFY."
00705GOTO 560
00710PRINT "INPUT 1,2,3,4,5 OR 6. RESPECIFY"
00715CHAIN "RSTRT"
00720IF I2 <> 16 THEN 740
00725IF I1=75 THEN 735
00730CHAIN "CMOD89"
00735CHAIN "CMOD75"
00740PRINT L$
00745PRINT " OPTION 5. MULTIPLE LINEAR COMBINATIONS"
00750PRINT " "
00755PRINT "THIS OPTION ALLOWS YOU TO EXAMINE THE MULTIVARIATE JOINT, "
00760PRINT "CONDITIONAL, AND MARGINAL T DISTRIBUTION OF THE VARIOUS"
00765PRINT "LINEAR COMBINATIONS YOU SPECIFY. THESE LINEAR COMBINATIONS"
00770PRINT "DEFINE NEW PARAMETERS AND MUST BE LINEARLY INDEPENDENT."
00775PRINT "YOU MUST ENTER ";M;" LINEAR COMBINATIONS."
00780R0=M
00785MAT I=ZER(M)
00790MAT Z=ZER(M,M)
00795FOR I=1 TO R0
00800PRINT "ENTER COEFFICIENTS FOR COMBINATION ";I;" SEPARATED BY COMMAS."
00805MAT INPUT I
00810FOR J=1 TO M
00815 Z(I,J)=I(J)
00820NEXT J
00825NEXT I
00830C0=M
00835MAT L=ZER(R0,C0)
00840MAT V=ZER(R0)
00850MAT A=ZER(C0)
00855MAT T=ZER(C0,R0)
00860MAT L=Z
00865GOSUB 9500
00870IF L0 <> 0 THEN 905
00875PRINT "LINEAR COMBINATIONS AREN'T FULL RANK."
00880PRINT "TO REDEFINE THEM TYPE '1'"
00885PRINT "ELSE TYPE '0'.";
00890GOSUB 9000
00895IF O1=1 THEN 780
00900GOTO 390
00905MAT L=Z
00910PRINT L$
00915MAT M=ZER(M)
00920PRINT "YOU MUST NOW SPECIFY HOW EACH PARAMETER IS TO BE TREATED IN THIS"
00925PRINT "ANALYSIS. YOU HAVE THREE OPTIONS: YOU CAN INDICATE IT IS A"
00930PRINT "PARAMETER ABOUT WHICH A PROBABILITY STATEMENT WILL BE MADE, "
00935PRINT "YOU CAN CONDITIONALIZE ON SOME VALUE, OR YOU CAN MARGINALIZE"
00940PRINT "WITH RESPECT TO IT."
00945PRINT " "
00950PRINT "FOR EACH PARAMETER:"
00955PRINT " TO MAKE A STATEMENT ABOUT IT TYPE, '1'."
00960PRINT " TO CONDITIONALIZE ON SOME VALUE OF IT, TYPE '2'."
00965PRINT " TO MARGINALIZE WITH RESPECT TO IT, TYPE '3'."
00970PRINT " "
00975FOR I=1 TO M
00980PRINT "PARAMETER ";I;
00985GOSUB 9000
00990O1=INT(O1)
00995IF O1>3 THEN 1065
01000IF O1<1 THEN 1065
01005M(I)=O1
01010IF O1 <> 2 THEN 1075
01015IF I2=15 THEN 1045
01020IF O1 <> 2 THEN 1075
01025PRINT "MARGINAL DISTRIBUTION FOR PARAMETER ";I;" IS DEFINED BY:"
01030:DF = ###### MEAN = ######.## STD. DEV. = ######.###
01035PRINT USING 1030,G,F(I),S(I)
01040PRINT " "
01045PRINT "ENTER THE CONDITIONING VALUE OF PARAMETER ";I;
01050GOSUB 9000
01055A(I)=O1
01060GOTO 1075
01065PRINT "INPUT MUST BE 1, 2, OR 3."
01070GOTO 980
01075NEXT I
01080R9=M
01085FOR I=1 TO M
01090IF M(I)=1 THEN 1120
01095NEXT I
01100PRINT "YOU MUST MAKE A STATMENT ABOUT AT LEAST ONE PARAMETER."
01105PRINT "TYPE '1' TO CONTINUE.";
01110GOSUB 9000
01115GOTO 910
01120IF I2=15 THEN 1126
01121RESTORE#7
01122INPUT#7,M
01123INPUT#7,S1
01124INPUT#7,G
01125INPUT#7,R9
01126 SCRATCH#7
01127PRINT#7,M
01128PRINT#7,S1
01129PRINT#7,G
01130PRINT#7,R9
01132FORI=1TOM
01133PRINT#7,N(I)
01134NEXTI
01135FOR I=1TOM
01136PRINT#7,B(I)
01138NEXTI
01139FORI=1TOR9
01140FORJ=1TOM
01141PRINT#7,L(I,J)
01142NEXTJ
01143NEXTI
01144FORI=1TOR9
01145PRINT#7,M(I)
01146NEXTI
01147FORI=1TOR9
01148PRINT#7,A(I)
01149NEXTI
01150CHAIN"CMOD70"
01151FORI=1TOM
01152FORJ=1TOM
01155E(I,J)=E(I,J)/(S2/G0)
01160NEXT J
01165NEXT I
01170MAT T=INV(E)
01172 MAT E=T
01175GOTO 1470
01180 MAT T=INV(E)
01182 MAT E=T
01185PRINT " "
01190PRINT " OPTION 4. EXAMINATION OF CONTOURS"
01195S9=0
01200PRINT " "
01205MAT V=ZER(M)
01210PRINT "ENTER THE ";M;" POINTS FOR WHICH YOU WISH TO SEE THE CONTOUR."
01215MAT INPUT V
01220MAT C=ZER(M)
01225FOR I=1 TO M
01230C(I)=V(I)-F(I)
01235NEXT I
01240MAT G=ZER(1,M)
01245MAT T=ZER(1,M)
01250MAT G=TRN(C)
01255MAT T=G*E
01265MAT U=T*C
01270F=U(0)/((S2/G0)*M)
01275PRINT L$
01280PRINT "THE POINT DEFINED BY"
01285FOR I=1 TO M
01290PRINT "PARAMETER ";I;" = ";V(I)
01295NEXT I
01300IF M>1 THEN 1335
01305S9=1
01310G=G0/2
01315M0=F(1)
01320S0=E(1,1)
01325GOSUB 1800
01330GOTO 1375
01335A=M/2
01340B=G0/2
01345J1=0
01350J2=F*A/(B+F*A)
01355M9=G0
01360GOSUB 5220
01365GOSUB 2280
01370G0=M9
01375IF P>0 THEN 1385
01380P=0
01385IF P <= .99 THEN 1400
01390PRINT "IS NOT IN THE 99% HIGHEST DENSITY CONTOUR."
01395GOTO 1410
01400:IS IN THE ##% HIGHEST DENSITY CONTOUR.
01405PRINT USING 1400,P*100
01410PRINT "TO ENTER ANOTHER SET TYPE '1'"
01415PRINT "ELSE '0'";
01420GOSUB 9000
01425IF O1=1 THEN 1205
01430MAT T=INV(E)
01432 MAT E=T
01435GOTO 390
01440IF I3 <> 16 THEN 1460
01445IF I9=1 THEN 1455
01450CHAIN "CMOD89"
01455CHAIN "CMOD75"
01460SCRATCH#8
01461 PRINT # 8 ,0
01465CHAIN "CMOD69"
01470MAT S=ZER(M)
01475MAT M=ZER(M)
01480MAT M=F
01485MAT T=INV(E)
01486 MAT E=T
01490FOR I=1 TO M
01495S(I)=SQR(S2*E(I,I)/(G0-2))
01500NEXT I
01505G=G0
01510RESTORE#3
01511 INPUT# 3,O1
01515IF O1=70 THEN 390
01520GOTO 910
01525REM
01530PRINT L$
01535PRINT "YOU HAVE THE FOLLOWING OPTIONS AVAILABLE FOR EXAMINATION OF"
01540PRINT H$
01545PRINT
01550PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
01555PRINT " 1. HIGHEST DENSITY REGION"
01560PRINT " 2. PROBABILITY LESS THAN SOME VALUE."
01565PRINT " 3. PROBABILITY BETWEEN TWO VALUES."
01570PRINT " 4. EXIT"
01575GOSUB 9000
01580IF O1 <> 4 THEN 1590
01585RETURN
01590IF O1=1 THEN 1615
01595IF O1=2 THEN 1760
01600IF O1=3 THEN 1875
01605PRINT "MUST BE 1,2,3 OR 4. RESPECIFY."
01610GOTO 1575
01615PRINT L$
01620PRINT " HIGHEST DENSITY REGIONS"
01625PRINT " "
01630GOTO 1675
01635PRINT H$
01640PRINT " "
01645:MEAN ######.## DEGREES OF FREEDOM ######
01650PRINT USING 1645,M0,G
01655:STANDARD DEV. ####.### SCALE PARAMETER ########.##
01660PRINT USING 1655,S0,S0*S0*(G-2)
01665PRINT "----------------------------------------------------------------"
01670RETURN
01675GOSUB 1635
01680PRINT "TO EXIT ROUTINE TYPE '0' WHEN ASKED FOR INPUT."
01685PRINT "INPUT P% AS NUMBER FROM 1 THROUGH 99."
01690PRINT "P%=";
01695GOSUB 9000
01700IF O1=0 THEN 1525
01705IF O1>99 THEN 1720
01710IF O1<1 THEN 1720
01715GOTO 1730
01720PRINT "P% MUST BE AT LEAST 1% AND NOT GREATER THAN 99%. REENTER."
01725GOTO 1695
01730P0=O1/200+.5
01735GOSUB 2055
01740: ##.#% HDR = #####.## TO ######.##
01745F0=J2*SQR(S0*S0*(G-2)/G)
01750PRINT USING 1740,O1,M0-F0,M0+F0
01755GOTO 1690
01760PRINT L$
01765PRINT " PROBABILITY LESS THAN SOME VALUE"
01770PRINT " "
01775GOSUB 1635
01780PRINT "TO EXIT ROUTINE TYPE '7777' WHEN ASKED FOR INPUT."
01785PRINT "INPUT VALUE";
01790GOSUB 9000
01795IF O1=7777 THEN 1525
01800X3=O1
01805Y0=ABS(M0-X3)
01810J2=Y0/S0/SQR((G-2)/G)
01815J1=0
01820GOSUB 6000
01825P=P-.5
01830IF X3<M0 THEN 1845
01835P=.5+P
01840GOTO 1850
01845P=.5-P
01850IF S9=0 THEN 1860
01855RETURN
01860: PROBABILITY LESS THAN ######.## = #.##
01865PRINT USING 1860,O1,P
01870GOTO 1785
01875PRINT L$
01880PRINT " PROBABILITY BETWEEN TWO VALUES"
01885PRINT
01890GOSUB 1635
01895PRINT "TO EXIT ROUTINE TYPE '0,0' WHEN ASKED FOR INPUT."
01900PRINT "INPUT TWO VALUES SEPARATED BY A COMMA. SMALLEST FIRST."
01905GOSUB 9050
01910IF O1 <> 0 THEN 1925
01915IF O2 <> 0 THEN 1925
01920GOTO 1525
01925IF O1 <= O2 THEN 1940
01930PRINT "SMALLER VALUE MUST BE ENTERED FIRST. RESPECIFY."
01935GOTO 1905
01940X3=O1
01945X4=O2
01950Y0=ABS(M0-X3)
01955J2=Y0/S0/SQR((G-2)/G)
01960J1=0
01965GOSUB 6000
01970P=P-.5
01975IF X3<M0 THEN 1990
01980P3=.5+P
01985GOTO 1995
01990P3=.5-P
01995Y0=ABS(M0-X4)
02000J2=Y0/S0/SQR((G-2)/G)
02005J1=0
02010GOSUB 6000
02015P4=P-.5
02020IF X4<M0 THEN 2035
02025P4=.5+P4
02030GOTO 2045
02035P4=.5-P4
02040: PROB(######.## < T < ######.##)=#.##
02045PRINT USING 2040,X3,X4,P4-P3
02050GOTO 1900
02055REM %ILE FINDER
02060P5=2
02065P6=2
02070P2=P0
02075GOSUB 2205
02080IF ABS(P-P0)<.0001 THEN 2200
02085IF P>P0 THEN 2120
02090E5=J2
02095P5=P
02100IF P6 <> 2 THEN 2150
02105P2=P2+.001
02110GOSUB 2205
02115GOTO 2080
02120E6=J2
02125P6=P
02130IF P5 <> 2 THEN 2150
02135P2=P2-.001
02140GOSUB 2205
02145GOTO 2080
02150J2=.5*(E6+E5)
02155GOSUB 6000
02160IF ABS(P-P0)<.0001 THEN 2200
02165IF P>P0 THEN 2185
02170P5=P
02175E5=J2
02180GOTO 2150
02185E6=J2
02190P6=P
02195GOTO 2150
02200RETURN
02205P3=P2
02210IF P2 <= .5 THEN 2220
02215P2=1-P2
02220A1=SQR(LOG(1/P2/P2))
02225A2=2.51552+.802853*A1+.010328*A1*A1
02230A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
02235A2=A1-A2
02240J2=SQR(G*EXP(A2*(G-5/6)*A2/(G-2/3+.1/G)/(G-2/3+.1/G))-G)
02245GOSUB 6000
02250IF P3 <= .5 THEN 2265
02255P2=P3
02260GOTO 2275
02265J2=-J2
02270P=1-P
02275RETURN
02280REM BETA CDF
05025 DIM W(16),O(16)
05035IF A+B>85 THEN 5280
05040P=0
05045C6=0
05050IF A>1 THEN 5080
05055C6=A
05060C7=B
05065A=C7
05070B=C6
05075J2=1-J2
05080D0=(J2-J1)*.5
05085D1=(J1+J2)*.5
05090FOR I1=1 TO 16
05095D9=D0*O(I1)+D1
05100IF D9=0 THEN 5115
05105IF D9=1 THEN 5115
05106D9=LOG(D9)*(A-1)+LOG(1-D9)*(B-1)
05108IF D9<-80 THEN 5115
05110P=P+W(I1)*EXP(D9)
05115NEXT I1
05120P=P*F0
05125P=P*D0
05130IF C6=0 THEN 5155
05135A=C6
05140B=C7
05145P=1-P
05150J2=1-J2
05155RETURN
05160FOR I1=1 TO 16
05165READ W(I1),O(I1)
05170NEXT I1
05175DATA 2.71525E-02,-.989401
05180DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
05185DATA .124629,-.755404,.149596,-.617876
05190DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
05195DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
05200DATA .149596,.617876,.124629,.755404
05205DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
05210DATA .989401
05215RETURN
05220G9=A+B
05225GOSUB 5850
05230F0=G0
05235G9=A
05240GOSUB 5850
05245F0=F0-G0
05250G9=B
05255GOSUB 5850
05260F0=F0-G0
05265IF A+B>85 THEN 5275
05270F0=EXP(F0)
05275RETURN
05280W1=(B*J2)**(1/3)
05285W2=(A*(1-J2))**(1/3)
05290GOSUB 5325
05295I1=P
05300W1=(B*J1)**(1/3)
05305W2=(A*(1-J1))**(1/3)
05310GOSUB 5325
05315P=I1-P
05320RETURN
05325Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A)
05330GOSUB 8000
05335RETURN
05850REM LOG GAMMA
05853REM INPUT G9
05860G5=G9
05863IF G9 <= 1.E+30 THEN 5872
05866G0=1.E+38
05869RETURN
05872IF G9>1.E-09 THEN 5881
05875G0=0
05878RETURN
05881IF G9<1.E+10 THEN 5890
05884G0=G9*(LOG(G9)-1)
05887RETURN
05890G6=1
05893IF 18<G5 THEN 5905
05896G6=G6*G5
05899G5=G5+1
05902GOTO 5893
05905R8=1/G5/G5
05908G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
05911C1=8.33333E-02
05914C2=2.77778E-03
05917C3=7.93651E-04
05920C4=5.95238E-04
05923G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
05926RETURN
06000REM T CDF
06003L0=LOG(1+J2*J2/G)
06004IF L0>.000001 THEN 6006
06005L0=0
06006P=.5
06007IF J2=0 THEN 6035
06008IF J2<6 THEN 6013
06009P=1
06010GOTO 6035
06011IF J2>12 THEN 6009
06013Y3=J2
06015Y3=(G-2/3+.1/G)*SQR(L0/(G-5/6))
06016GOSUB 8000
06017GOTO 6035
06035RETURN
08000REM NORMAL CDF
08050Y4=ABS(Y3)
08060X1=X
08070X=Y3
08080T=1/(1+.231642*Y4)
08090D=.398942*EXP(-X*X/2)
08100C1=1.33027
08110C2=1.82126
08120C3=1.78148
08130C4=.356564
08140C5=.319382
08150P=1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5)
08160IF X >= 0 THEN 8180
08170P=1-P
08180X=X1
08190RETURN
09000REM
09005INPUT O1
09010IF O1=-9999 THEN 9020
09015RETURN
09020CHAIN "RSTRT"
09050INPUT O1,O2
09055IF O2=-9999 THEN 9020
09060GOTO 9010
09500L0=0
09505L1=1
09510FOR L3=1 TO R0
09520L4=L3
09525L5=L3+1
09530FOR L6=L5 TO R0
09535IF ABS(L(L4,L3)) >= ABS(L(L6,L3)) THEN 9545
09540L4=L6
09545NEXT L6
09550IF ABS(L(L4,L3))>10 ** (-7) THEN 9560
09552IF R0<C0 THEN 9655
09555RETURN
09560FOR L6=L3 TO C0
09565L9=L(L3,L6)
09570L(L3,L6)=L(L4,L6)
09575L(L4,L6)=L9
09580NEXT L6
09585L9=L(L3,L3)
09590L(L3,L3)=1
09595FOR L6=L5 TO C0
09600L(L3,L6)=L(L3,L6)/L9
09605NEXT L6
09610FOR L6=1 TO R0
09615IF L6=L3 THEN 9640
09620L9=L(L6,L3)
09625FOR L8=L3 TO C0
09630L(L6,L8)=L(L6,L8)-L9*L(L3,L8)
09635NEXT L8
09640NEXT L6
09645NEXT L3
09650IF R0 >= C0 THEN 9675
09655FOR L1=1 TO R0
09660IF L(R0,L1) <> 0 THEN 9675
09665NEXT L1
09670RETURN
09675L0=1
09680RETURN
09999END