00100	TITLE FITS    -  A LEAST SQUARES CUBIC FIT  -  NOVEMBER 1972.
00200	
00300	COMMENT/
00400	
00500	   Y = (((B0*X + B1)*X + B2)*X + B3).
00600	
00700		  SY = B3*N   +  B2*SX +  B1*X2  +  B0*X3.
00800		  XY = B3*SX  +  B2*X2 +  B1*X3  +  B0*X4.
00900		 X2Y = B3*X2  +  B2*X3 +  B1*X4  +  B0*X5.
01000		 X3Y = B3*X3  +  B2*X4 +  B1*X5  +  B0*X6.
01100	/
01200	
01300	;COEFFICIENTS OF THE FOUR EQUATIONS.
01400	
01500		ROW1:	BLOCK 6
01600		ROW2:	BLOCK 6
01700		ROW3:	BLOCK 6
01800		ROW4:	BLOCK 6
01900	
02000	;Use Accumulators for accumulating.
02100		ACCUMULATORS{X,Y,SX,SY,X2,XY,X3,X2Y,X4,X3Y,X5,X6}
02200			I←1
02300	
02400	;CUBFIT(A,B,N)
02500	; - halfword array of N points (SX,y) given in A.
02600	; - returns four coefficients in real array B.
02700	
02800	SUBR(CUBFIT)
02900	
03000		dac 12,AC12
03100		lac I,ARG3
03200		lacn ARG1↔SOS
03300		hrlm 0,I		;loop count and pointer.
03400		lac[xwd 4,5]↔setz 4,↔blt 15	;clear accumulators.
     

00100	;ACCUMULATE PRODUCTS OF THE COORDINATES OF THE DATA POINTS.
00200	
00300	L1:	NIP X,(I)	↔	FSC X,233-12
00400		NAP Y,(I)	↔	FSC Y,233-12
00500	
00600		FADR SX,X
00700		FADR SY,Y
00800		FMPR Y,X	↔	FADR XY,Y
00900		FMPR Y,X	↔	FADR X2Y,Y
01000		FMPR Y,X	↔	FADR X3Y,Y
01100	
01200		LAC X
01300		FMPR X↔FADR X2,
01400		FMPR X↔FADR X3,
01500		FMPR X↔FADR X4,
01600		FMPR X↔FADR X5,
01700		FMPR X↔FADR X6,
01800	
01900	L2:	AOBJN I,L1
02000	
02100	;SETUP THE COEFFICIENTS.
02200	
02300		LAC 0,ARG1↔AOS↔FSC 0,233
02400		MOVEI 1,ROW1
02500		MOVEI 2,ROW2
02600		DAC  0,(1)1↔DAC SX,(1)2↔DAC X2,(1)3↔DAC X3,(1)4↔DAC SY,(1)5
02700		DAC SX,(2)1↔DAC X2,(2)2↔DAC X3,(2)3↔DAC X4,(2)4↔DAC XY,(2)5
02800		MOVEI 3,ROW3
02900		MOVEI 4,ROW4
03000		DAC X2,(3)1↔DAC X3,(3)2↔DAC X4,(3)3↔DAC X5,(3)4↔DAC X2Y,(3)5
03100		DAC X3,(4)1↔DAC X4,(4)2↔DAC X5,(4)3↔DAC X6,(4)4↔DAC X3Y,(4)5
     

00100	;Solve four equations in four unknowns by Gaussian Elimination.
00200	
00300	;TRIANGULARIZATION.
00400		;Accumulators 1,2,3,4 are row pointers.
00500		R←5	;Reciprocal of the Pivot.
00600		M←6	;Multiplier of the row.
00700		SETZM PARITY#
00800	
00900	;Find column-1 pivot for 4 by 4 matrix.
01000	
01100	T1:	lacm (1)1↔lacm M,(2)1↔caml M↔jrst .+3↔exch 1,2↔setcmm PARITY
01200		lacm (1)1↔lacm M,(3)1↔caml M↔jrst .+3↔exch 1,3↔setcmm PARITY
01300		lacm (1)1↔lacm M,(4)1↔caml M↔jrst .+3↔exch 1,4↔setcmm PARITY
01400	
01500	;Reciprocal of the Pivot.
01600	
01700		movsi R,(1.0)↔fdvr R,(1)1
01800	
01900	;Zero column-1 elements of rows 2, 3, & 4.
02000	
02100		lac M,(2)1↔fmpr M,R↔	 SETZM (2)1
02200			lacn M↔fmpr (1)2↔fadrm (2)2
02300			lacn M↔fmpr (1)3↔fadrm (2)3
02400			lacn M↔fmpr (1)4↔fadrm (2)4
02500			lacn M↔fmpr (1)5↔fadrm (2)5
02600	
02700		lac M,(3)1↔fmpr M,R↔	 SETZM (3)1
02800			lacn M↔fmpr (1)2↔fadrm (3)2
02900			lacn M↔fmpr (1)3↔fadrm (3)3
03000			lacn M↔fmpr (1)4↔fadrm (3)4
03100			lacn M↔fmpr (1)5↔fadrm (3)5
03200	
03300		lac M,(4)1↔fmpr M,R↔	 SETZM (4)1
03400			lacn M↔fmpr (1)2↔fadrm (4)2
03500			lacn M↔fmpr (1)3↔fadrm (4)3
03600			lacn M↔fmpr (1)4↔fadrm (4)4
03700			lacn M↔fmpr (1)5↔fadrm (4)5
03800	
03900	;Normalize First Row.
04000	
04100		FMPRM R,(1)1
04200		fmprm R,(1)2
04300		fmprm R,(1)3
04400		fmprm R,(1)4
04500		fmprm R,(1)5
     

