[pdp8.git] / sw / f4 / FRTSRC / sqrt.ra

/
/
/	S  Q  R  T
/	-  -  -  -
/
/SUBROUTINE	SQRT(X)
/
/ VERSION 5A 4-27-77  PT
/
	SECT	SQRT
	JA	#SQRT
	0			/THE MANTISSA ND EXPOENT DIDDLE AREAS.
	0
SQRTEX,	0
	0
SQRT13,	0
	0
	13			/PHONEY EXPONENT PATCH.
/
	EXTERN	#ARGER
SQRTM1,	TRAP4	#ARGER
	TEXT	+SQRT  +
SQRTXR,	SETX	XRSQRT
	SETB	BPSQRT
BPSQRT,	F 0.0
XRSQRT,	F 0.0
SQRT1,	F 0.0
SQRT2,	F 0.0
SQRT3,	F 0.0
F1SQRT,	F 1.
F2SQRT,	F 2.
	ORG	10*3+BPSQRT
	FNOP
	JA SQRTXR
	0
SQTRTN,	JA	.
SQRTS1,	0			/IF BETWEEN 1/4 & 1/2
	3200
	0
	0			/IF BETWEEN 1/2 & 1
	2240
	0
/
SQRTS2,	7777			/IF BETWEEN 1/4 & 1/2
	2327
	7772
	7777			/IF BETWEEN 1/2 & 1
	3300
	0
	BASE	0
#SQRT,	STARTD
	FLDA	10*3
	FSTA	SQTRTN
	FLDA	0
	SETX	XRSQRT
	SETB	BPSQRT
	BASE	BPSQRT
	LDX	1,1
	FSTA	BPSQRT
	FLDA%	BPSQRT,1  /ADDR OF X
	FSTA	BPSQRT
	STARTF
	FLDA%	BPSQRT 	/GET X
	JEQ	SQTRTN	/IF =0 JUST RTN
	JLT	SQRTM1	/IF <0 THEN ERROR
	FSTA	SQRTEX+1	/SAVE NUMBER AWAY FOR A SECOND.
	FLDA	SQRT13		/GET A RIGHT ADJUSTED 13 IN THE FAC.
	FSTA	SQRTEX-2	/STORE AWAY RIGHT AHEAD OF THE EXPONENT.
	FLDA	SQRTEX		/NOW RETREIVE THE EXPONENT AS HIGH ORDER WORD.
	ALN	0		/CHOP OFF CRAP.
	JEQ	SQRTSC		/IS IT EXACTLY ZERO? IF SO, SPECIAL CASE.
	FNORM			/NORMALIZE IT.
	FSUB	F1SQRT		/NOW SUBTRACT ONE FROM IT.
	FDIV	F2SQRT		/CHOP IT IN HALF NOW.
	FSTA	SQRT1		/AND SAVE 1/2 EXP IN A TEMP.
	ALN	0		/NOW FIX THE EXPONENT.
	FNORM			/AND NORMALIZE IT TO REMOVE UNDESIRABLE BITS.
	FSUB	SQRT1		/NOW SUBTRACT OFF EXTRANEOUS BITS.
	FMUL	F2SQRT		/EXPAND IT AGAIN [FAC =0 OR -1], OR 0 TO +1
	JGE	.+3		/MAKE SURE ITS POSITIVE.
	FNEG			/NOW MAKE IT 0 IF NO BIT OR +1 IF BIT
SQRTBK,	ATX	1		/SAVE IN AN INDEX.
	FSUB	F1SQRT		/SUBTRACT ONE TO MAKE IT -1 IF NO BIT OR 0 IF BIT.
	ALN	0		/AND NOW SHIFT IT RIGHT.
	FSTA	SQRTEX-1	/AND SAVE IT OVER THE OLD EXPONENT.
	FLDA	SQRT1		/RECALL OLD PART
	ALN	0		/FIX IT UP, NOW.
	FSTA	SQRT1		/AND STORE IT BACK FOR LATER USE
/
/		SQRTEX IS NOW 1/4 <X< 1
/
	FLDA	SQRTEX+1	/RECALL NUMBER.
	FSTA	SQRT2		/SAVE IN A TEMP.
