--- /dev/null
+/
+/
+/ 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
projects / pdp8.git / blobdiff