Trailing-Edge
-
PDP-10 Archives
-
decuslib10-10
-
43,50517/sqrt.mac
There is 1 other file named sqrt.mac in the archive. Click here to see a list.
TITLE SQRT for RPGLIB V1
SUBTTL Square root and floating point routines
; SQRT for RPGLIB V1
;
;
; This module contains the square root routine along with other
; miscellaneous floating point routines.
;
; Copyright (C) 1976, Bob Currier and Cerritos College
; All rights reserved
;
;
SEARCH RPGPRM, UUOSYM, MACTEN
%%LBLP==:%%LBLP
DEBUG==:DEBUG
BIS==:BIS
SALL
TWOSEG
RELOC 400000
ENTRY SQRT. ; square root routine
ENTRY FLOT.1 ; float a single precision integer
ENTRY FLOT.2 ; float a double precision integer
ENTRY FIX. ; fix a double precision integer
PR==15
NM==5
AC0==0
AC1==1
G==T2
H==T3
;SQRT. Main square root routine
;
;Uses Newton's method
;
;
SQRT.: JUMPL AC1,SQRT.1 ; must be positive
MOVE G,AC1 ; get arg
FDVR G,[2.0] ; get a guess
SQ.1: MOVE H,AC1 ; get arg
FDVR H,G ; get H/G
FAD H,G ; get G+(H/G)
FDVR H,[2.0] ; get (G+(H/G))/2
MOVE T1,H ; get it where we can play with it
FSB T1,G ; get difference
TLNN T1,377000 ; get it to nearest 1e-127
JRST SQ.2 ; yes -
MOVE G,H ; no get new guess
JRST SQ.1 ; and loop
SQ.2: MOVE AC1,H ; get result into proper AC
POPJ PP, ; and exit
SQRT.1: PUSHJ PP,%%H.11## ; square root of negative number illegal
SETZ AC1, ; return 0 on continue
POPJ PP, ; exit
;FIX. Convert floating point to double precision fixed point
;
;Calling sequence:
;
; MOVE AC,[floating.point.number]
; MOVE 16,[Z AC,result.address]
; PUSHJ 17,FIX.
;
;
FIX.: LDB T1,[POINT 4,PARM,12] ; get the AC
MOVE T1,(T1) ; get the floating point number
TDNN T1,[XWD 777,777777] ; is it zero?
JRST CC3C2D ; yes -
MOVM T2,T1 ; no -
MULI T2,400 ; get exponent
MOVEI T4,0 ; zap
ASHC T3,-306(T2) ; shift to get integer
JUMPGE T1,CC3C2C ; jump if input was positive
SETCA T3, ; else negate the two word result
MOVNS T4 ; complement and increment
SKIPN T4 ; zero?
ADDI T3,1 ; yes - increment
SKIPGE T3
TLO T4,1B18 ; set sign bit if necessary
CC3C2C: MOVEM T3,0(PARM) ; stash result
MOVEM T4,1(PARM) ; both parts
POPJ PP, ; exit
CC3C2D: SETZM 0(PARM) ; input is zero - give zero result
SETZM 1(PARM) ; double precision
POPJ PP, ; exit
;FLOT.1 Convert single precision fixed point to floating point
;
;Calling sequence:
;
; MOVE 16,[Z AC,operand]
; PUSHJ 17,FLOT.1
;
;Exit with result in accumulator specified by AC field of 16
;
;
FLOT.1: SKIPN NM,0(PARM) ; is input zero?
JRST FLOT1A ; yes -
IDIVI NM,1B18 ; no - break it into two pieces
CAIE NM,0 ; high part zero?
TLC NM,254000 ; no - jam exponent
TLC NM+1,233000 ; jam exponent of low part
FADR NM,NM+1 ; add two halves
FLOT1A: LDB T2,[POINT 4,PARM,12] ; get AC field
MOVEM NM,0(T2) ; store result
POPJ PP, ; and exit
;FLOT.2 Convert double precision fixed point number to floating point
;
;Calling sequence:
;
; MOVE 16,[Z AC,ADDRESS.OF.HIGH.ORDER.HALF.OF.FIXED.POINT.NUMBER]
; PUSHJ 17,FLOT.2
;
;Exit with result in AC specified by AC field of 16
;
;
FLOT.2: MOVE PR,[1.0] ; set scaling factor
MOVE T2,1(PARM) ; get low order
SKIPL T1,0(PARM) ; get high order part - is it negative?
JRST FLOT2B ; no -
TLO T2,1B18 ; yes - be sure low part is negative also
SETCA T1, ; negate the two word value
MOVMS T2 ; get magnitude
TLZ T2,1B18 ; reset sign bit
JUMPN T2,FLOT2C ; zero?
AOJA T1,FLOT2C ; yes - bump and jump
FLOT2B: TLZ T2,1B18 ; make sure low part is positive
FLOT2C: TLNE T1,377776 ; is value > 2**54 ?
JRST FLOT2E ; yes -
ASHC T1,10 ; no - shift over 8 bits
ASH T2,-10 ; restore low part
SKIPE T1 ; set exponent in high part unless
TLC T1,266000 ; it is zero
SKIPE T2 ; set exponent in low part unless
TLC T2,233000 ; it is negative
FADR T1,T2 ; add two halves
FMP T1,PR ; multiply by scaling factor
SKIPGE 0(PARM) ; was input negative?
MOVNS T1 ; yes - negate result
LDB T2,[POINT 4,PARM,12] ; get AC field
MOVEM T1,0(T2) ; store result
POPJ PP, ; exit
FLOT2E: ADD PR,[1B8] ; number > 2**54 - increase scaling factor
ASHC T1,-1 ; reduce T1
JRST FLOT2C ; try again
END