/
	FMUL	SQRTS1,1	/MULTIPLY BY CORRECT CONSTANT.
	FADD	SQRTS2,1	/AND NOW ADD IN CORRECT CONSTANT.
/
/		NOTE: INITIAL APPROXIMATION DEPENDS ON WHETHER X IS 1/4<X<1/2 OR
/		     1/2<X<1
/
	FSTA	SQRT3		/SAVE IN A SECOND TEMP.
	FLDA	SQRT2		/RECALL INITIAL.
	FDIV	SQRT3		/CALCULATE X(0)/X(1)
	FADD	SQRT3		/X(1)+X(0)/X(1)
	FDIV	F2SQRT		/1/2(X(1)+X(0)/X(1))
	FSTA	SQRT3		/SAVE AGAIN. NOW X(2)
	FLDA	SQRT2		/RECALL ORIGINAL.
	FDIV	SQRT3		/X(0)/X(2)
	FADD	SQRT3		/X(2)+X(0)/X(2)
	FSTA	SQRTEX+1	/NOW STORE AWAY FOR FINAL EXPONENT DIDDLING.
/
	STARTD
/
	FCLA			/ZERO HIGH ORDER EXPONENT PART.
	FSTA	SQRTEX-1
	FLDA	SQRT1		/RECALL MODIFIED EXPONENT.
	FADDM	SQRTEX		/UPDATE FRACTIONAL EXPONENT.
/
	STARTF			/RETRUN TO FLOATING MODE.
/
	FLDA	SQRTEX+1	/PICK UP THE ANSWER.
	JA	SQTRTN		/AND RTN
/
SQRTSC,	FSUB	F1SQRT		/SPECIAL CASE FUDGE.
	FSTA	SQRT1		/SET EXPONENT ADD ON TO -1.
	FNEG			/AND SET ODD BIT ON.
	JA	SQRTBK		/AND GO BACK UP.
\f
Commit	Line	Data
7af5ad59 PH	1	/
	2	/
	3	/ S Q R T
	4	/ - - - -
	5	/
	6	/SUBROUTINE SQRT(X)
	7	/
	8	/ VERSION 5A 4-27-77 PT
	9	/
	10	SECT SQRT
	11	JA #SQRT
	12	0 /THE MANTISSA ND EXPOENT DIDDLE AREAS.
	13	0
	14	SQRTEX, 0
	15	0
	16	SQRT13, 0
	17	0
	18	13 /PHONEY EXPONENT PATCH.
	19	/
	20	EXTERN #ARGER
	21	SQRTM1, TRAP4 #ARGER
	22	TEXT +SQRT +
	23	SQRTXR, SETX XRSQRT
	24	SETB BPSQRT
	25	BPSQRT, F 0.0
	26	XRSQRT, F 0.0
	27	SQRT1, F 0.0
	28	SQRT2, F 0.0
	29	SQRT3, F 0.0
	30	F1SQRT, F 1.
	31	F2SQRT, F 2.
	32	ORG 10*3+BPSQRT
	33	FNOP
	34	JA SQRTXR
	35	0
	36	SQTRTN, JA .
	37	SQRTS1, 0 /IF BETWEEN 1/4 & 1/2
	38	3200
	39	0
	40	0 /IF BETWEEN 1/2 & 1
	41	2240
	42	0
	43	/
	44	SQRTS2, 7777 /IF BETWEEN 1/4 & 1/2
	45	2327
	46	7772
	47	7777 /IF BETWEEN 1/2 & 1
	48	3300
	49	0
	50	BASE 0
	51	#SQRT, STARTD
	52	FLDA 10*3
	53	FSTA SQTRTN
	54	FLDA 0
	55	SETX XRSQRT
	56	SETB BPSQRT
	57	BASE BPSQRT
	58	LDX 1,1
	59	FSTA BPSQRT
	60	FLDA% BPSQRT,1 /ADDR OF X
	61	FSTA BPSQRT
	62	STARTF
	63	FLDA% BPSQRT /GET X
	64	JEQ SQTRTN /IF =0 JUST RTN
