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