00100	;Find column-2 pivot for 3 by 3 matriX
00200	T2:	lacm(2)2↔lacm M,(3)2↔caml M↔jrst .+3↔exch 2,3↔setcmm PARITY
00300		lacm(2)2↔lacm M,(4)2↔caml M↔jrst .+3↔exch 2,4↔setcmm PARITY
00400	
00500	;Reciprocal of the pivot.
00600	
00700		movsi R,(1.0)↔fdvr R,(2)2
00800	
00900	;Zero column-2 elements of rows 3 & 4.
01000	
01100		lac M,(3)2↔fmpr M,R	↔SETZM (3)2
01200			lacn M↔fmpr (2)3↔fadrm (3)3
01300			lacn M↔fmpr (2)4↔fadrm (3)4
01400			lacn M↔fmpr (2)5↔fadrm (3)5
01500	
01600		lac M,(4)2↔fmpr M,R	↔SETZM (4)2
01700			lacn M↔fmpr (2)3↔fadrm (4)3
01800			lacn M↔fmpr (2)4↔fadrm (4)4
01900			lacn M↔fmpr (2)5↔fadrm (4)5
02000	
02100	;Normalize Second Row.
02200	
02300		FMPRM R,(2)2
02400		fmprm R,(2)3
02500		fmprm R,(2)4
02600		fmprm R,(2)5
     

00100	;Find column-3 pivot for 2 by 2 matriX
00200	
00300	T3:	lacm(3)3↔lacm M,(4)3↔caml M↔jrst .+3↔exch 3,4↔setcmm PARITY
00400	
00500	;Reciprocal of the pivot.
00600	
00700		movsi R,(1.0)↔fdvr R,(3)3
00800	
00900	;Zero column-3 element of row 4.
01000	
01100		lac M,(4)3↔fmpr M,R	↔SETZM (4)3
01200			lacn M↔fmpr (3)4↔fadrm (4)4
01300			lacn M↔fmpr (3)5↔fadrm (4)5
01400	
01500	;Normalize Row-3.
01600	
01700		FMPRM R,(3)3
01800		fmprm R,(3)4
01900		fmprm R,(3)5
02000	
02100	;Find column-4 pivot for 1 by 1 matriX  (trivail).
02200	;Reciprocal of the pivot.
02300	
02400	T4:	movsi R,(1.0)↔fdvr R,(4)4
02500	
02600	;Normalize Row-4.
02700	
02800		FMPRM R,(4)4
02900		fmprm R,(4)5
     

00100	;Back Substitution.
00200		
00300	T5:	lacn (4)5↔fmpr (3)4↔fadrm (3)5
00400		lacn (4)5↔fmpr (2)4↔fadrm (2)5
00500		lacn (4)5↔fmpr (1)4↔fadrm (1)5
00600		lacn (3)5↔fmpr (2)3↔fadrm (2)5
00700		lacn (3)5↔fmpr (1)3↔fadrm (1)5
00800		lacn (2)5↔fmpr (1)2↔fadrm (1)5
00900	
01000	;Move results to Vector B.
01100		lac  5,ARG2	;B vector pointer.
01200		skipe PARITY↔jrst T7
01300	
01400	T6:	lac(1)5↔dac(5)3
01500		lac(2)5↔dac(5)2
01600		lac(3)5↔dac(5)1
01700		lac(4)5↔dac(5)0
01800		lac 12,AC12↔pop3j
01900		
02000	T7:	lacn(1)5↔dac(5)3
02100		lacn(2)5↔dac(5)2
02200		lacn(3)5↔dac(5)1
02300		lacn(4)5↔dac(5)0
02400		lac 12,AC12↔pop3j
02500	END