65	JLT SQRTM1 /IF <0 THEN ERROR
66	FSTA SQRTEX+1 /SAVE NUMBER AWAY FOR A SECOND.
67	FLDA SQRT13 /GET A RIGHT ADJUSTED 13 IN THE FAC.
68	FSTA SQRTEX-2 /STORE AWAY RIGHT AHEAD OF THE EXPONENT.
69	FLDA SQRTEX /NOW RETREIVE THE EXPONENT AS HIGH ORDER WORD.
70	ALN 0 /CHOP OFF CRAP.
71	JEQ SQRTSC /IS IT EXACTLY ZERO? IF SO, SPECIAL CASE.
72	FNORM /NORMALIZE IT.
73	FSUB F1SQRT /NOW SUBTRACT ONE FROM IT.
74	FDIV F2SQRT /CHOP IT IN HALF NOW.
75	FSTA SQRT1 /AND SAVE 1/2 EXP IN A TEMP.
76	ALN 0 /NOW FIX THE EXPONENT.
77	FNORM /AND NORMALIZE IT TO REMOVE UNDESIRABLE BITS.
78	FSUB SQRT1 /NOW SUBTRACT OFF EXTRANEOUS BITS.
79	FMUL F2SQRT /EXPAND IT AGAIN [FAC =0 OR -1], OR 0 TO +1
80	JGE .+3 /MAKE SURE ITS POSITIVE.
81	FNEG /NOW MAKE IT 0 IF NO BIT OR +1 IF BIT
82	SQRTBK, ATX 1 /SAVE IN AN INDEX.
83	FSUB F1SQRT /SUBTRACT ONE TO MAKE IT -1 IF NO BIT OR 0 IF BIT.
84	ALN 0 /AND NOW SHIFT IT RIGHT.
85	FSTA SQRTEX-1 /AND SAVE IT OVER THE OLD EXPONENT.
86	FLDA SQRT1 /RECALL OLD PART
87	ALN 0 /FIX IT UP, NOW.
88	FSTA SQRT1 /AND STORE IT BACK FOR LATER USE
89	/
90	/ SQRTEX IS NOW 1/4 <X< 1
91	/
92	FLDA SQRTEX+1 /RECALL NUMBER.
93	FSTA SQRT2 /SAVE IN A TEMP.
94	/
95	FMUL SQRTS1,1 /MULTIPLY BY CORRECT CONSTANT.
96	FADD SQRTS2,1 /AND NOW ADD IN CORRECT CONSTANT.
97	/
98	/ NOTE: INITIAL APPROXIMATION DEPENDS ON WHETHER X IS 1/4<X<1/2 OR
99	/ 1/2<X<1
100	/
101	FSTA SQRT3 /SAVE IN A SECOND TEMP.
102	FLDA SQRT2 /RECALL INITIAL.
103	FDIV SQRT3 /CALCULATE X(0)/X(1)
104	FADD SQRT3 /X(1)+X(0)/X(1)
105	FDIV F2SQRT /1/2(X(1)+X(0)/X(1))
106	FSTA SQRT3 /SAVE AGAIN. NOW X(2)
107	FLDA SQRT2 /RECALL ORIGINAL.
108	FDIV SQRT3 /X(0)/X(2)
109	FADD SQRT3 /X(2)+X(0)/X(2)
110	FSTA SQRTEX+1 /NOW STORE AWAY FOR FINAL EXPONENT DIDDLING.
111	/
112	STARTD
113	/
114	FCLA /ZERO HIGH ORDER EXPONENT PART.
115	FSTA SQRTEX-1
116	FLDA SQRT1 /RECALL MODIFIED EXPONENT.
117	FADDM SQRTEX /UPDATE FRACTIONAL EXPONENT.
118	/
119	STARTF /RETRUN TO FLOATING MODE.
120	/
121	FLDA SQRTEX+1 /PICK UP THE ANSWER.
122	JA SQTRTN /AND RTN
123	/
124	SQRTSC, FSUB F1SQRT /SPECIAL CASE FUDGE.
125	FSTA SQRT1 /SET EXPONENT ADD ON TO -1.
126	FNEG /AND SET ODD BIT ON.
127	JA SQRTBK /AND GO BACK UP.
128	\f