BEGIN	COMMENT     Title : Multiple Linear Regression Analysis,
		   Author : Marten van Gelderen,
		  Version : 5H(246), Simula 67, DEC-10, Date 800216,
		Copyright : Foundation Mathematical Centre, Amsterdam;

	COMMENT ************************************* external procedures ********************************;

	EXTERNAL REAL PROCEDURE Maxreal;
	EXTERNAL BOOLEAN PROCEDURE Dotypeout;
	EXTERNAL REF (Infile) PROCEDURE Findinfile;
	EXTERNAL REF (Printfile) PROCEDURE Findprintfile;
	EXTERNAL REF (Directfile) PROCEDURE Finddirectfile;
	EXTERNAL LONG REAL PROCEDURE Scanreal, Lmin, Lmax;
	EXTERNAL INTEGER PROCEDURE Maxint, Scanint, Imin, Imax;
	EXTERNAL TEXT PROCEDURE Upcase, Compress, Rest, Scanto, Filename;

	INTEGER scalelimit, optionlimit; scalelimit:= 3; optionlimit:= 10; Lowten ('#'); Simset

BEGIN	INTEGER ARRAY hashprime [1 : scalelimit]; BOOLEAN ARRAY optionwanted [1 : optionlimit];

	COMMENT ************************************* error messages *************************************;

	BOOLEAN erroneous;

	PROCEDURE overflow (string); NAME string; TEXT string;
	INSPECT printstream DO
	BEGIN	Eject (Line + 5); Outtext (string); Outtext (" overflow");
		IF tablescale LT scalelimit THEN
		BEGIN	Outtext ("    (message)"); Outimage;
			Eject (Line + 5); Outtext ("Rerun of job :"); Outint (jobnumber, 3); Outimage;
			Eject (1); tablescale:= tablescale + 1; stacklimit:= stacklimit * 2; GOTO rerun;
		END; Outimage; erroneous:= TRUE; GOTO endjob;
	END *** overflow ***;

	PROCEDURE errortext (errnum); INTEGER errnum;
	INSPECT printstream DO
	BEGIN	Eject (Line + 1); Outtext ("Error : "); INSPECT errorstream DO
		BEGIN	Locate (errnum - errnum // 10 + 11); Inimage;
			printstream . Outtext (Image . Strip); printstream . Outimage;
		END OTHERWISE Outint (errnum, 2); erroneous:= TRUE;
	END *** errortext ***;

	PROCEDURE error (errnum); INTEGER errnum;
	INSPECT printstream DO
	BEGIN	REF (atomcell) atom; errortext (errnum);
		FOR atom:- lastname, lastnumber, lastsymbol DO
		BEGIN Setpos (Pos + 13); printitem (atom); END; Outimage;
	END *** error ***;

	PROCEDURE error1 (errnum, term); INTEGER errnum, term;
	INSPECT printstream DO
	BEGIN	errortext (errnum); Outint (term, 14); Outimage; END;

	PROCEDURE error2 (errnum, wrong, case); INTEGER errnum; REAL wrong, case;
	INSPECT printstream DO
	BEGIN	errortext (errnum); Outfix (wrong, 5, 20); Outfix (case, 5, 20); Outimage; GOTO endjob; END;

	COMMENT ************************************* file handling **************************************;

	INTEGER imagelength, errorlength;

	REF (Infile) inputstream;	REF (Printfile) printstream;
	REF (Infile) datastream;	REF (Printfile) outputstream;	   REF (Directfile) errorstream;

	TEXT PROCEDURE enterspec (stream); NAME stream; TEXT stream;
	BEGIN Outtext (stream); Breakoutimage; Inimage; enterspec:- Sysin . Image . Strip; END;

	REF (Infile) PROCEDURE getinpfile (filespec, stream); NAME stream; TEXT filespec, stream;
	IF filespec == NOTEXT THEN getinpfile:- NONE ELSE INSPECT Findinfile (filespec) DO
	getinpfile:- THIS Infile OTHERWISE getinpfile:- getinpfile (enterspec (stream), stream);

	REF (Printfile) PROCEDURE getprifile (filespec, stream); NAME stream; TEXT filespec, stream;
	IF filespec == NOTEXT THEN getprifile:- NONE ELSE INSPECT Findprintfile (filespec) DO
	getprifile:- THIS Printfile OTHERWISE getprifile:- getprifile (enterspec (stream), stream);

	PROCEDURE enterfilespecs;
	BEGIN	Eject (Line + 1); Outtext ("Enter file specifications"); Outimage;
		inputstream  :- getinpfile (enterspec ("Inputstream  : "), "Inputstream  ? ");
		printstream  :- getprifile (enterspec ("Printstream  : "), "Printstream  ? ");
		datastream   :- getinpfile (enterspec ("Datastream   : "), "Datastream   ? ");
		outputstream :- getprifile (enterspec ("Outputstream : "), "Outputstream ? ");
		errorstream  :- Finddirectfile ("Hlp: Mulerr.hlp", FALSE);
		Eject (Line + 1); Outtext ("Executing ..."); Outimage; Outimage;
	END *** enterfilespecs ***;

	PROCEDURE readinpspec (filespec); TEXT filespec;
	IF filespec . More AND filespec . Strip NE Filename (inputstream) . Strip THEN
	BEGIN	IF inputstream =/= Sysin THEN inputstream . Close;
		inputstream:- getinpfile (filespec . Strip, "Inputstream  ? "); runable:= FALSE;
		INSPECT inputstream DO Open (Blanks (imagelength)) OTHERWISE inputstream:- Sysin;
	END *** readinpspec ***;

	PROCEDURE readprispec (filespec); TEXT filespec; IF filespec . More THEN
	BEGIN	IF printstream =/= Sysout THEN printstream . Close;
		printstream:- getprifile (filespec . Strip, "Printstream  ? ");
		INSPECT printstream DO Open (Blanks (imagelength)) OTHERWISE printstream:- Sysout;
	END *** readprispec ***;

	PROCEDURE readdatspec (filespec); TEXT filespec; IF filespec . More THEN
	BEGIN	INSPECT datastream DO Close; datastream:- getinpfile (filespec . Strip, "Datastream   ? ");
		INSPECT datastream DO Open (Blanks (imagelength)); datarecord:= 0; datafileread:= FALSE;
	END *** readdatspec ***;

	PROCEDURE readoutspec (filespec); TEXT filespec; IF filespec . More THEN
	BEGIN	INSPECT outputstream DO Close; outputstream:- getprifile (filespec . Strip, "Outputstream ? ");
		INSPECT outputstream DO Open (Blanks (imagelength));
	END *** readoutspec ***;

	PROCEDURE readfilespecs;
	BEGIN	runable:= TRUE; IF lastchar EQ '/' AND lastline . More THEN
		BEGIN	TEXT scanline;
			lastline:- Copy (Rest (lastline)); scanline:- Scanto (lastline, '=');
			readprispec (Scanto (scanline, ';')); readoutspec (Scanto (scanline, ';'));
			readinpspec (Scanto (lastline, ';')); readdatspec (Scanto (lastline, ';'));
			IF NOT runable THEN readline (inputstream);
		END;
	END *** readfilespecs ***;

	PROCEDURE openfiles;
	BEGIN	INSPECT inputstream  DO Open (Blanks (imagelength)) OTHERWISE inputstream:- Sysin;
		INSPECT printstream  DO Open (Blanks (imagelength)) OTHERWISE printstream:- Sysout;
		INSPECT datastream   DO Open (Blanks (imagelength)); Sysout . Linesperpage (-1);
		INSPECT outputstream DO Open (Blanks (imagelength)); Sysout . Image:- Blanks (imagelength);
		INSPECT errorstream  DO Open (Blanks (errorlength));
	END *** openfiles ***;

	PROCEDURE closefiles;
	BEGIN	INSPECT inputstream DO Close;	INSPECT printstream  DO Close;
		INSPECT datastream  DO Close;	INSPECT outputstream DO Close;	INSPECT errorstream DO Close;
	END *** closefiles ***;

	COMMENT ************************************* text & data storage ********************************;

	INTEGER bufferlength, resetindex;

	Link CLASS textbuffer; BEGIN TEXT line; END;

	CLASS textstorage;
	BEGIN	INTEGER size;
		REF (textbuffer) pointer; REF (Head) base;

		PROCEDURE reset;
		BEGIN pointer:- base . First; endtext:= FALSE; endline:= TRUE; END;

		PROCEDURE store (line); TEXT line;
		BEGIN	pointer:- NEW textbuffer;
			pointer . Into (base); pointer . line:- Copy (line);
			pointer . line . Setpos (line . Pos - (IF endline THEN 0 ELSE 1));
			size:= size + (IF endline THEN 0 ELSE 1);
		END *** store ***;

		TEXT PROCEDURE nextline;
		BEGIN	nextline:- lastline:- pointer . line;
			pointer:- pointer . Suc; endtext:= pointer == NONE;
		END *** nextline ***;

		size:= 0; pointer:- NONE; base:- NEW Head;
	END *** textstorage ***;

	Link CLASS databuffer;
	BEGIN LONG REAL ARRAY buffer [1 : bufferlength]; END;

	CLASS datastorage;
	BEGIN	INTEGER index, size;
		REF (databuffer) pointer; REF (Head) base;

		PROCEDURE store (number); LONG REAL number;
		BEGIN	IF index GT bufferlength THEN
			BEGIN	pointer:- NEW databuffer;
				pointer . Into (base); index:= 1;
			END; pointer . buffer [index]:= number;
			index:= index + 1; size:= size + 1;
		END *** store ***;

		PROCEDURE bufferror; INSPECT printstream DO
		BEGIN errortext (3); Outimage; GOTO endjob; END;

		size:= 0; index:= bufferlength + 1;
		pointer:- NONE; base:- NEW Head;
	END *** datastorage ***;

	datastorage CLASS inputstorage;
	BEGIN	PROCEDURE reset;
		BEGIN pointer:- base . First; index:= 0; END;

		LONG REAL PROCEDURE nextnumber;
		BEGIN	IF index LT bufferlength THEN index:= index + 1
			ELSE BEGIN pointer:- pointer . Suc; index:= 1; END;
			IF pointer == NONE THEN bufferror;
			nextnumber:= pointer . buffer [index];
		END *** nextnumber ***;

		PROCEDURE skipnumbers (count); INTEGER count;
		BEGIN	index:= index + count; WHILE index GE bufferlength DO
			BEGIN	index:= index - bufferlength; pointer:- pointer . Suc;
				IF pointer == NONE THEN bufferror;
			END;
		END *** skipnumbers ***;
	END *** inputstorage ***;

	datastorage CLASS rightstorage;
	BEGIN	INTEGER address;

		LONG REAL PROCEDURE lastnumber;
		BEGIN	IF index GT 1 THEN index:= index - 1 ELSE
			BEGIN pointer:- pointer . Pred; index:= bufferlength; END;
			IF pointer == NONE THEN bufferror;
			lastnumber:= pointer . buffer [index];
		END *** lastnumber ***;
	END *** rightstorage ***;

	rightstorage CLASS leftstorage;
	BEGIN	INTEGER repeats;

		PROCEDURE fixrepeats;
		BEGIN	repeats:= size - repeats;
			store (repeats); repeats:= size;
		END *** fixrepeats ***;

		repeats:= 0;
	END *** leftstorage ***;

	rightstorage CLASS designstorage;
	BEGIN	LONG REAL sum, squares, min, max;

		PROCEDURE reset;
		BEGIN pointer:- base . Last; index:= resetindex; END;

		sum:= squares:= 0; min:= Maxreal; max:= - Maxreal;
	END *** designstorage ***;

	COMMENT ************************************* atom cells *****************************************;

	CLASS atomcell;
	BEGIN TEXT item, entry; INTEGER index, type; END;

	atomcell CLASS identifier;;	atomcell CLASS number;;
	atomcell CLASS functionname;;	atomcell CLASS optionname;;

	REF (atomcell)	lastatom, lastname, lastnumber, lastsymbol,
			addsign, subsign, mulsign, divsign, idisign, ttpsign, withsign,
			equal, comma, lparent, rparent, lbracket, rbracket, smaller, greater, semicolon;

	COMMENT ************************************* item & atom hashing ********************************;

	INTEGER maxrank; BOOLEAN declared; REF (hashing) atomhash;

	CLASS hashing (hashprime); INTEGER hashprime;
	BEGIN	INTEGER twinprime, key, maxkey, vacant, hashlimit, index, shift; TEXT entry;

		REF (atomcell) ARRAY hashtable [1 : hashprime];

		PROCEDURE hasher;
		BEGIN	key:= maxkey - maxrank; WHILE entry . More DO
			key:= Mod (entry . Pos + key * Rank (entry . Getchar), maxkey);
		END *** hasher ***;

		BOOLEAN PROCEDURE present;
		BEGIN	entry:- Upcase (Compress (Copy (lastitem), ' ')); hasher;
			index:= Mod (key, hashprime) + 1; shift:= Mod (key, twinprime) + 1;
			declared:= FALSE; WHILE NOT declared AND hashtable [index] =/= NONE DO
			IF hashtable [index] . entry EQ entry THEN
			declared:= TRUE ELSE index:= Mod (index + shift, hashprime) + 1;
			present:= declared;
		END *** present ***;

		REF (atomcell) PROCEDURE insert (newcell); REF (atomcell) newcell;
		BEGIN	IF vacant LT hashlimit THEN overflow ("Hash table"); vacant:= vacant - 1;
			newcell . item:- lastitem; newcell . entry:- entry; newcell . index:= index;
			insert:- hashtable [index]:- newcell; newcell . entry . Setpos (1);
		END *** insert ***;

		REF (atomcell) PROCEDURE retrieve; retrieve:- hashtable [index];

		REF (atomcell) PROCEDURE reserve (newcell, item, type);
		VALUE item; REF (atomcell) newcell; TEXT item; INTEGER type;
		BEGIN	REF (atomcell) atomtype; lastitem:- item;
			atomtype:- reserve:- IF present THEN retrieve ELSE insert (newcell);
			atomtype . type:= type;
		END *** reserve ***;

		twinprime:= hashprime - 2; maxkey:= Maxint // maxrank;
		vacant:= hashprime; hashlimit:= Entier (0.2 * hashprime);
		FOR index:= 1 STEP 1 UNTIL hashprime DO hashtable [index]:- NONE;

		addsign   :- reserve (NEW atomcell, "+", 0);
		subsign   :- reserve (NEW atomcell, "-", 0);
		mulsign   :- reserve (NEW atomcell, "*", 0);
		divsign   :- reserve (NEW atomcell, "/", 0);
		idisign   :- reserve (NEW atomcell, ":", 0);
		ttpsign   :- reserve (NEW atomcell, "^", 0);
		withsign  :- reserve (NEW atomcell, "&", 0);
		equal     :- reserve (NEW atomcell, "=", 0);
		comma     :- reserve (NEW atomcell, ",", 0);
		lparent   :- reserve (NEW atomcell, "(", 0);
		rparent   :- reserve (NEW atomcell, ")", 0);
		lbracket  :- reserve (NEW atomcell, "[", 0);
		rbracket  :- reserve (NEW atomcell, "]", 0);
		smaller   :- reserve (NEW atomcell, "<", 0);
		greater   :- reserve (NEW atomcell, ">", 0);
		semicolon :- reserve (NEW atomcell, ";", 0);
		reserve (NEW functionname, "Abs", 1) . index:= 1;
		reserve (NEW functionname, "Sign", 2) . index:= 1;
		reserve (NEW functionname, "Sqrt", 3) . index:= 1;
		reserve (NEW functionname, "Sin", 4) . index:= 1;
		reserve (NEW functionname, "Cos", 5) . index:= 1;
		reserve (NEW functionname, "Tan", 6) . index:= 1;
		reserve (NEW functionname, "Ln", 7) . index:= 1;
		reserve (NEW functionname, "Log", 8) . index:= 1;
		reserve (NEW functionname, "Exp", 9) . index:= 1;
		reserve (NEW functionname, "Entier", 10) . index:= 1;
		reserve (NEW functionname, "Round", 11) . index:= 1;
		reserve (NEW functionname, "Mod", 12) . index:= 2;
		reserve (NEW functionname, "Min", 13) . index:= 2;
		reserve (NEW functionname, "Max", 14) . index:= 2;
		reserve (NEW functionname, "Arcsin", 15) . index:= 1;
		reserve (NEW functionname, "Arccos", 16) . index:= 1;
		reserve (NEW functionname, "Arctan", 17) . index:= 1;
		reserve (NEW functionname, "Sinh", 18) . index:= 1;
		reserve (NEW functionname, "Cosh", 19) . index:= 1;
		reserve (NEW functionname, "Tanh", 20) . index:= 1;
		reserve (NEW functionname, "Indicator", 21) . index:= 3;
		reserve (NEW optionname, "Transformed data matrix", 1);
		reserve (NEW optionname, "Correlation matrix", 2);
		reserve (NEW optionname, "Residual analysis", 3);
		reserve (NEW optionname, "No regression analysis", 4);
		reserve (NEW optionname, "Process submodels", 5);
		reserve (NEW optionname, "Print input data", 6);
		reserve (NEW optionname, "No input data rewind", 7);
		reserve (NEW optionname, "Save original model", 8);
		reserve (NEW optionname, "Test reduced model", 9);
		reserve (NEW optionname, "Missing values", 10);
	END *** hashing ***;

	COMMENT ************************************* text & data handling *******************************;

	BOOLEAN endtext, endline, dataread, dataskipped, datafileread;
	INTEGER datarecord; TEXT lastline, lastitem; CHARACTER lastchar;
	REF (textstorage) currenttext, modeltext, inputtext, optiontext; REF (inputstorage) inputdata;

	PROCEDURE readtext (sourcetext);
	NAME sourcetext; REF (textstorage) sourcetext;
	BEGIN	sourcetext:- NEW textstorage; INSPECT sourcetext DO
		BEGIN	WHILE NOT endtext DO
			BEGIN store (lastline); readline (inputstream); END;
			IF size EQ 0 THEN sourcetext:- NONE;
		END;
	END *** readtext ***;

	PROCEDURE readline (stream); REF (Infile) stream;
	INSPECT stream DO
	BEGIN	IF stream == Sysin THEN BEGIN Outchar ('*'); Breakoutimage; END; Inimage;
		lastline:- Image . Strip; WHILE nextchar LT ' ' AND NOT endline DO;
		endtext:= lastchar EQ '"' OR Endfile;
	END *** readline ***;

	PROCEDURE readdata (stream); REF (Infile) stream;
	BEGIN	INTEGER start; LONG REAL nextreal; erroneous:= FALSE;
		inputdata:- NEW inputstorage; INSPECT inputdata DO
		BEGIN	WHILE NOT endtext DO
			BEGIN	WHILE NOT endline DO
				BEGIN	IF IF integral (lastchar) THEN TRUE ELSE numerical (lastchar) THEN
					BEGIN	start:= lastline . Pos - 1; lastline . Setpos (start);
						nextreal:= Scanreal (lastline); IF start EQ lastline . Pos
						THEN error1 (4, size + 1) ELSE store (nextreal);
						lastline . Setpos (lastline . Pos + 1);
					END; WHILE nextchar LT ' ' AND NOT endline DO;
				END; readline (stream);
			END; dataskipped:= erroneous;
			dataread:= size GT 0 AND NOT dataskipped;
			IF dataread THEN BEGIN store (size); reset; END;
		END; INSPECT datastream DO IF Endfile AND datafileread THEN
		BEGIN Close; datastream:- NONE; END ELSE datafileread:= stream == datastream;
	END *** readdata ***;

	PROCEDURE readdatafile;
	INSPECT datastream DO
	BEGIN	readline (datastream);
		WHILE (endtext OR endline) AND NOT Endfile DO
		BEGIN	IF endtext THEN
			BEGIN	datarecord:= datarecord + 1; endtext:= FALSE;
				IF readkeyword EQ "EO" THEN readline (datastream);
			END; IF endline THEN readline (datastream);
		END; readdata (datastream); datarecord:= datarecord + 1;
	END *** readdatafile ***;

	TEXT PROCEDURE readkeyword;
	BEGIN	readkeyword:- Upcase (Copy (nextitem (2)));
		WHILE nextchar NE '"' AND NOT endline DO;
		WHILE nextchar LT ' ' AND NOT endline DO;
	END *** readkeyword ***;

	PROCEDURE readjob;
	BEGIN	TEXT keyword; readline (inputstream);
		runable:= FALSE; WHILE NOT runable DO
		BEGIN	echotext; keyword:- readkeyword; endtext:= FALSE;
			IF keyword EQ "RU" THEN readfilespecs ELSE
			IF keyword EQ "MO" THEN readtext (modeltext) ELSE
			IF keyword EQ "IN" THEN readtext (inputtext) ELSE
			IF keyword EQ "OP" THEN readtext (optiontext) ELSE
			IF keyword EQ "DA" THEN readdata (inputstream) ELSE
			IF keyword EQ "HE" THEN printhelp ELSE
			IF keyword EQ "EX" THEN GOTO exit;
		END; jobnumber:= jobnumber + 1;
	END *** readjob ***;

	PROCEDURE printtext;
	INSPECT printstream DO INSPECT currenttext DO
	BEGIN	reset; Eject (Line + 5);
		WHILE NOT endtext DO BEGIN Outtext (nextline); Outimage; END; reset;
	END *** printtext ***;

	PROCEDURE printdata;
	INSPECT printstream DO INSPECT inputdata DO
	BEGIN	INTEGER i; Eject (Line + 5); Outtext ("""Data"""); Outimage; Eject (Line + 1);
		reset; FOR i:= size - 1 STEP - 1 UNTIL 1 DO Outfix (nextnumber, 3, 12); Outimage;
	END *** printdata ***;

	PROCEDURE printitem (atom); REF (atomcell) atom;
	INSPECT printstream DO INSPECT atom DO Outtext (item . Sub (1, Imin (item . Length, 8)));

	PROCEDURE printhelp;
	BEGIN	TEXT filespec; filespec:- Copy ("Hlp: Mulreg.hlp");
		IF lastchar EQ '/' THEN filespec . Sub (9, 3):= nextitem (3);
		INSPECT printstream DO INSPECT Findinfile (filespec) DO
		BEGIN	Open (Blanks (imagelength)); Inimage; WHILE NOT Endfile DO
			BEGIN Outtext (Image . Strip); Outimage; Inimage; END; Close;
			IF printstream =/= Sysout THEN Eject (1); Dotypeout (Sysout);
			INSPECT Sysout DO BEGIN Outtext ("End of helptext"); Outimage; END;
		END OTHERWISE INSPECT Sysout DO
		BEGIN Outtext ("No help information available"); Outimage; END;
		readline (inputstream);
	END *** printhelp ***;

	INTEGER PROCEDURE field (w); INTEGER w;
	field:= IF w EQ 0 THEN 1 ELSE Entier (Ln (Abs (w)) / 2.3) + (IF w LT 0 THEN 2 ELSE 1);

	PROCEDURE echotext;
	INSPECT printstream DO
	WHILE NOT endtext OR inputstream . Endfile DO
	BEGIN	WHILE NOT endtext DO
		BEGIN	IF inputstream == Sysin AND printstream == Sysout THEN
			BEGIN	IF lastline =/= NOTEXT THEN Outchar ('?')
				ELSE Outtext ("For help type: ""Help""");
			END; Outtext (lastline); Outimage; readline (inputstream);
		END; INSPECT inputstream DO IF Endfile THEN
		BEGIN Close; inputstream:- Sysin; readline (inputstream); END;
	END *** echotext ***;

	CHARACTER PROCEDURE nextchar;
	BEGIN	endline:= NOT lastline . More; lastchar:= ' ';
		WHILE lastchar EQ ' ' AND lastline . More DO lastchar:= lastline . Getchar;
		nextchar:= lastchar;
	END *** nextchar ***;

	BOOLEAN PROCEDURE integral (char); CHARACTER char;
	integral:= Digit (char) OR char EQ '+' OR char EQ '-';

	BOOLEAN PROCEDURE numerical (char); CHARACTER char;
	numerical:= char EQ '.' OR char EQ '#';

	TEXT PROCEDURE nextitem (length); INTEGER length;
	BEGIN	WHILE nextchar LT ' ' AND NOT endline DO; nextitem:-
		lastline . Sub (lastline . Pos - 1, Imin (length, lastline . Length - lastline . Pos + 2));
	END *** nextitem ***;

	REF (atomcell) PROCEDURE nextatom;
	INSPECT currenttext DO INSPECT atomhash DO
	BEGIN	WHILE endline AND NOT endtext DO BEGIN nextline; nextchar; END;
		IF endline AND endtext THEN BEGIN nextatom:- lastatom; END ELSE
		BEGIN	INTEGER start, width; start:= lastline . Pos - 1;
			IF Letter (lastchar) THEN
			BEGIN	WHILE Letter (nextchar) OR Digit (lastchar) DO;
				width:= lastline . Pos - start - 1;
				IF endline THEN width:= width + 1;
				lastitem:- lastline . Sub (start, width) . Strip;
				nextatom:- lastatom:- lastname:-
				IF present THEN retrieve ELSE insert (NEW identifier);
			END ELSE
			IF IF Digit (lastchar) THEN TRUE ELSE numerical (lastchar) THEN
			BEGIN	IF Digit (lastchar) THEN WHILE Digit (nextchar) DO;
				IF lastchar EQ '.' THEN
				BEGIN	IF NOT Digit (nextchar) THEN error (5)
					ELSE WHILE Digit (nextchar) DO;
				END;
				IF lastchar EQ '#' THEN
				BEGIN	IF NOT integral (nextchar) THEN error (6)
					ELSE WHILE Digit (nextchar) DO;
				END;
				width:= lastline . Pos - start - 1;
				IF endline THEN width:= width + 1;
				lastitem:- lastline . Sub (start, width) . Strip;
				nextatom:- lastatom:- lastnumber:-
				IF present THEN retrieve ELSE insert (NEW number);
			END ELSE
			BEGIN	nextchar; lastitem:- lastline . Sub (start, 1);
				nextatom:- lastatom:- lastsymbol:-
				IF present THEN retrieve ELSE insert (NEW atomcell);
			END; WHILE lastchar LT ' ' AND NOT endline DO nextchar;
		END;
	END *** nextatom ***;

	COMMENT ************************************* system variables ***********************************;

	INTEGER stk, add, sub, mul, div, idi, ttp, neg, eql, rep, dep, res, tst, ini,
		frs, fls, smc, tcv, tav, rav, sav, lav, rtn, jmp, fix, skp, fun, savedf, jobnumber,
		termcounter, leftcounter, rightcounter, designcounter, weightcounter, programcounter,
		tablescale, programscale, stacklimit, tablelimit, programspace, stackspace, stacklength;
	LONG REAL precision, Ln10, maxobs, maxexp, savess, savems; BOOLEAN runable, savemodel, savefit;

	PROCEDURE initializesystem;
	BEGIN	Outtext ("Multiple Linear Regression Analysis"); Outimage;
		hashprime [1]:= 271; hashprime [2]:= 523; hashprime [3]:= 1033;
		jobnumber:= datarecord:= 0; tablescale:= 1; programscale:= 1; stacklimit:= 32;
		bufferlength:= 128; imagelength:= 132; errorlength:= 75; maxrank:= 131;
		precision:= &&-18; Ln10:= 2.30258509299404572&&0; maxobs:= &&19; maxexp:= 88;
		stk:= tcv:= 1; add:= tav:= 2; sub:= rav:= 3; mul:= sav:= 4; div:= lav:= 5;
		idi:= rtn:= 6; ttp:= jmp:= 7; neg:= fix:= 8; eql:= skp:= 9; rep:= fun:= 10;
		dep:= 11; res:= 12; tst:= 13; ini:= 14; frs:= 15; fls:= 16; smc:= 17;
		modeltext:- inputtext:- optiontext:- NONE; inputdata:- NONE;
		dataread:= dataskipped:= datafileread:= savemodel:= FALSE;
	END *** initializesystem ***;

	PROCEDURE initializerun;
	BEGIN	lastatom:- lastname:- lastnumber:- lastsymbol:- NONE;
		programcounter:= 1; leftcounter:= rightcounter:= 0; erroneous:= FALSE;
		tablelimit:= hashprime [tablescale]; programspace:= tablelimit * programscale;
		atomhash:- NEW hashing (tablelimit); stacklength:= stackspace:= tablelimit + stacklimit;
	END *** initializerun ***;

	initializesystem; enterfilespecs; openfiles;

  job : readjob;  rerun : initializerun;

BEGIN	INTEGER ARRAY program [1 : programspace]; LONG REAL ARRAY stack [1 : stackspace];

	PROCEDURE load (instruction, address); INTEGER instruction, address;
	IF programcounter LE programspace THEN
	BEGIN	program [programcounter]:= address * 100 + instruction;
		programcounter:= programcounter + 1;
	END ELSE overflow ("Program space");

	COMMENT ************************************* model statement ************************************;

	INTEGER modelinstructions, repeatinstruction, idtype;

	PROCEDURE modelstatement;
	INSPECT modeltext DO
	BEGIN	load (tst, 0); load (fls, 0); lefthandpart; load (dep, 0);
		IF lastatom == equal THEN nextatom ELSE error (11);
		IF lastatom == addsign THEN nextatom; righthandpart; load (rtn, 1);
	END OTHERWISE error (10);

	PROCEDURE lefthandpart;
	BEGIN	idtype:= 2; expression; load (rep, 0); load (rtn, 2);
		repeatinstruction:= programcounter - 1; load (smc, 0);
		idtype:= 1; IF lastatom == withsign THEN
		BEGIN nextatom; expression; END ELSE load (tcv, 1);
	END *** lefthandpart ***;

	PROCEDURE righthandpart;
	BEGIN	term; load (res, 0); termcounter:= 1;
		WHILE lastatom == addsign DO
		BEGIN nextatom; term; load (res, 0); termcounter:= termcounter + 1; END;
	END *** righthandpart ***;

	PROCEDURE expression;
	BEGIN	IF lastatom == subsign THEN
		BEGIN nextatom; term; load (neg, 0); END ELSE
		BEGIN IF lastatom == addsign THEN nextatom; term; END;
		nextterm;
	END *** expression ***;

	PROCEDURE nextterm;
	IF lastatom == addsign THEN
	BEGIN	load (stk, 0); nextatom; term; load (add, 0);
		nextterm;
	END ELSE
	IF lastatom == subsign THEN
	BEGIN	load (stk, 0); nextatom; term; load (sub, 0);
		nextterm;
	END *** nextterm ***;

	PROCEDURE term; BEGIN factor; nextfactor; END;

	PROCEDURE nextfactor;
	IF lastatom == mulsign THEN
	BEGIN	load (stk, 0); nextatom; factor; load (mul, 0);
		nextfactor;
	END ELSE
	IF lastatom == divsign THEN
	BEGIN	load (stk, 0); nextatom; factor; load (div, 0);
		nextfactor;
	END ELSE
	IF lastatom == idisign THEN
	BEGIN	load (stk, 0); nextatom; factor; load (idi, 0);
		nextfactor;
	END *** nextfactor ***;

	PROCEDURE factor; BEGIN primary; nextprimary; END;

	PROCEDURE nextprimary;
	WHILE lastatom == ttpsign DO
	BEGIN load (stk, 0); nextatom; primary; load (ttp, 0); END;

	PROCEDURE primary;
	IF lastatom == lparent THEN
	BEGIN	nextatom; expression;
		IF lastatom == rparent THEN nextatom ELSE error (12);
	END ELSE INSPECT lastatom
	WHEN   identifier DO arithname
	WHEN       number DO arithvalue
	WHEN functionname DO functiondesignator (lastatom)
	WHEN   optionname DO BEGIN error (13); nextatom; END
	OTHERWISE BEGIN error (14); nextatom; END;

	PROCEDURE arithname;
	INSPECT lastatom DO IF idtype EQ 6 THEN
	BEGIN	IF IF declared THEN type EQ 6 ELSE FALSE THEN
		arithvalue ELSE BEGIN error (15); nextatom; END;
	END ELSE BEGIN arithvalue; type:= idtype; END;

	PROCEDURE arithvalue;
	INSPECT lastatom DO BEGIN load (tav, index); nextatom; END;

	PROCEDURE functiondesignator (atom); REF (atomcell) atom;
	BEGIN	INTEGER params; params:= 0;
		IF nextatom == lparent THEN
		BEGIN	nextatom; parameterlist (params);
			IF lastatom == rparent THEN nextatom ELSE error (16);
		END; INSPECT atom DO
		IF params EQ index THEN load (fun, type) ELSE error (17);
	END *** functiondesignator ***;

	PROCEDURE parameterlist (params); NAME params; INTEGER params;
	BEGIN	expression; params:= 1;
		WHILE lastatom == comma DO
		BEGIN load (stk, 0); nextatom; expression; params:= params + 1; END;
	END *** parameterlist ***;

	COMMENT ************************************* input statement ************************************;

	INTEGER inputinstructions; BOOLEAN three, four, five;

	PROCEDURE inputstatement;
	INSPECT inputtext DO
	BEGIN	inputpart;
		WHILE lastatom == comma DO BEGIN nextatom; inputpart; END;
	END OTHERWISE error (20);

	PROCEDURE inputpart;
	IF lastatom == smaller THEN
	BEGIN	idtype:= 6; nextatom; expression;
		IF lastatom == greater THEN nextatom ELSE error (21); control;
	END ELSE INSPECT lastatom
	WHEN   identifier DO inputname
	WHEN       number DO inputvalue
	WHEN functionname DO BEGIN functiondesignator (lastatom); control; END
	WHEN   optionname DO BEGIN error (22); nextatom; control; END
	OTHERWISE description (FALSE);

	PROCEDURE inputname;
	INSPECT lastatom DO IF declared THEN inputvalue
	ELSE BEGIN load (sav, index); nextatom; type:= 6; END;

	PROCEDURE inputvalue;
	INSPECT lastatom DO BEGIN load (tav, index); nextatom; control; END;

	PROCEDURE control;
	IF lastatom == mulsign THEN
	BEGIN	INTEGER address; address:= programcounter;
		load (0, 0); load (stk, 0); nextatom; description (TRUE);
		program [address]:= programcounter * 100 + jmp;
	END ELSE load (eql, 0);

	PROCEDURE description (repeat); BOOLEAN repeat;
	BEGIN	INTEGER address; address:= programcounter;
		IF lastatom == lparent THEN
		BEGIN	nextatom; inputstatement;
			IF lastatom == rparent THEN nextatom ELSE error (23);
			IF repeat THEN load (rtn, address);
		END ELSE
		IF lastatom == lbracket THEN
		BEGIN	three:= four:= five:= FALSE;
			nextatom; variablelist;
			IF lastatom == rbracket THEN nextatom ELSE error (24);
			IF repeat THEN mixed (address)
			ELSE IF five THEN load (fix, address);
		END ELSE
		BEGIN error (25); nextatom; END;
	END *** description ***;

	PROCEDURE variablelist;
	BEGIN	variable;
		WHILE nextatom == comma DO BEGIN nextatom; variable; END;
	END *** variablelist ***;

	PROCEDURE variable;
	INSPECT lastatom WHEN identifier DO
	BEGIN	IF declared THEN
		BEGIN	IF type EQ 1 OR type EQ 3 THEN
			BEGIN	IF type EQ 1 THEN rightcounter:= rightcounter + 1;
				load (rav, index); type:= 3; three:= TRUE;
			END ELSE
			IF type EQ 2 OR type EQ 5 THEN
			BEGIN	IF type EQ 2 THEN leftcounter:= leftcounter + 1;
				load (lav, index); type:= 5; five:= TRUE;
			END ELSE
			IF type EQ 4 THEN load (sav, index) ELSE error (26);
		END ELSE
		BEGIN load (sav, index); type:= 4; four:= TRUE; END;
	END OTHERWISE error (27);

	PROCEDURE mixed (address); INTEGER address;
	IF NOT (three OR five) THEN skip (address) ELSE
	IF NOT (three OR four) THEN BEGIN load (rtn, address); load (fix, address); END
	ELSE BEGIN IF five THEN load (fix, address); load (rtn, address); END;

	PROCEDURE skip (address); INTEGER address;
	BEGIN	INTEGER count; count:= programcounter - address;
		programcounter:= address - 1; load (skp, count);
	END *** skip ***;

	COMMENT ************************************* option statement ***********************************;

	INTEGER optioninstructions, optionnumber, submodel;

	PROCEDURE optionstatement;
	BEGIN	INTEGER i; FOR i:= 1 STEP 1 UNTIL optionlimit DO optionwanted [i]:= FALSE;
		INSPECT optiontext DO IF lastatom =/= semicolon THEN
		BEGIN	option; WHILE lastatom == comma DO BEGIN nextatom; option; END;
		END ELSE currenttext:- optiontext:- NONE;
	END *** optionstatement ***;

	PROCEDURE option;
	BEGIN	INSPECT lastatom WHEN optionname DO
		BEGIN	optionnumber:= type; optionwanted [type]:= TRUE; END
		WHEN number DO
		BEGIN	optionnumber:= Scanint (entry); entry . Setpos (1);
			IF 1 LE optionnumber AND optionnumber LE optionlimit
			THEN optionwanted [optionnumber]:= TRUE ELSE error (30);
		END OTHERWISE error (31); IF nextatom == lparent THEN
		BEGIN	nextatom; specifier;
			WHILE nextatom == comma DO BEGIN nextatom; specifier; END;
			IF lastatom == rparent THEN nextatom ELSE error (32);
		END;
	END *** option ***;

	PROCEDURE specifier;
	INSPECT lastatom WHEN number DO IF optionnumber EQ 5 THEN
	BEGIN	submodel:= Abs (Scanint (entry)); entry . Setpos (1);
		IF submodel LE termcounter - 1 THEN load (submodel, 0) ELSE error (33);
	END ELSE IF optionnumber EQ 10 THEN
	BEGIN	stack [stacklength]:= Scanreal (entry); stacklength:= stacklength - 1;
		IF entry . Pos EQ 1 THEN error (33); entry . Setpos (1);
	END ELSE error (34) OTHERWISE error (35);

	COMMENT ************************************* translation ****************************************;

	PROCEDURE translate (sourcetext, statement, restofstatement, instructions); NAME instructions;
	REF (textstorage) sourcetext; PROCEDURE statement, restofstatement; INTEGER instructions;
	BEGIN	currenttext:- sourcetext; printtext; nextatom; statement; INSPECT currenttext DO
		BEGIN	IF lastatom =/= semicolon AND NOT (endline AND endtext) THEN error (36);
			WHILE lastatom =/= semicolon AND NOT (endline AND endtext) DO restofstatement;
			IF lastatom =/= semicolon THEN error (37);
		END; instructions:= programcounter - 1;
	END *** translate ***;

	PROCEDURE prepare;
	BEGIN	IF rightcounter EQ 0 THEN BEGIN rightcounter:= 1; error1 (40, 1); END;
		IF leftcounter EQ 0 THEN
		BEGIN leftcounter:= programcounter:= 1; load (tcv, 1); load (ini, 0); END;
		designcounter:= termcounter + 1; weightcounter:= termcounter + 2;
	END *** prepare ***;

	translate (modeltext, modelstatement, expression, modelinstructions);
	translate (inputtext, inputstatement, inputstatement, inputinstructions);
	translate (optiontext, optionstatement, optionstatement, optioninstructions);
	IF erroneous THEN GOTO endjob; prepare;

BEGIN	INTEGER stackpointer, termpointer, casepointer, repeats, repeatsize, designsize;
	LONG REAL weight, weightroot, weightsize, repeatroot, repeatsum, repeatsquares;

	REF (designstorage) ARRAY designdata [1 : weightcounter];
BEGIN	REF  (rightstorage) ARRAY  rightdata [1 :  rightcounter];
	REF   (leftstorage) ARRAY   leftdata [1 :   leftcounter];

	COMMENT ************************************* checking *******************************************;

	PROCEDURE checkmodel;
	BEGIN	INTEGER count, param, term, right, left,
		instruction, address, last, next; BOOLEAN constant;
		param:= term:= right:= left:= 0; constant:= FALSE;
		FOR count:= 3 STEP 1 UNTIL inputinstructions DO
		BEGIN	instruction:= program [count]; address:= instruction // 100;
			IF instruction GT 100 AND Mod (instruction, 100) EQ tav THEN
			BEGIN	INSPECT atomhash . hashtable [address] WHEN identifier DO
				IF type EQ 1 THEN
				BEGIN	IF term GE 1 THEN
					BEGIN	param:= param + 1; stack [address]:= 1;
						designdata [term] . address:= address;
						next:= program [count + 1];
						IF last EQ res AND next EQ res THEN constant:= TRUE;
						IF last EQ res EQV next EQ mul THEN error1 (41, term);
					END ELSE error1 (42, term);
				END ELSE
				IF type EQ 3 AND stack [address] EQ 0 THEN
				BEGIN	right:= right + 1; stack [address]:= right;
					rightdata [right] . address:= address;
				END ELSE
				IF type EQ 5 AND stack [address] EQ 0 THEN
				BEGIN	left:= left + 1; stack [address]:= left;
					leftdata [left] . address:= address;
				END ELSE
				IF type EQ 2 THEN error1 (43, term)
				WHEN number DO IF entry . Pos EQ 1 THEN
				BEGIN	stack [address]:= Scanreal (entry);
					IF entry . Pos EQ 1 THEN error1 (44, term);
				END;
			END ELSE
			IF instruction EQ res THEN
			BEGIN	IF last NE mul THEN
				BEGIN	IF constant THEN constant:= FALSE
					ELSE error1 (45, term)
				END; IF param GT 1 THEN error1 (46, term)
				ELSE IF param EQ 0 THEN error1 (47, term);
				param:= 0; term:= term + 1;
			END ELSE
			IF instruction EQ dep THEN BEGIN instruction:= res; term:= 1; END;
			last:= instruction;
		END;
		IF erroneous THEN GOTO endjob;
		IF NOT (dataread OR datafileread) THEN readdatafile;
		IF NOT dataread THEN INSPECT printstream DO
		BEGIN	errortext (IF dataskipped THEN 2 ELSE 1); INSPECT datastream DO
			BEGIN	printstream . Setpos (14); Outtext ("In record ");
				Outint (datarecord, field (datarecord));
				Outtext (" of file : "); Outtext (Filename (datastream) . Strip);
			END; Outimage; GOTO endjob;
		END;
		IF NOT optionwanted [7] THEN INSPECT inputdata DO
		BEGIN IF optionwanted [6] THEN printdata; reset; END;
		stackpointer:= tablelimit + 1; programcounter:= modelinstructions + 1;
	END *** checkmodel ***;

	PROCEDURE checkinput;
	BEGIN	INTEGER i; REAL size, check;
		IF NOT optionwanted [7] THEN INSPECT inputdata DO
		IF nextnumber NE size - 1 THEN INSPECT printstream DO
		BEGIN	Eject (Line + 5); Outtext ("Not all input data (in this");
			Outtext (" record) is processed    (message)"); Outimage;
		END;
		size:= leftdata [1] . size; IF size EQ 0 THEN error2 (70, size, 0);
		FOR i:= 2 STEP 1 UNTIL leftcounter DO
		BEGIN check:= leftdata [i] . size; IF size NE check THEN error2 (71, size, check); END;
		size:= rightdata [1] . size; IF size EQ 0 THEN error2 (72, size, 0);
		FOR i:= 2 STEP 1 UNTIL rightcounter DO
		BEGIN check:= rightdata [i] . size; IF size NE check THEN error2 (73, size, check); END;
		stack [stackpointer]:= size; stackpointer:= stackpointer + 1;
		repeatsize:= 0; weightsize:= 0; programcounter:= 1; casepointer:= size + 1;
	END *** checkinput ***;

	COMMENT ************************************* execution ******************************************;

	PROCEDURE execute (instructions); INTEGER instructions;
	BEGIN	INTEGER i, j, instruction, address; LONG REAL F, G; BOOLEAN missingvalue, missingcase;
		PROCEDURE checkmissing;
		FOR j:= stacklength + 1 STEP 1 UNTIL stackspace DO IF stack [j] EQ G THEN missingvalue:= TRUE;
		SWITCH macro:=	stk, add, sub, mul, div, idi, ttp, neg,
				eql, rep, dep, res, tst, ini, frs, fls, smc;
		SWITCH macro2:= tcv, tav, rav, sav, lav, rtn, jmp, fix, skp, fun;
	  next: instruction:= program [programcounter]; address:= instruction // 100;
		programcounter:= programcounter + 1;
		GOTO IF instruction LT 100 THEN macro  [instruction]
					   ELSE macro2 [Mod (instruction, 100)];
	   stk: IF stackpointer GT stacklength THEN overflow ("Stack space");
		stack [stackpointer]:= F; stackpointer:= stackpointer + 1; GOTO exit;
	   add: stackpointer:= stackpointer - 1; F:= stack [stackpointer] + F; GOTO exit;
	   sub: stackpointer:= stackpointer - 1; F:= stack [stackpointer] - F; GOTO exit;
	   mul: stackpointer:= stackpointer - 1; F:= stack [stackpointer] * F; GOTO exit;
	   div: IF F EQ 0 THEN error2 (50, termpointer, casepointer);
		stackpointer:= stackpointer - 1; F:= stack [stackpointer] / F; GOTO exit;
	   idi: IF F EQ 0 THEN error2 (51, termpointer, casepointer);
		stackpointer:= stackpointer - 1; F:= stack [stackpointer] // F; GOTO exit;
	   ttp: stackpointer:= stackpointer - 1; G:= stack [stackpointer]; I:= F;
		IF G EQ 0 AND F LT 0 THEN error2 (54, F, casepointer);
		IF G LT 0 AND F NE I THEN error2 (55, F, casepointer);
		F:= IF F EQ I THEN G ** I ELSE G ** F; GOTO exit;
	   neg: F:= - F; GOTO exit;
	   eql: G:= inputdata . nextnumber;
		IF Abs (F - G) GT precision * Abs (F) THEN error2 (74, F, G); GOTO exit;
	   rep: IF Abs (F) GT maxobs THEN error2 (52, 0, casepointer);
		INSPECT designdata [designcounter] DO
		BEGIN IF F LT min THEN min:= F; IF F GT max THEN max:= F; END; repeats:= repeats + 1;
		repeatsum:= repeatsum + F; repeatsquares:= repeatsquares + F * F; GOTO exit;
	   dep: IF F LE precision THEN error2 (56, F, casepointer);
		weight:= F; weightroot:= Sqrt (weight); repeatroot:= Sqrt (repeats);
		repeatsize:= repeatsize + repeats; weightsize:= weightsize + repeats * weight;
		designdata [weightcounter] . store (repeatroot * weightroot);
		INSPECT designdata [designcounter] DO
		BEGIN	store (repeatsum / repeatroot * weightroot);
			sum:= sum + repeatsum * weight; squares:= squares + repeatsquares * weight;
		END; termpointer:= 1; GOTO exit;
	   res: IF Abs (F) GT maxobs THEN error2 (53, termpointer, casepointer);
		INSPECT designdata [termpointer] DO
		BEGIN	store (F * repeatroot * weightroot);
			IF F LT min THEN min:= F; IF F GT max THEN max:= F;
			G:= F * repeats * weight; sum:= sum + G; squares:= squares + F * G;
		END; termpointer:= termpointer + 1; GOTO exit;
	   tst: F:= leftdata [1] . lastnumber; FOR i:= 2 STEP 1 UNTIL leftcounter DO
		BEGIN G:= leftdata [i] . lastnumber; IF F NE G THEN error2 (75, F, G); END;
	   ini: stack [stackpointer]:= F; stackpointer:= stackpointer + 1;
		repeats:= 0; repeatsum:= repeatsquares:= 0; casepointer:= casepointer - 1;
	   frs: missingvalue:= FALSE; FOR i:= 1 STEP 1 UNTIL rightcounter DO INSPECT rightdata [i] DO
		BEGIN G:= stack [address]:= lastnumber; IF optionwanted [10] THEN checkmissing; END;
		missingcase:= missingvalue; GOTO exit;
	   fls: missingvalue:= FALSE; FOR i:= 1 STEP 1 UNTIL leftcounter DO INSPECT leftdata [i] DO
		BEGIN G:= stack [address]:= lastnumber; IF optionwanted [10] THEN checkmissing; END;
		IF missingcase OR missingvalue THEN programcounter:= repeatinstruction; GOTO exit;
	   smc: IF missingcase OR repeats EQ 0 THEN programcounter:= modelinstructions; GOTO exit;
	   tcv: F:= address; GOTO exit;
	   tav: F:= stack [address]; GOTO exit;
	   rav: rightdata [stack [address]] . store (inputdata . nextnumber); GOTO exit;
	   sav: stack [address]:= inputdata . nextnumber; GOTO exit;
	   lav: leftdata [stack [address]] . store (inputdata . nextnumber); GOTO exit;
	   rtn: F:= stack [stackpointer - 1]; IF F LT 1.5 THEN stackpointer:= stackpointer - 1 ELSE
		BEGIN stack [stackpointer - 1]:= F - 1; programcounter:= address; END; GOTO exit;
	   jmp: IF Abs (F - Entier (F + precision * Abs (F))) GT 2 * precision * Abs (F) THEN error2
		(76, F, IF stackpointer EQ tablelimit + 1 THEN 0 ELSE stack [stackpointer - 1] - 1);
		IF F LT 0.5 THEN programcounter:= address; GOTO exit;
	   fix: FOR i:= programcounter - 2 STEP - 1 UNTIL address DO
		IF Mod (program [i], 100) EQ 5 THEN leftdata [stack [program [i] // 100]] . fixrepeats; GOTO exit;
	   skp: inputdata . skipnumbers (F * address); GOTO exit;
	   fun: BEGIN	SWITCH function:= fabs, fsign, fsqrt, fsin, fcos, ftan, fln, flog, fexp,
			fent, frnd, fmod, fmin, fmax, fasin, facos, fatan, fsinh, fcosh, ftanh, find;
			GOTO function [address];
		  fabs: F:= Abs (F); GOTO exit;
		 fsign: F:= Sign (F); GOTO exit;
		 fsqrt: IF F LT 0 THEN error2 (60, F, casepointer);
			F:= Sqrt (F); GOTO exit;
		  fsin: F:= Sin (F); GOTO exit;
		  fcos: F:= Cos (F); GOTO exit;
		  ftan: F:= Tan (F); GOTO exit;
		   fln: IF F LE 0 THEN error2 (61, F, casepointer);
			F:= Ln (F); GOTO exit;
		  flog: IF F LE 0 THEN error2 (62, F, casepointer);
			F:= Ln (F) / Ln10; GOTO exit;
		  fexp: IF F GT maxexp THEN error2 (63, F, casepointer);
			F:= IF F LT - maxexp THEN 0 ELSE Exp (F); GOTO exit;
		  fent: F:= Entier (F); GOTO exit;
		  frnd: F:= Entier (F + 0.5); GOTO exit;
		  fmod: stackpointer:= stackpointer - 1;
			F:= Mod (stack [stackpointer], F); GOTO exit;
		  fmin: stackpointer:= stackpointer - 1;
			F:= Lmin (stack [stackpointer], F); GOTO exit;
		  fmax: stackpointer:= stackpointer - 1;
			F:= Lmax (stack [stackpointer], F); GOTO exit;
		 fasin: IF Abs (F) GT 1 THEN error2 (64, F, casepointer);
			F:= Arcsin (F); GOTO exit;
		 facos: IF Abs (F) GT 1 THEN error2 (65, F, casepointer);
			F:= Arccos (F); GOTO exit;
		 fatan: F:= Arctan (F); GOTO exit;
		 fsinh: IF Abs (F) GT maxexp THEN error2 (66, F, casepointer);
			F:= Sinh (F); GOTO exit;
		 fcosh: IF Abs (F) GT maxexp THEN error2 (67, F, casepointer);
			F:= Cosh (F); GOTO exit;
		 ftanh: F:= Tanh (F); GOTO exit;
		  find: stackpointer:= stackpointer - 2; G:= stack [stackpointer + 1];
			F:= IF stack [stackpointer] - G LE precision * Abs (G) AND
			    G - F LE precision * Abs (G) THEN 1 ELSE 0; GOTO exit;
		END *** fun ***;
	  exit: IF programcounter LE instructions THEN GOTO next;
	END *** execute ***;

	PROCEDURE initializestorage;
	BEGIN	INTEGER i;
		FOR i:= 1 STEP 1 UNTIL weightcounter DO designdata [i]:- NEW designstorage;
		FOR i:= 1 STEP 1 UNTIL  rightcounter DO  rightdata [i]:- NEW  rightstorage;
		FOR i:= 1 STEP 1 UNTIL   leftcounter DO   leftdata [i]:- NEW   leftstorage;
	END *** initializestorage ***;

	initializestorage;
	checkmodel; execute (inputinstructions);
	checkinput; execute (modelinstructions);
	designsize:= designdata [1] . size;

END *** execution ***;

BEGIN	INTEGER ARRAY  pivot [1 : termcounter];
	LONG REAL ARRAY diag [1 : termcounter], resprod [1 : designcounter],
			data [1 : designsize, 1 : weightcounter], depvar [1 : designsize];

	COMMENT ************************************* print layout ***************************************;

	INTEGER low, upp, max, pagenr, partnr, matrixwidth;

	PROCEDURE printname (col); INTEGER col;
	INSPECT printstream DO
	IF col EQ weightcounter THEN
	BEGIN	IF repeatsize EQ weightsize THEN
		Outtext (" repeats") ELSE Outtext ("  weight");
	END ELSE IF col EQ designcounter THEN Outtext ("dep.var.")
	ELSE printitem (atomhash . hashtable [designdata [col] . address]);

	PROCEDURE underline (width); INTEGER width;
	INSPECT printstream DO BEGIN WHILE Pos LE width DO Outchar ('='); Outimage; END;

	PROCEDURE dashedline;
	INSPECT printstream DO
	BEGIN	Eject (Line + 1); Outtext ("----------------------------------");
		Outtext ("----------------------------------------------------");
		Outtext ("-------------------------"); Outimage; Eject (Line + 1);
	END *** dashedline ***;

	PROCEDURE heading;
	INSPECT printstream DO
	BEGIN	INTEGER j; Setpos (14);
		FOR j:= low STEP 1 UNTIL upp DO
		BEGIN printname (j); Setpos (Pos // 12 * 12 + 14); END;
		Outimage; Eject (Line + 1);
	END *** heading ***;

	PROCEDURE newpage1;
	INSPECT printstream DO
	BEGIN	Eject (1); Outtext ("Transformed data matrix");
		IF designsize GT 50 OR designcounter GT 10 THEN
		BEGIN	pagenr:= pagenr + 1;
			Outtext ("  -  page"); Outint (pagenr, 3);
		END; Outimage; underline (23); Eject (Line + 1);
		Outtext ("obs.no."); heading;
	END *** newpage1 ***;

	PROCEDURE newpart (string); NAME string; TEXT string;
	INSPECT printstream DO
	BEGIN	IF Line GT 49 - termcounter + low
		THEN Eject (1) ELSE Eject (Line + 5);
		Outtext ("Correlation matrix of the "); Outtext (string);
		IF termcounter GT 10 AND designcounter GT 10 THEN
		BEGIN	partnr:= partnr + 1;
			Outtext ("  -  part"); Outint (partnr, 3);
		END; Outimage; underline (35); Eject (Line + 1); heading;
	END *** newpart ***;

	PROCEDURE newpage2;
	INSPECT printstream DO
	BEGIN	Eject (1); Outtext ("Residual analysis");
		IF designsize GT 50 THEN
		BEGIN	pagenr:= pagenr + 1;
			Outtext ("  -  page"); Outint (pagenr, 3);
		END; Outimage; underline (17); Setpos (95);
		Outtext ("standardized         studentized"); Outimage;
		Outtext ("obs.no.        observation        fitted value");
		Outtext ("    standard deviation          residual");
		Outtext ("            residual            residual");
		Outimage; Eject (Line + 1);
	END *** newpage2 ***;

	COMMENT ************************************* Dekker matrix routines *****************************;

	LONG REAL PROCEDURE vecvec (l, u, a, b);
	INTEGER l, u; LONG REAL ARRAY a, b;
	BEGIN	INTEGER k; LONG REAL s; s:= 0;
		FOR k:= l STEP 1 UNTIL u DO s:= a [k] * b [k] + s;
		vecvec:= s;
	END *** vecvec ***;

	LONG REAL PROCEDURE matvec (l, u, i, a, b);
	INTEGER l, u, i; LONG REAL ARRAY a, b;
	BEGIN	INTEGER k; LONG REAL s; s:= 0;
		FOR k:= l STEP 1 UNTIL u DO s:= a [i,k] * b [k] + s;
		matvec:= s;
	END *** matvec ***;

	LONG REAL PROCEDURE tamvec (l, u, i, a, b);
	INTEGER l, u, i; LONG REAL ARRAY a, b;
	BEGIN	INTEGER k; LONG REAL s; s:= 0;
		FOR k:= l STEP 1 UNTIL u DO s:= a [k,i] * b [k] + s;
		tamvec:= s;
	END *** tamvec ***;

	LONG REAL PROCEDURE tammat (l, u, i, j, a, b);
	INTEGER l, u, i, j; LONG REAL ARRAY a, b;
	BEGIN	INTEGER k; LONG REAL s; s:= 0;
		FOR k:= l STEP 1 UNTIL u DO s:= a [k,i] * b [k,j] + s;
		tammat:= s;
	END *** tammat ***;

	PROCEDURE elmcol (l, u, i, j, a, b, x);
	INTEGER l, u, i, j; LONG REAL x; LONG REAL ARRAY a, b;
	FOR l:= l STEP 1 UNTIL u DO a [l,i]:= a [l,i] + b [l,j] * x;

	PROCEDURE elmveccol (l, u, i, a, b, x);
	INTEGER l, u, i; LONG REAL x; LONG REAL ARRAY a, b;
	FOR l:= l STEP 1 UNTIL u DO a [l]:= a [l] + b [l,i] * x;

	PROCEDURE ichcol (l, u, i, j, a);
	INTEGER l, u, i, j; LONG REAL ARRAY a;
	BEGIN	LONG REAL r;
		FOR l:= l STEP 1 UNTIL u DO
		BEGIN r:= a [l,i]; a [l,i]:= a [l,j]; a [l,j]:= r; END;
	END *** ichcol ***;

	PROCEDURE ichrow (l, u, i, j, a);
	INTEGER l, u, i, j; LONG REAL ARRAY a;
	BEGIN	LONG REAL r;
		FOR l:= l STEP 1 UNTIL u DO
		BEGIN r:= a [i,l]; a [i,l]:= a [j,l]; a [j,l]:= r; END;
	END *** ichrow ***;

	PROCEDURE ichrowcol (l, u, i, j, a);
	INTEGER l, u, i, j; LONG REAL ARRAY a;
	BEGIN	LONG REAL r;
		FOR l:= l STEP 1 UNTIL u DO
		BEGIN r:= a [i,l]; a [i,l]:= a [l,j]; a [l,j]:= r; END;
	END *** ichrowcol ***;

	PROCEDURE chlinv2 (a, n); INTEGER n; LONG REAL ARRAY a;
	BEGIN	LONG REAL r; INTEGER i, j, i1;
		LONG REAL ARRAY u [1 : n];
		FOR i:= n STEP - 1 UNTIL 1 DO
		BEGIN	r:= 1 / a [i,i]; i1:= i + 1;
			FOR j:= i1 STEP 1 UNTIL n DO u [j]:= a [i,j];
			FOR j:= n STEP - 1 UNTIL i1 DO
			a [i,j]:= - (tamvec (i1, j, j, a, u) + matvec (j + 1, n, j, a, u)) * r;
			a [i,i]:= (r - matvec (i1, n, i, a, u)) * r;
		END;
	END *** chlinv2 ***;

	INTEGER PROCEDURE lsqdec (a, n, m, tol, aid, ci);
	INTEGER n, m; LONG REAL tol; LONG REAL ARRAY a, aid; INTEGER ARRAY ci;
	BEGIN	INTEGER j, k, kpiv;
		LONG REAL beta, sigma, norm, w, eps, akk, aidk;
		LONG REAL ARRAY sum [1 : m];
		norm:= 0; lsqdec:= m;
		FOR k:= 1 STEP 1 UNTIL m DO
		BEGIN	w:= sum [k]:= tammat (1, n, k, k, a, a);
			IF w GT norm THEN norm:= w;
		END;
		w:= Sqrt (norm); eps:= tol * w;
		FOR k:= 1 STEP 1 UNTIL m DO
		BEGIN	sigma:= sum [k]; kpiv:= k;
			FOR j:= k + 1 STEP 1 UNTIL m DO
			IF sum [j] GT sigma THEN
			BEGIN sigma:= sum [j]; kpiv:= j; END;
			IF kpiv NE k THEN
			BEGIN sum [kpiv]:= sum [k]; ichcol (1, n, k, kpiv, a); END;
			ci [k]:= kpiv; akk:= a [k,k];
			sigma:= tammat (k, n, k, k, a, a); w:= Sqrt (sigma);
			aidk:= aid [k]:= IF akk LT 0 THEN w ELSE - w;
			IF w LT eps THEN BEGIN lsqdec:= k - 1; GOTO exit; END;
			beta:= 1 / (sigma - akk * aidk); a [k,k]:= akk - aidk;
			FOR j:= k + 1 STEP 1 UNTIL m DO
			BEGIN	elmcol (k, n, j, k, a, a, - beta * tammat (k, n, k, j, a, a));
				sum [j]:= sum [j] - a [k,j] ** 2;
			END;
		END;	exit:
	END *** lsqdec ***;

	PROCEDURE lsqsol (a, n, m, aid, ci, b);
	INTEGER n, m; LONG REAL ARRAY a, aid, b; INTEGER ARRAY ci;
	BEGIN	INTEGER k, cik; LONG REAL w;
		FOR k:= 1 STEP 1 UNTIL m DO
		elmveccol (k, n, k, b, a, tamvec (k, n, k, a, b) / (aid [k] * a [k,k]));
		FOR k:= m STEP - 1 UNTIL 1 DO
		b [k]:= (b [k] - matvec (k + 1, m, k, a, b)) / aid [k];
		FOR k:= m STEP - 1 UNTIL 1 DO
		BEGIN	cik:= ci [k]; IF cik NE k THEN
			BEGIN w:= b [k]; b [k]:= b [cik]; b [cik]:= w; END;
		END;
	END *** lsqsol ***;

	PROCEDURE lsqinv (a, m, aid, c);
	INTEGER m; LONG REAL ARRAY a, aid; INTEGER ARRAY c;
	BEGIN	INTEGER i, ci; LONG REAL w;
		FOR i:= 1 STEP 1 UNTIL m DO a [i,i]:= aid [i];
		chlinv2 (a, m);
		FOR i:= m STEP - 1 UNTIL 1 DO
		BEGIN	ci:= c [i]; IF ci NE i THEN
			BEGIN	ichcol (1, i - 1, i, ci, a);
				ichrowcol (i + 1, ci - 1, i, ci, a);
				ichrow (ci + 1, m, i, ci, a);
				w:= a [i,i]; a [i,i]:= a [ci,ci]; a [ci,ci]:= w;
			END;
		END;
	END *** lsqinv ***;

	COMMENT ************************************* distribution functions *****************************;

	REAL PROCEDURE Phi (x); REAL x;
	BEGIN	REAL y, z, w; y:= Abs (x) / 2;
		IF y GE 3 THEN z:= 1 ELSE IF y LT 1 THEN
		BEGIN	w:= y * y;
			z:= (((((((( 0.000124818987 * w
			    - 0.001075204047) * w + 0.005198775019) * w
			    - 0.019198292004) * w + 0.059054035642) * w
			    - 0.151968751364) * w + 0.319152932694) * w
			    - 0.531923007300) * w + 0.797884560593) * y * 2
		END ELSE
		BEGIN	y:= y - 2;
			z:= ((((((((((((( - 0.000045255659 * y
			    + 0.000152529290) * y - 0.000019538132) * y
			    - 0.000676904986) * y + 0.001390604284) * y
			    - 0.000794620820) * y - 0.002034254874) * y
			    + 0.006549791214) * y - 0.010557625006) * y
			    + 0.011630447319) * y - 0.009279453341) * y
			    + 0.005353579108) * y - 0.002141268741) * y
			    + 0.000535310849) * y + 0.999936657524
		END;
		Phi:= IF x GT 0 THEN (z + 1) / 2 ELSE (1 - z) / 2
	END *** Phi ***;

	REAL PROCEDURE Fisher (x, df1, df2); REAL x; INTEGER df1, df2;
	IF x LE 0 THEN Fisher:= 0 ELSE
	IF x GE 10000 THEN Fisher:= 1 ELSE
	IF df1 GT 1000 AND df2 GT 1000 THEN
	Fisher:= Phi ((x ** (1/3) * (1 - 2/df2/9) + 2/df1/9 - 1) /
		 Sqrt (x ** (2/3) *      2/df2/9  + 2/df1/9)) ELSE
	BEGIN	INTEGER a, b, i, j; REAL w, y, z, zk, d, p;
		a:= df1 // 2 * 2 - df1 + 2; b:= df2 // 2 * 2 - df2 + 2;
		w:= x * df1 / df2; z:= 1 / (1 + w);
		IF a EQ 1 THEN
		BEGIN	IF b EQ 1 THEN
			BEGIN	p:= Sqrt (w); y:= 0.318309886184;
				d:= y * z / p; p:= Arctan (p) * 2 * y;
			END ELSE
			BEGIN p:= Sqrt (w * z); d:= 0.5 * p * z / w; END;
		END ELSE
		IF b EQ 1 THEN
		BEGIN p:= Sqrt (z); d:= 0.5 * z * p; p:= 1 - p; END
		ELSE BEGIN d:= z * z; p:= w * z; END;
		y:= 2 * w / z;
		IF a EQ 1 THEN
		BEGIN	FOR j:= b + 2 STEP 2 UNTIL df2 DO
			BEGIN	d:= (a / (j - 2) + 1) * d * z;
				p:= d * y / (j - 1) + p;
			END;
		END ELSE
		BEGIN	zk:= z ** ((df2 - 1) // 2); d:= d * zk * df2 / b;
			p:= p * zk + w * z * (zk - 1) / (z - 1);
		END;
		y:= w * z; z:= 2 / z; b:= df2 - 2;
		FOR i:= a + 2 STEP 2 UNTIL df1 DO
		BEGIN	j:= i + b; d:= y * d * j / (i - 2);
			p:= p - z * d / j;
		END;
		Fisher:= IF p LT 0 THEN 0 ELSE IF p GT 1 THEN 1 ELSE p;
	END *** Fisher ***;

	COMMENT ************************************* transformed data matrix ****************************;

	PROCEDURE transformed_data_matrix;
	IF optionwanted [1] AND submodel EQ 0 THEN
	BEGIN	INTEGER i, j, datacounter; datacounter:=
		IF designsize EQ weightsize THEN designcounter ELSE weightcounter;
		INSPECT outputstream DO
		BEGIN	Outint (designsize, 12); Outint (datacounter, 12); Outimage;
			FOR i:= 1 STEP 1 UNTIL designsize DO
			BEGIN	weightroot:= data [i, weightcounter];
				FOR j:= 1 STEP 1 UNTIL datacounter DO IF j NE weightcounter
				THEN Outreal (data [i,j] / weightroot, 16, 22)
				ELSE Outreal (weightroot * weightroot, 16, 22); Outimage;
			END;
		END;
		pagenr:= 0; upp:= 0;
		INSPECT printstream DO
		FOR low:= upp + 1 WHILE upp LT datacounter DO
		BEGIN	upp:= IF datacounter - low GT matrixwidth
			THEN low + matrixwidth ELSE datacounter; newpage1;
			FOR i:= 1 STEP 1 UNTIL designsize DO
			BEGIN	Outint (i, 6); Setpos (10);
				weightroot:= data [i, weightcounter];
				FOR j:= low STEP 1 UNTIL upp DO IF j NE weightcounter
				THEN Outfix (data [i,j] / weightroot, 3, 12)
				ELSE Outfix (weightroot * weightroot, 3, 12); Outimage;
				IF Mod (i, 10) EQ 0 AND i LT designsize THEN
				BEGIN IF Mod (i, 50) EQ 0 THEN newpage1 ELSE Eject (Line + 1); END;
			END;
		END;
	END *** transformed_data_matrix ***;

	COMMENT ************************************* control information ********************************;

	INTEGER minterm, constterm; LONG REAL standdev;
	BOOLEAN constant, depvarconstant, oneconstant, moreconstants;

	PROCEDURE control_information;
	INSPECT printstream DO
	BEGIN	INTEGER j; constterm:= 0;
		constant:= depvarconstant:= oneconstant:= moreconstants:= FALSE;
		Eject (1); Outtext ("Control information");
		IF submodel GT 0 THEN
		BEGIN Outtext ("  -  submodel"); Outint (submodel, 3); END;
		Outimage; underline (19); Eject (Line + 1);
		Outtext ("transformed variable"); Outimage;
		Outtext ("denoted by parameter          mean           ");
		Outtext ("standard deviation                minimum    ");
		Outtext ("              maximum"); Outimage; Eject (Line + 1);
		FOR j:= termcounter + 1 STEP 1 UNTIL designcounter - 1 DO
		BEGIN	Setpos (2); printname (j); Setpos (14);
			Outtext ("omitted"); Outimage;
		END;
	 	IF submodel GT 0 THEN Eject (Line + 1);
		FOR j:= 1 STEP 1 UNTIL termcounter, designcounter DO
		INSPECT designdata [j] DO
		BEGIN	Setpos (2); printname(j); Setpos (12);
			standdev:= squares - sum * sum / weightsize;
			standdev:= IF standdev LT precision * squares OR weightsize LE 1
				   THEN 0 ELSE Sqrt (standdev / (weightsize - 1));
			Outfix (sum / weightsize, 6, 25); Outfix (standdev, 6, 25);
			Outfix (min, 6, 25); Outfix (max, 6, 25); Outimage;
			IF standdev EQ 0 THEN
			BEGIN	depvarconstant:= j EQ designcounter;
				oneconstant:= termcounter EQ 1; constterm:= j;
				IF NOT depvarconstant THEN
				BEGIN moreconstants:= constant; constant:= TRUE; END;
			END;
		END;
		Eject (Line + 1); Outtext ("Number of observations :  ");
		Outint (designsize, Imax (field (designsize), 3)); Outimage;
		minterm:= IF constterm EQ 1 THEN 2 ELSE 1;
		IF NOT (constant OR optionwanted [4]) THEN
		BEGIN	Eject (Line + 5); Outtext ("There is no constant independent");
			Outtext (" variable in the transformed (sub)model    (message)"); Outimage;
		END;
	END *** control_information ***;

	COMMENT ************************************* correlation matrix *********************************;

	PROCEDURE correlation_matrix;
	IF optionwanted [2] AND submodel EQ 0 THEN
	BEGIN	INTEGER i, j, pj; LONG REAL resdevj, correlation;
		LONG REAL ARRAY resdev [1 : designcounter],
		correl [1 : (designcounter + 1) * designcounter // 2];
		FOR j:= 1 STEP 1 UNTIL designcounter DO
		INSPECT designdata [j] DO
		BEGIN	resdevj:= squares - sum * sum / weightsize;
			resdevj:= resdev [j]:= IF resdevj LT precision * squares
				  THEN 0 ELSE Sqrt (resdevj * weightsize);
			pj:= (j - 1) * j // 2; correl [pj + j]:= 1;
			FOR i:= j - 1 STEP - 1 UNTIL 1 DO correl [pj + i]:=
			IF resdev [i] EQ 0 OR resdevj EQ 0 THEN Maxreal ELSE
			(tammat (1, designsize, i, j, data, data) * weightsize
			 - designdata [i] . sum * sum) / (resdev [i] * resdevj);
		END;
		INSPECT outputstream DO
		BEGIN	Outint (designcounter, 12); Outimage;
			FOR j:= 1 STEP 1 UNTIL designcounter DO
			BEGIN	pj:= (j - 1) * j // 2;
				FOR i:= 1 STEP 1 UNTIL j DO
				Outreal (correl [pj + i], 16, 22);
			END; Outimage;
		END;
		partnr:= 0; upp:= 0;
		INSPECT printstream DO
		FOR low:= upp + 1 WHILE upp LT designcounter DO
		BEGIN	upp:= IF designcounter - low GT matrixwidth
			THEN low + matrixwidth ELSE designcounter;
			newpart ("variables");
			FOR j:= low STEP 1 UNTIL designcounter DO
			BEGIN	Setpos (2); printname (j); Setpos (10);
				pj:= (j - 1) * j // 2; max:= Imin (j, upp);
				FOR i:= low STEP 1 UNTIL max DO
				BEGIN	correlation:= correl [pj + i];
					IF Abs (correlation) LT 1.0000005
					THEN Outfix (correlation, 6, 12)
					ELSE Outtext ("        *   ");
				END; Outimage;
			END;
		END;
	END *** correlation_matrix ***;

	COMMENT ************************************* no regression analysis *****************************;

	PROCEDURE no_regression_analysis;
	IF optionwanted [4] OR designsize LE termcounter
	OR depvarconstant OR oneconstant OR moreconstants THEN
	INSPECT printstream DO
	BEGIN	Eject (Line + 5); Outtext ("No regression analysis");
		IF NOT optionwanted [4] THEN Outtext (" possible") ELSE
		BEGIN Outtext (" by option"); Outimage; GOTO endjob; END; Outimage;
		IF depvarconstant THEN
		BEGIN	Eject (Line + 2); Outtext ("The dependent");
			Outtext (" variable in the transformed model is constant"); Outimage;
		END;
		IF oneconstant THEN
		BEGIN	Eject (Line + 2); Outtext ("The only independent");
			Outtext (" variable in the transformed model is constant"); Outimage;
		END;
		IF moreconstants THEN
		BEGIN	Eject (Line + 2); Outtext ("There are several constant");
			Outtext (" independent variables in the transformed model"); Outimage;
		END;
		IF designsize LE termcounter THEN
		BEGIN	Eject (Line + 2); Outtext ("The number of respondents");
			Outtext (" is less than or equal to the number of");
			Outtext (" independent variables in the transformed model"); Outimage;
		END;
		erroneous:= TRUE; GOTO endjob;
	END *** no_regression_analysis ***;

	COMMENT ************************************* least squares **************************************;

	INTEGER dftotal, dfmean, dfreg, dfres, dflack, dferror, dfred;
      LONG REAL sstotal, ssmean, ssreg, ssres, sslack, sserror, ssred, explained, adjusted,
		msreg, msres, mslack, mserror, msred, frpar, frmean, frreg, frlack, frred;
	BOOLEAN	perfectfit, lackoffit, reduction;

	PROCEDURE least_squares;
	BEGIN	INTEGER rank;
		rank:= lsqdec (data, designsize, termcounter, precision, diag, pivot);
		IF rank LT termcounter THEN INSPECT printstream DO
		BEGIN	Eject (Line + 5); Outtext ("No regression analysis possible"); Outimage;
			Eject (Line + 2); Outtext ("The rank of the transformed data matrix is ");
			Outint (termcounter - rank, field (termcounter - rank));
			Outtext (" less than the number of independent variables in the transformed model");
			Outimage; erroneous:= TRUE; GOTO endjob;
		END;
		lsqsol (data, designsize, termcounter, diag, pivot, depvar);
		lsqinv (data, termcounter, diag, pivot);
		perfectfit:= lackoffit:= FALSE;
		dftotal:= repeatsize; sstotal:= designdata [designcounter] . squares;
		dfreg:= termcounter; ssreg:= vecvec (1, termcounter, depvar, resprod);
		dfres:= dftotal - dfreg; ssres:= sstotal - ssreg;
		IF ssres LE precision * sstotal THEN BEGIN ssres:= 0; perfectfit:= TRUE; END;
		IF constant THEN
		BEGIN	ssmean:= designdata [designcounter] . sum ** 2 / weightsize;
			dfmean:= 1; dfreg:= dfreg - dfmean; ssreg:= ssreg - ssmean;
		END ELSE BEGIN dfmean:= 0; ssmean:= 0; END;
		explained:= 1 - ssres / (sstotal - ssmean);
		 adjusted:= 1 - ssres / (sstotal - ssmean) * (dftotal - dfmean) / dfres;
		msreg:= ssreg / dfreg; msres:= ssres / dfres;
		frmean:= IF perfectfit THEN Maxreal ELSE ssmean / msres;
		frreg:= IF perfectfit THEN Maxreal ELSE msreg / msres;
		dferror:= repeatsize - designsize;
		IF dferror GT 0 THEN
		BEGIN	sserror:= sstotal - resprod [designcounter];
			IF sserror LE precision * sstotal THEN BEGIN sserror:= 0; lackoffit:= TRUE; END;
			sslack:= ssres - sserror; IF sslack LE precision * ssres THEN sslack:= 0;
			dflack:= dfres - dferror; mslack:= sslack / dflack; mserror:= sserror / dferror;
			frlack:= IF perfectfit OR lackoffit THEN 0 ELSE mslack / mserror;
		END;
		IF submodel EQ 0 AND (optionwanted [5] OR optionwanted [8]) THEN
		BEGIN	savemodel:= TRUE; savefit:= perfectfit;
			savedf:= dfres; savess:= ssres; savems:= msres;
		END;
		reduction:= IF (submodel GT 0 OR optionwanted [9]) AND savemodel
			    THEN dfres GT savedf AND ssres GT savess ELSE FALSE;
		IF reduction THEN
		BEGIN	dfred:= dfres - savedf; ssred:= ssres - savess; msred:= ssred / dfred;
			frred:= IF savefit THEN Maxreal ELSE msred / savems;
		END;
	END *** least_squares ***;

	COMMENT ************************************* multiple correlation *******************************;

	PROCEDURE multiple_correlation;
	IF adjusted GT precision THEN
	INSPECT printstream DO
	BEGIN	Eject (Line + 5); Outtext ("Multiple correlation coefficient");
		Outfix (Sqrt (explained), 6, 13); Outtext ("    (adjusted");
		Outfix (Sqrt ( adjusted), 6, 11); Outchar (')'); Outimage; underline (32);
	END *** multiple_correlation ***;

	COMMENT ************************************* explained variation ********************************;

	PROCEDURE explained_variation;
	INSPECT printstream DO
	BEGIN	Eject (Line + 5); Outtext ("Proportion of variation explained");
		Outfix (explained, 6, 12); Outtext ("    (adjusted");
		Outfix ( adjusted, 6, 11); Outchar (')'); Outimage; underline (33);
	END *** explained_variation ***;

	COMMENT ************************************* standard error deviation ***************************;

	PROCEDURE standard_error_deviation;
	INSPECT printstream DO
	BEGIN	Eject (Line + 5); Outtext ("Standard deviation of the error term");
		Outfix (Sqrt (msres), 6, 22); Outimage; underline (36);
	END *** standard_error_deviation ***;

	COMMENT ************************************* regression parameters ******************************;

	PROCEDURE regression_parameters;
	INSPECT printstream DO
	BEGIN	INTEGER j;
		INSPECT outputstream DO
		BEGIN	IF submodel EQ 0 THEN
			BEGIN Outint (subselection + 1, 12); Outimage; END;
			IF submodel EQ 0 OR subselection GT 0 THEN
			BEGIN	Outint (termcounter, 12); Outimage;
				FOR j:= 1 STEP 1 UNTIL termcounter DO
				Outreal (depvar [j], 16, 22); Outimage;
			END;
		END;
		Eject (1); Outtext ("Regression parameters"); Outimage;
		underline (21); Setpos (101); Outtext ("right  tail"); Outimage;
		Outtext ("parameter                 estimate           ");
		Outtext ("standard deviation              F - ratio    ");
		Outtext ("          probability"); Outimage; Eject (Line + 1);
		FOR j:= 1 STEP 1 UNTIL termcounter DO
		BEGIN	Setpos (2); printname (j); Setpos (12); standdev:= Sqrt (data [j,j] * msres);
			frpar:= IF perfectfit THEN Maxreal ELSE depvar [j] ** 2 / standdev / standdev;
			Outfix (depvar [j], 10, 25); Outfix (standdev, 10, 25); Outfix (frpar, 6, 25);
			Outfix (1 - Fisher (frpar, 1, dfres), 6, 25); Outimage; Eject (Line + 1);
		END;
	END *** regression_parameters ***;

	COMMENT ************************************* covariance matrix **********************************;

	PROCEDURE covariance_matrix;
	IF optionwanted [2] AND (submodel EQ 0 OR subselection GT 0) THEN
	BEGIN	INTEGER i, j; LONG REAL datajj;
		INSPECT outputstream DO
		BEGIN	Outint (termcounter, 12); Outimage;
			FOR j:= 1 STEP 1 UNTIL termcounter DO
			FOR i:= 1 STEP 1 UNTIL j DO
			Outreal (data [i,j], 16, 22); Outimage;
		END;
		partnr:= 0; upp:= 0;
		INSPECT printstream DO
		FOR low:= upp + 1 WHILE upp LT termcounter DO
		BEGIN	upp:= IF termcounter - low GT matrixwidth
			THEN low + matrixwidth ELSE termcounter;
			newpart ("estimates");
			FOR j:= low STEP 1 UNTIL termcounter DO
			BEGIN	Setpos (2); printname (j); Setpos (10);
				datajj:= data [j,j]; max:= Imin (j, upp);
				FOR i:= low STEP 1 UNTIL max DO
				Outfix (data [i,j] / Sqrt (data [i,i] * datajj), 6, 12);
				Outimage;
			END;
		END;
	END *** covariance_matrix ***;

	COMMENT ************************************* analysis of variance *******************************;

	PROCEDURE analysis_of_variance;
	INSPECT printstream DO
	BEGIN	Eject (Line + 5); Outtext ("Analysis of variance"); Outimage;
		underline (20); Eject (Line + 1); Outtext ("source of");
		Setpos (101); Outtext ("right  tail"); Outimage;
		Outtext ("variation          df          sum of squares");
		Outtext ("           mean square             F - ratio ");
		Outtext ("          probability"); Outimage; dashedline;
		Outtext ("total"); Setpos (16); Outint (dftotal, 6);
		Outfix (sstotal, 6, 24); Outimage; dashedline;
		IF constant THEN
		BEGIN	Outtext ("mean"); Setpos (16); Outint (dfmean, 6);
			Outfix (ssmean, 6, 24); Outfix (ssmean, 6, 22); Outfix (frmean, 6, 22);
			Outfix (1 - Fisher (frmean, dfmean, dfres), 6, 22); Outimage;
		END;
		Outtext ("regression"); Setpos (16); Outint (dfreg, 6);
		Outfix (ssreg, 6, 24); Outfix (msreg, 6, 22); Outfix (frreg, 6, 22);
		Outfix (1 - Fisher (frreg, dfreg, dfres), 6, 22); Outimage;
		Outtext ("residual"); Setpos (16); Outint (dfres, 6);
		Outfix (ssres, 6, 24); Outfix (msres, 6, 22); Outimage; dashedline;
		IF dferror GT 0 THEN
		BEGIN	Outtext (" lack of fit"); Setpos (16); Outint (dflack, 6);
			Outfix (sslack, 6, 24); Outfix (mslack, 6, 22); Outfix (frlack, 6, 22);
			Outfix (1 - Fisher (frlack, dflack, dferror), 6, 22); Outimage;
			Outtext (" pure error"); Setpos (16); Outint (dferror, 6);
			Outfix (sserror, 6, 24); Outfix (mserror, 6, 22); Outimage; dashedline;
		END;
		IF reduction THEN
		BEGIN	Outtext (" reduction"); Setpos (16); Outint (dfred, 6);
			Outfix (ssred, 6, 24); Outfix (msred, 6, 22); Outfix (frred, 6, 22);
			Outfix (1 - Fisher (frred, dfred, savedf), 6, 22); Outimage; dashedline;
		END;
	END *** analysis_of_variance ***;

	COMMENT ************************************* null hypotheses ************************************;

	PROCEDURE null_hypotheses;
	INSPECT printstream DO
	BEGIN	INTEGER j;
		Eject (Line + 1); Outtext ("regression null hypothesis :  ");
		FOR j:= 1 STEP 1 UNTIL termcounter DO
		BEGIN	IF j NE constterm THEN
			BEGIN printname (j); Outtext (" = "); END;
			IF Mod (j, 8) EQ 0 AND j LT termcounter THEN
			BEGIN Outimage; Setpos (31); END;
		END;
		Outchar ('0'); IF submodel GT 0 THEN Outtext ("  (in the reduced model)"); Outimage;
		IF submodel GT 0 AND reduction THEN
		BEGIN	Eject (Line + 1); Outtext (" reduction null hypothesis :  ");
			FOR j:= termcounter + 1 STEP 1 UNTIL designcounter - 1 DO
			BEGIN	printname (j); Outtext (" = ");
				IF Mod (j - termcounter, 8) EQ 0 AND j LT designcounter - 1 THEN
				BEGIN Outimage; Setpos (31); END;
			END;
			Outtext ("0  (in the original model)"); Outimage;
		END;
		IF (submodel GT 0 OR optionwanted [9]) AND NOT reduction THEN
		BEGIN	Eject (Line + 5); IF NOT savemodel
			THEN Outtext ("No reduction test possible, the original model was not saved")
			ELSE Outtext ("No increase in residual degrees of freedom or sum of squares");
			Outtext ("    (message)"); Outimage;
		END;
	END *** null_hypotheses ***;

	COMMENT ************************************* residual analysis **********************************;

	PROCEDURE residual_analysis;
	IF optionwanted [3] AND (submodel EQ 0 OR subselection GT 0) THEN
	INSPECT printstream DO
	BEGIN	INTEGER i, j, studnum; LONG REAL ARRAY resrow, resdev [1 : termcounter];
		LONG REAL observation, fittedvalue, residual, residualsum,
		standres, studdev, studres, maxstud, dfstud, frstud, studprob, sqrtmsres;
		resetdesigndata; pagenr:= 0; newpage2;
		studnum:= 0; residualsum:= 0; maxstud:= - Maxreal;
		sqrtmsres:= Sqrt (msres * (designsize - termcounter) / designsize);
		FOR i:= 2 STEP 1 UNTIL termcounter DO
		FOR j:= i - 1 STEP - 1 UNTIL 1 DO data [i,j]:= data [j,i];
		INSPECT outputstream DO BEGIN Outint (designsize, 12); Outimage; END;
		FOR i:= 1 STEP 1 UNTIL designsize DO
		BEGIN	FOR j:= 1 STEP 1 UNTIL termcounter DO
			resrow [j]:= designdata [j] . lastnumber;
			weightroot := data [i, weightcounter];
			observation:= data [i, designcounter] / weightroot;
			fittedvalue:= vecvec (1, termcounter, resrow, depvar) / weightroot;
			FOR j:= 1 STEP 1 UNTIL termcounter DO
			resdev [j]:= tamvec (1, termcounter, j, data, resrow);
			standdev:= Sqrt (vecvec (1, termcounter, resrow, resdev) * msres) / weightroot;
			residual:= observation - fittedvalue;
			residualsum:= residualsum + residual * weightroot * weightroot;
			standres:= IF perfectfit THEN 0 ELSE residual / sqrtmsres;
			studdev:= msres - standdev * standdev;
			studres:= IF studdev LE precision * msres THEN 0 ELSE residual / Sqrt (studdev);
			IF studres * studres GT maxstud THEN
			BEGIN studnum:= i; maxstud:= studres * studres; END;
			Outint (i, 6); Outfix (observation, 6, 20); Outfix (fittedvalue, 6, 20);
			Outfix (standdev, 6, 20); Outfix (residual, 6, 20);
			Outfix (standres, 6, 20); Outfix (studres, 6, 20); Outimage;
			IF Mod (i, 10) EQ 0 AND i LT designsize THEN
			BEGIN IF Mod (i, 50) EQ 0 THEN newpage2 ELSE Eject (Line + 1); END;
			INSPECT outputstream DO
			BEGIN	Outreal (observation, 16, 22); Outreal (fittedvalue, 16, 22);
				Outreal (standdev, 16, 22); Outreal (residual, 16, 22);
				Outreal (standres, 16, 22); Outreal (studres, 16, 22); Outimage;
			END;
		END;
		dfstud:= designsize - termcounter - 1; frstud:= dfstud + 1 - maxstud;
		studprob:= IF dfstud EQ 0 THEN 1 ELSE IF frstud LE precision THEN 0
		ELSE designsize * (1 - Fisher (maxstud * dfstud / frstud, 1, dfstud));
		Eject (Line + 1); Setpos (56);
		Outtext ("sum of residuals :"); Outfix (residualsum, 6, 13); Outimage;
		Eject (Line + 1); Setpos (20 - field (studnum));
		Outtext ("Upper bound for the right tail probability of ");
		Outtext ("the largest absolute studentized residual (no. ");
		Outint (studnum, field (studnum)); Outtext (") :"); Setpos (118);
		Outfix (Lmin (studprob, 1&&0), 6, 9); Outimage;
	END *** residual_analysis ***;

	COMMENT ************************************* process submodels **********************************;

	INTEGER subselection, submodel, notprocessed;

	PROCEDURE resetdesigndata;
	BEGIN	INTEGER j;
		FOR j:= 1 STEP 1 UNTIL termcounter DO designdata [j] . reset;
	END *** resetdesigndata ***;

	PROCEDURE process_submodels;
	BEGIN	IF optionwanted [5] THEN
		BEGIN	IF subselection EQ 0 THEN
			BEGIN	submodel:= submodel + 1;
				termcounter:= termcounter - 1;
				IF termcounter GE minterm THEN
				BEGIN resetdesigndata; GOTO process; END;
				notprocessed:= minterm - 1;
			END ELSE
			FOR programcounter:= programcounter + 1
			WHILE programcounter LE optioninstructions DO
			BEGIN	submodel:= program [programcounter];
				termcounter:= designcounter - 1 - submodel;
				IF termcounter GE minterm THEN
				BEGIN resetdesigndata; GOTO process; END;
				notprocessed:= notprocessed + 1;
			END;
		END;
		IF notprocessed GT 0 THEN INSPECT printstream DO
		BEGIN	Eject (Line + 5); Outtext ("Transformed submodels with as only independent");
			Outtext (" variable a constant, are not processed    (message)"); Outimage;
			Eject (Line + 5); Outtext ("Number of submodels not processed :");
			Outint (notprocessed, 3); Outimage;
		END;
	END *** process_submodels ***;

	COMMENT ************************************* regression analysis ********************************;

	PROCEDURE initialize_regression;
	BEGIN	INTEGER i, j;
		FOR i:= 1 STEP 1 UNTIL designsize DO
		BEGIN	FOR j:= IF submodel EQ 0 THEN weightcounter ELSE termcounter
			STEP - 1 UNTIL 1 DO data [i,j]:= designdata [j] . lastnumber;
			depvar [i]:= data [i, designcounter];
		END;
		IF submodel EQ 0 THEN
		FOR j:= 1 STEP 1 UNTIL designcounter DO
		resprod [j]:= tamvec (1, designsize, j, data, depvar);
	END *** initialize_regression ***;

	resetindex:= designdata [1] . index;
	matrixwidth:= (imagelength - 21) // 12;
	subselection:= optioninstructions - inputinstructions;
	programcounter:= inputinstructions; submodel:= notprocessed:= 0;

process : initialize_regression;	transformed_data_matrix;	control_information;
	  correlation_matrix;		no_regression_analysis;		least_squares;
	  multiple_correlation;		explained_variation;		standard_error_deviation;
	  regression_parameters;	covariance_matrix;		analysis_of_variance;
	  null_hypotheses;		residual_analysis;		process_submodels;

END *** regression ***;

END *** designdata ***;

END *** translator ***;

	COMMENT ************************************* final messages *************************************;

endjob : INSPECT outputstream DO
	 BEGIN Outtext ("""Eor"" of job :"); Outint (jobnumber, 3); Outimage; END;

	 IF printstream =/= Sysout THEN
	 BEGIN	IF erroneous THEN BEGIN Outtext ("Erroneous"); Outimage; END;
	 	Outtext ("End of job :"); Outint (jobnumber, 3); Outimage;
	 END;

	 INSPECT printstream DO
	 BEGIN	Eject (Line + 4); Outimage; Dotypeout (Sysout);
		Outtext ("End of job :"); Outint (jobnumber, 3); Outimage;
		IF printstream =/= Sysout THEN Eject (1); erroneous:= FALSE; GOTO job;
	 END;

  exit : Eject (Line + 1); Dotypeout (Sysout);
	 Outtext ("End of Multiple Linear Regression Analysis"); Outimage; closefiles;

END *** environment ***;

END *** inputsystem ***;