Google
 

Trailing-Edge - PDP-10 Archives - fortv11 - mthgdb.mac
There are 9 other files named mthgdb.mac in the archive. Click here to see a list.

	SEARCH MTHPRM
	TV	MTHGDB	GFLOAT DOUBLE PRECISION ROUTINES ,2(4020)

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SUBTTL	REVISION HISTORY

COMMENT \

***** Begin Revision History *****

1100	CKS	--------
	Creation.

1464	DAW	-------
	Put LERR's in new format.

1673	CDM	9-Sept-81
	Added functions (GDIM, GINT, GNINT, GPROD IGNINT).

1721	CDM	15-Sept-81
	Added fix to GDIM to prevent overflow for GDIM(-INF,+INF).

1724	BL	17-Sep-81
	Clean up error messages.

***** Begin Mathlib *****

3055	CKS/JLC	17-Mar-82
	Fix underflow/overflow error message in GEXP2.

3200	JLC	10-May-82
	Mathlib integration. Change all LERRs to $LCALLs. Change
	TWOSEG/RELOC to SEGMENT macros.

3201	RVM	20-May-82
	Fix problem that made an argument of plus or minus one
	to GASIN and GACOS illegal.

3206	AHM	8-Jun-82
	Remove CKS's older version  of GINT. in  favor of CDM's  newer
	version.

3215	RVM	20-Aug-82
	Change a CAIGE to a CAIG in GINT.  Numbers less than 1.00D0
	were not being set to 0.00D0.

***** Begin Version 1A *****

3233	CKS	23-Mar-83
	Rewrote GMOD to remove restrictions on the magnitudes of
	its arguments.

3243	BCM	29-Mar-83
	Get and clear PC flags with new GETFLG and RESFLG macros.  Fixes
	always clearing underflow in MTHTRP.  Only applies to JOV branches.

***** End Revision History *****

***** Begin Version 2 *****

4002	JLC	3-Aug-83
	Move GEXP. to after GTAN., which calls GEXP.
	Make 1.0 a special case for GEXP3.

4011	JLC	4-May-84
	Add code for MOD, AMOD, DMOD, and GMOD for 2nd arg=0.

4012	MRB	6-Jun-84
	Change name of routines GTODA and DTOGA to GTODA. and DTOGA.
	to allow flagger. These symbols will be defined in FORCOM.

4145	MRB	17-Aug-84
	Change GABS., GMAX1., GMIN1. & GSIGN. to <nnn.+0> to get 
	macro to generate the correct polish stuff.

4020	DCE	17-Dec-85
	The GFLOATING exponentiation routine uses an incorrect method of
	getting the exponent of the magnitude of certain intermediate 
	GFLOATING numbers. A MOVM on the first word of the double word 
	value is used followed by stripping out the exponent field. 
	This does not work for all negative GFLOATING numbers where the 
	mantissa part in the first word is all zero (the first 
	non-zero bit of the mantissa is in the second word). Due to the 
	fact that these errors occur in an intermediate calculation, 
	it is impossible to specify the class of input numbers which 
	will result in this problem.
\
	PRGEND
TITLE	GACOS	ARC SINE AND ARC COSINE FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 19, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GACOS
EXTERN	GACOS.
GACOS=GACOS.
PRGEND
TITLE	GASIN	ARC SINE AND ARC COSINE FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 19, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GASIN
EXTERN	GASIN.
GASIN=GASIN.
PRGEND
TITLE	GASIN.	ARC SINE AND ARC COSINE FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JUNE 19, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DASIN(X) AND DACOS(X) ARE CALCULATED AS FOLLOWS

;  LET R(G) = (G*(RP1+G*(RP2+G*(RP3+G*(RP4+G*RP5)))))/(Q0+G*(Q1
;               +G*(Q2+G*(Q3+G*(Q4+G)))))
;      (RP1,RP2,RP3,RP4,RP5,Q0,Q1,Q2,Q3,Q4, AND Q5 ARE GIVEN BELOW)

;  LET S = Y + Y*R(G)

;  FOR SITUATIONS 1. AND 3. BELOW,
;	G = Y**2, WHERE Y = ABS(X)

;  FOR SITUATIONS 2. AND 4., G = (1.0 - ABS(X))/2.0
;	AND Y = -2.0*SQRT(G)

;  LET W = ABS(X) AND TERM THE RESULT T.
;  THEN, CONSIDER THE FOUR SITUATIONS:

;	1.  DASIN, FOR W <= 0.5.  T = S, AND T IS NEGATED 
;		FOR NEGATIVE X
;	2.  DASIN, FOR W > 0.5.  T = PI/2 + S, AND T IS NEGATED
;		FOR NEGATIVE X
;	3.  DACOS, FOR W <= 0.5.  T = PI/2 - S IF X > 0
;				    = PI/2 + S IF X < 0
;	4.  DACOS FOR W > 0.5.  T = -S IF X > 0
;				  = PI + S IF X < 0.

;THE RANGE OF DEFINITION FOR DASIN/DACOS IS ABS(X) <= 1.0. AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE.  DASIN/DACOS WILL BE SET
;  TO + MACHINE INFINITY.

;REQUIRED (CALLED) ROUTINES:  DSQRT

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GASIN
;	OR
;  PUSHJ	P,GACOS

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GACOS,.)	;ENTRY TO GACOS ROUTINE
	HRRZI	T0,2		;SET FLAG TO TWO
	MOVEM	T0,FLAG		;MOVE FLAG TO MEMORY
	JRST	GETX		;GO TO GETX
	
	HELLO	(GASIN,.)	;ENTRY TO GASIN ROUTINE
	SETZM	FLAG		;SET FLAG TO ZEROES
GETX:	DMOVE	T0,@(L)		;OBTAIN X
	JUMPGE	T0,XPOS		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;Y = -X
XPOS:	CAMGE	T0,HALF		;IF Y IS .LT. 1/2
	  JRST	LE		;GO TO LE
	CAME	T0,HALF		;IF Y IS .GT. 1/2
	  JRST	GT		;GO TO GT
	JUMPE	T1,LE		;
GT:	CAMLE	T0,ONE		;IF Y IS .GT. ONE
	  JRST	MSTK		;GO TO MSTK
	CAME	T0,ONE		;IF Y IS .LT. ONE
	  JRST	ALG		;GO TO ALG
	JUMPE	T1,ALG		;[3201] CHECK SECOND WORD
MSTK:	$LCALL	AOI
;LERR	(LIB,%,<DASIN or DACOS: ABS(arg) .GT. 1.0; result = +infinity>)
	HRLOI	T0,377777	;RESULT = 
	SETO	T1,		;+MACHINE INFINITY
	GOODBY	(1)		;RETURN

ALG:	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	MOVN	T5,FLAG		;GET A COPY OF - FLAG
	HRRZI	T4,2		;SET T4 TO 2
	ADD	T5,T4		;I = 2-FLAG
	HRLZI	T2,200140	;SET T2,T3 TO
	SETZ	T3,		;ONE
	GFSB	T2,T0		;G = 1-Y
	EXTEND	T2,[GFSC -1]	;* 1/2
	DMOVEM	T2,TEMP		;SAVE A COPY OF G
	FUNCT	GSQRT.,<TEMP>	;Y = GSQRT(G)
	EXTEND	T0,[GFSC 1]	;* 2
	DMOVN	T0,T0		;Y = -Y
	MOVEM	T5,I		;MOVE I TO MEMORY
	JRST	GOTY		;GO TO GOTY
LE:	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	MOVE	T5,FLAG		;I = FLAG
	MOVEM	T5,I		;MOVE I TO MEMORY
	CAMGE	T0,TWOM30	;IF Y < TWOM30
	JRST	GOTRES		;GO TO GOTRES
	DMOVE	T2,T0		;SAVE A COPY OF Y
	GFMP	T2,T2		;G = Y*Y
GOTY:	DMOVE	T4,T2		;SAVE A COPY OF G
	GFAD	T4,Q4		;XDEN = G + Q4
	GFMP	T4,T2		;*G
	GFAD	T4,Q3		;+Q3
	GFMP	T4,T2		;*G
	GFAD	T4,Q2		;+Q2
	GFMP	T4,T2		;*G
	GFAD	T4,Q1		;+Q1
	GFMP	T4,T2		;*G
	GFAD	T4,Q0		;+Q0
	DMOVEM	T2,TEMP		;SAVE A COPY OF G
	GFMP	T2,RP5		;XNUM = G*RP5
	GFAD	T2,RP4		;+RP4
	GFMP	T2,TEMP		;*G
	GFAD	T2,RP3		;+RP3
	GFMP	T2,TEMP		;*G
	GFAD	T2,RP2		;+RP2
	GFMP	T2,TEMP		;*G
	GFAD	T2,RP1		;+RP1
	GFMP	T2,TEMP		;*G
	GFDV	T2,T4		;RESULT = XNUM/XDEN
	GFMP	T2,T0		;*Y
	JRST	GETI		;GO TO GETI
GOTRES:	DMOVE	T2,AHI		;ZERO T2
GETI:	MOVE	T5,I		;GET A COPY OF I
	SKIPE	FLAG		;IF FLAG IS .NE. 0
	  JRST	NE		;GO TO NE
	GFAD	T0,ALO(T5)
	GFAD	T0,T2
	GFAD	T0,AHI(T5)
	SKIPGE	@(L)		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;RESULT = -RESULT
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN

NE:	SKIPGE	@(L)		;IF X IS NEGATIVE
	  JRST 	NEGX		;GO TO NEGX
	DMOVN	T0,T0		;RESULT = -RESULT
	GFAD	T0,ALO(T5)
	GFSB	T0,T2
	GFAD	T0,AHI(T5)
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN

NEGX:	GFAD	T0,BLO(T5)
	GFAD	T0,T2
	GFAD	T0,BHI(T5)
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN

RP1:    DOUBLE  577211206522,276137263562       ;-.27368494524164255994D+2      
RP2:    DOUBLE  200671152471,260255221215       ;.57208227877891731407D+2       
RP3:    DOUBLE  577130237232,262622254077       ;-.39688862997504877339D+2      
RP4:    DOUBLE  200450470273,047303444713       ;.10152522233806463645D+2       
RP5:    DOUBLE  577723321022,120636063337       ;-.69674573447350646411D+0      
Q0:     DOUBLE  576726744776,016507406363       ;-.16421096714498560795D+3      
Q1:     DOUBLE  201164111170,200753410270       ;.41714430248260412556D+3       
Q2:     DOUBLE  576620210610,035225367030       ;-.38186303361750149284D+3      
Q3:     DOUBLE  201045571744,262621615746       ;.15095270841030604719D+3       
Q4:     DOUBLE  577220264274,210147474061       ;-.23823859153670238830D+2      
AHI:	DOUBLE	0,0
A1HI:   DOUBLE  200162207732,242102643021       ;PI/2                           
BHI:    DOUBLE  200262207732,242102643021       ;PI                             
B1HI:   DOUBLE  200162207732,242102643021       ;PI/2                           
ALO:	DOUBLE	0,0				;
A1LO:	DOUBLE	170551423063,024270033407	;NEXT 59 BITS OF PI/2
BLO:	DOUBLE	170651423063,024270033407	;NEXT 59 BITS OF PI
B1LO:	DOUBLE	170551423063,024270033407	;NEXT 59 BITS OF PI/2
TWOM30:	174340000000				;HIGH ORDER PART OF 2**(-30)
ONE:	200140000000				;1.0
HALF:	200040000000				;1/2

	SEGMENT	DATA

I:	0
FLAG:	0
TEMP:	DOUBLE	0,0
	PRGEND
TITLE	GATAN2	ARC TAN FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 9, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GATAN2
EXTERN	GATN2.
GATAN2=GATN2.
PRGEND
TITLE	GATAN	ARC TAN FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	Chris Smith/CKS

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GATAN
EXTERN	GATAN.
GATAN=GATAN.
PRGEND
TITLE	GATAN.	ARC TAN FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.



;DATAN(X) is computed as follows:
;
;If X < 0, compute DATAN(|X|) below, then DATAN(X) = -DATAN(-X).
;
;If X > 0, use the identity
;
;	DATAN(X) = DATAN(XHI) + DATAN(Z)
;	Z = (X - XHI) / (1 + X*XHI)
;
;where XHI is chosen so that  |Z| <= tan(pi/32).
;
;XHI is chosen to be exactly representable as a double precision number,
;and so Z can be calculated without loss of significance.
;
;DATAN(XHI) is found by table lookup.  It is stored as ATANHI + ATANLO to
;provide guard bits for the final addition to DATAN(Z).
;
;DATAN(Z) is evaluated by means of a polynomial approximation from Hart et al.
;(formula 4904).
;
;If X < tan(pi/32), DATAN(X) = DATAN(Z).
;If X > 1/tan(pi/32), DATAN(X) = pi/2 - DATAN(1/X).
;
;If tan(pi/32) < X < 1/tan(pi/32), an appropriate XHI is obtained by indexing
;into a table.  The table tells which XHI to use for various ranges of X.
;The index into the table is formed from the low 3 exponent bits and the high
;3 fraction bits of X.


	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	T6=6			;MORE AC NAMES
	T7=7
	T8=10



	HELLO (GATAN,.)		;DATAN ENTRY

	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	PUSH	P,T6

	DMOVE	T0,@(L)		;GET ARGUMENT X
	MOVEM	T0,SGNFLG	;SAVE ARGUMENT SIGN FOR RESULT
	CAIGE	T0,0		;GET |X|
	  DMOVN	T0,T0

	CAML	T0,MAXX		;IS X LARGE?
	  JRST	LARGEX		;YES, GO USE ATAN(X) = PI/2 - ATAN(1/X)
	SETZ	T2,		;T2 WILL GET OFFSET INTO XHI TABLES
	CAMG	T0,MINX		;IS X SMALL ENOUGH THAT NO ARG REDUCTION IS REQUIRED?
	  JRST	CALC		;YES, JOIN CALCULATION BELOW

	MOVE	T3,T0		;GET HIGH WORD OF X
	LSHC	T2,^D12		;GET EXPONENT, SHIFT HIGH FRACTION BIT
				;(ALWAYS 1) INTO SIGN BIT OF T3
	ASHC	T2,3		;GET THREE FRACTION BITS, LEAVING
				;THE 1 BEHIND
	HLRZ	T2,OFFSET-20000+24(T2) ;GET OFFSET INTO XHI TABLES

	DMOVE	T3,T0		;GET A COPY OF X
	GFSB	T0,XHI(T2)	;GET X-XHI
	GFMP	T3,XHI(T2)	;    X*XHI
	GFAD	T3,ONE		;    1 + X*XHI
	GFDV	T0,T3		;    (X-XHI) / (1 + X*XHI)

;Here T0 has the reduced argument Z with |Z| <= tan(pi/32).  T2 has the
;index into the ATAN(XHI) tables ATANHI and ATANLO.  SGNFLG is set negative
;if the result should be negated.

CALC:	DMOVE	T3,T0		;GET |Z|
	CAIGE	T3,0
	  DMOVN	T3,T3
	CAMG	T3,EPS		;IS IT SMALL ENOUGH THAT ATAN(Z) = Z?
	  JRST	SMALLX		;YES, BE FAST, AVOID UNDERFLOW
	GFMP	T3,T3		;GET Z**2
	DMOVE	T5,P06		;GET P(Z**2)
	GFMP	T5,T3
	GFAD	T5,P05
	GFMP	T5,T3
	GFAD	T5,P04
	GFMP	T5,T3
	GFAD	T5,P03
	GFMP	T5,T3
	GFAD	T5,P02
	GFMP	T5,T3
	GFAD	T5,P01
	GFMP	T3,T5
	GFMP	T3,T0		; * Z
	GFAD	T0,T3		; + Z = ATAN(Z)

SMALLX:	GFAD	T0,ATANLO(T2)	;  + ATAN(XHI) LOW
	GFAD	T0,ATANHI(T2)	;  + ATAN(XHI) HI   = DATAN(X)
	SKIPG	SGNFLG		;ATTACH SIGN TO RESULT
	  DMOVN	T0,T0
RET:	POP	P,T6		;RETURN
	POP	P,T5
	POP	P,T4
	POP	P,T3
	POP	P,T2
	POPJ	P,

LARGEX:	DMOVE	T3,T0		;GET -1/X
	DMOVN	T0,ONE
	GFDV	T0,T3
	MOVEI	T2,PI2OFFS	;GET OFFSET OF PI/2
	JRST	CALC		;GO COMPUTE PI/2 + ATAN(-1/X)
SUBTTL DATAN2 (Y,X)


;To compute DATAN2(Y,X), let U = |Y| and V = |X|, and compute DATAN(U/V).
;Then find DATAN2(Y,X) based on the signs of Y and X as follows:
;
;	 X	 Y	 DATAN2(Y/X)
;	
;	pos	pos	  DATAN(U/V)
;	pos	neg	 -DATAN(U/V)
;	neg	pos	-(DATAN(U/V) - pi)
;	neg	neg	  DATAN(U/V) - pi
;
;The add of -pi is combined with the add of DATAN(XHI) which is the last step
;of the DATAN algorithm.
;
;The reduced argument for DATAN is Z = (U/V - XHI) / (1 + U/V * XHI).
;This is rewritten as (U - V*XHI) / (V + U*XHI).  To accurately calculate the
;numerator, find VHI and VLO with 
;
;	V = VHI + VLO 
;	VHI has at most 27 significant bits
;	VLO has at most 35 significant bits
;
;and choose XHI with at most 13 significant bits.  Then VHI*XHI and VLO*XHI can
;be exactly represented as double precision numbers, and the numerator is
;
;	U - V*XHI = (U - VHI*XHI) - VLO*XHI



	HELLO (GATN2.,GATAN2)	;DATAN2 ENTRY

	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	PUSH	P,T6

	DMOVE	T0,@0(L)	;GET Y
	GFDV	T0,@1(L)	;GET Y/X
	  JFCL	EXCEP		;OVERFLOW AND UNDERFLOW CAN OCCUR
	DMOVEM	T0,SGNFLG	;RESULT SHOULD BE MULTIPLIED BY SGN(Y/X)

	DMOVE	T3,T0		;GET |Y/X|, ATAN ARG
	CAIGE	T3,0
	  DMOVN	T3,T3
	CAML	T3,MAXX		;IS |Y/X| LARGE?
	  JRST	LARGE2		;YES, GO USE ATAN(Y/X) = PI/2 - ATAN(X/Y)
	SETZ	T2,		;GET OFFSET INTO ATAN TABLES
	CAMG	T3,MINX		;IS |Y/X| SMALL ENOUGH TO USE POLYNOMIAL DIRECTLY?
	  JRST	[DMOVE T0,T3	;YES, DO SO
		 JRST SMALL2]
	LSHC	T2,^D12		;GET INDEX INTO OFFSET TABLES
	ASHC	T2,3
	HLRZ	T2,OFFSET-20000+24(T2) ;GET INDEX INTO XHI TABLES

	PUSH	P,T7		;SAVE MORE REGISTERS
	PUSH	P,T8

	DMOVE	T0,@0(L)	;GET |Y| = U
	CAIGE	T0,0
	  DMOVN	T0,0
	DMOVEM	T0,USAVE	;SAVE FOR LATER
	DMOVE	T3,@1(L)	;GET |X| = V
	CAIGE	T3,0
	  DMOVN	T3,T3
	MOVE	T5,T3		;GET A COPY OF V
	DMOVE	T7,T3		;GET ANOTHER
	SETZ	T6,		;GET HIGH 27 BITS OF V = VHI
	GFSB	T7,T5		;GET LOW 35 BITS OF V = VLO
	GFMP	T5,XHI(T2)	;GET V*XHI = VHI * XHI
	GFMP	T7,XHI(T2)	;	    + VLO * XHI
	GFSB	T0,T5		;GET (U - VHI*XHI)
	GFSB	T0,T7		;		   - VLO*XHI
	DMOVE	T5,USAVE	;GET U BACK
	GFMP	T5,XHI(T2)	;GET U * XHI
	GFAD	T3,T5		;GET V + U*XHI
	GFDV	T0,T3		;GET (U - V*XHI) / (V + U*XHI)

	POP	P,T8		;RESTORE REGISTERS
	POP	P,T7

SMALL2:	SKIPGE	@1(L)		;IF SECOND ARG (X) IS NEGATIVE
	  ADDI	T2,MPIOFFS	;  ADD -PI TO RESULT
	JRST	CALC		;GO GET ATAN AND RETURN

LARGE2:	DMOVN	T0,@1(L)	;GET -X/Y
	GFDV	T0,@0(L)
	MOVMM	T0,SGNFLG	;SET SGNFLG POSITIVE
	MOVEI	T2,PI2OFFS	;ADD PI/2 TO RESULT IF FIRST ARG (Y) IS POSITIVE
	SKIPGE	@0(L)
	  MOVEI	T2,MPI2OFFS	;ADD -PI/2 TO RESULT IF FIRST ARG (Y) IS NEGATIVE
	JRST	CALC		;GO COMPUTE +/- PI/2 + ATAN(-X/Y)

EXCEP:	SKIPN	@1(L)		;CHECK FOR DIVIDE BY 0
	  JRST	DIVCHK		;IF DIVIDE CHECK, GO CHECK FOR ATAN2(0,0)
	JUMPN	T0,OVER		;IF OVERFLOW, RESULT IS +/- PI/2
	SKIPL	@1(L)		;IF UNDERFLOW, CHECK SECOND ARGUMENT
	 $LCALL	RUN,RET
;	  LERR	(LIB,%,<DATAN2: result underflow>,,RET)
				;IF SECOND ARG (X) POSITIVE, RESULT UNDERFLOWS
	DMOVN	T0,MPI		;ELSE RESULT IS PI WITH SIGN OF FIRST ARG
	JRST	YSIGN		;GO ATTACH SIGN AND RETURN

DIVCHK:	SKIPE	@0(L)		;CHECK FOR BOTH ARGS 0
	  JRST	OVER		;ATAN2(NONZERO,0) IS SAME AS OVERFLOW
	$LCALL	BAZ
;LERR	(LIB,%,<DATAN2: both arguments are zero, result=0.0>)
	SETZB	T0,T1		;RETURN 0
	JRST	RET

OVER:	DMOVE	T0,PI2		;OVERFLOW, RESULT IS PI/2 WITH SIGN OF FIRST ARG
YSIGN:	SKIPGE	@0(L)		;ATTACH SIGN OF FIRST ARGUMENT (Y)
	  DMOVN	T0,T0
	JRST	RET		;RETURN
SUBTTL TABLES

;This table is indexed by the low 3 exponent bits and the high 3 fraction bits
;of X, where MINX < X < MAXX.  It gives the offset into XHI, ATANHI, and ATANLO
;of a suitable XHI.

DEFINE OFFS (X) <
 IRP X,<XWD 2*X,X>
>

	RADIX 10
OFFSET:	OFFS	<1,1,1,2,2,3,3,4,4,4,5,5,5,6,6,7,7,7,8,8,8,9,9,9>
	OFFS	<10,10,10,10,11,11,11,12,12,12,12,12,13,13,13,13>
	OFFS	<13,13,14,14,14,14,14,14,14,14,14,14,14,14,14>
	RADIX 8

PI2OFFS=2*^D15			;OFFSET INTO ATANHI AND ATANLO OF PI/2
MPIOFFS=2*^D16			;OFFSET OF -PI
MPI2OFFS=2*^D31			;OFFSET OF -PI/2

XHI: 	EXP    000000000000,0		; .0000000000000000000    not used
 	EXP    177565774000,0		; .1054534912109375000    X1
 	EXP    177640774000,0		; .1288757324218750000    X2
 	EXP    177647774000,0		; .1562194824218750000    
 	EXP    177661764000,0		; .1952209472656250000    
 	EXP    177677740000,0		; .2497558593750000000    
 	EXP    177747754000,0		; .3121948242187500000    
 	EXP    177761720000,0		; .3898925781250000000    
 	EXP    177777630000,0		; .4984130859375000000    
 	EXP    200051574000,0		; .6522216796875000000    
 	EXP    200067404000,0		; .8673095703125000000    
 	EXP    200145344000,0		; 1.170166015625000000    
 	EXP    200164510000,0		; 1.645019531250000000    
 	EXP    200251104000,0		; 2.570800781250000000    
 	EXP    200352700000,0		; 5.359375000000000000    X14

ATANHI:	EXP  000000000000,000000000000		; .0000000000000000000    ATAN(0)
 	EXP  177565626152,025171771712		; .1050651824695016828    ATAN(X1)
 	EXP  177640637315,230510354743		; .1281692623534198424    
 	EXP  177647527650,022607523014		; .1549669515088296697    
 	EXP  177661266130,275334150302		; .1927961221331547523    
 	EXP  177676517562,252343141216		; .2447488705187413410    
 	EXP  177746567507,360775545336		; .3026068193134274480    
 	EXP  177757453662,234657557072		; .3717628303965550499    
 	EXP  177773136266,273706162457		; .4623772720674932798    
 	EXP  200044771623,334166732152		; .5779354489017924541    
 	EXP  200055563263,207236241201		; .7144577222199199901    
 	EXP  200067214043,155315244501		; .8636495725719767220    
 	EXP  200140622720,225243500416		; 1.024591515455200379    
 	EXP  200146311711,011245646124		; 1.199822549393342833    
 	EXP  200154271467,254276104335		; 1.386328658261967262    ATAN(X14)
PI2: 	EXP  200162207732,242102643022		; 1.570796326794896620    PI/2
MPI: 	EXP  577515570045,135675134756		;-3.141592653589793241    -PI
 	EXP  577517324610,256420774615		;-3.036527471120291553    -PI + ATAN(X1)
 	EXP  577517622022,067321553615		;-3.013423391236373393    -PI + ATAN(X2)
 	EXP  577520155437,337025522117		;-2.986625702080963569    
 	EXP  577520643352,351552743372		;-2.948796531456638489    
 	EXP  577521515034,210413303027		;-2.896843783071051899    
 	EXP  577522447016,133774711512		;-2.838985834276365791    
 	EXP  577523535433,261363112666		;-2.769829823193238186    
 	EXP  577525103674,065265753224		;-2.679215381522299960    
 	EXP  577526766412,124732723411		;-2.563657204688000783    
 	EXP  577531124722,077544605217		;-2.427134931369873246    
 	EXP  577533433056,071160406077		;-2.277943081017816514    
 	EXP  577536101415,250416775165		;-2.117001138134592862    
 	EXP  577601672023,305040140060		;-1.941770104196450408    
 	EXP  577607651602,150070376271		;-1.755263995327825979    -PI + ATAN(X14)
 	EXP  577615570045,135675134756		;-1.570796326794896620    -PI/2

ATANLO:	EXP  000000000000,000000000000		; .0000000000000000000    
 	EXP  607611362655,250215702700		;-.9237013676958846587E-19
 	EXP  170143152717,266776243666		; .5964602757438210869E-19
 	EXP  170154077372,225515747423		; .7474896823762959724E-19
 	EXP  607735160270,370510376242		;-.2946026702232522009E-19
 	EXP  610004420274,014224300071		;-.2518569148307390217E-19
 	EXP  170347024456,077470314101		; .2645467902793737610E-18
 	EXP  607510312505,015573123116		;-.1883944550948088893E-18
 	EXP  170254520527,135456552170		; .1513056980995441703E-18
 	EXP  170452555761,155354735174		; .5788933265131145596E-18
 	EXP  170453236423,212017737506		; .5869551379299690115E-18
 	EXP  607321026727,320125051737		;-.6363620490158901181E-18
 	EXP  170443623615,021735205127		; .4850262996822095826E-18
 	EXP  607222115425,120272454763		;-.1242727478712024599E-17
 	EXP  607226174644,043107635453		;-.1131804334593366021E-17
 	EXP  607223046146,050560067016		;-.1217705177797396539E-17
 	EXP  170654731631,327217710762		; .2435410355594793078E-17
 	EXP  607123161307,104424247040		;-.2427449340111014898E-17
 	EXP  607106015110,124777656057		;-.3142794913755447875E-17
 	EXP  170471157106,257651240651		; .7754358478556155783E-18
 	EXP  170674303434,273125014744		; .3273311826560871401E-17
 	EXP  170570727666,236602064744		; .1542862926123315625E-17
 	EXP  170543470577,076355704763		; .9652336698973597403E-18
 	EXP  607111346356,050107456126		;-.2957154527430437100E-17
 	EXP  170577335536,232205477162		; .1719354315705933698E-17
 	EXP  607336325130,312454401102		;-.4551432698457065547E-18
 	EXP  607127601336,271623700703		;-.2181804934405659197E-17
 	EXP  607101137417,313245123351		;-.3405122121351518328E-17
 	EXP  170665676575,033607152207		; .2920436655277002656E-17
 	EXP  170554001110,376732276727		; .1192682876882768478E-17
 	EXP  170560060327,321547457416		; .1303606021001427053E-17
 	EXP  170554731631,327217710762		; .1217705177797396539E-17
 
;COEFFICIENTS OF APPROXIMATION POLYNOMIAL ATAN(X) = X*P(X**2)

P06:	EXP 177546176241,173026455171	; .074700604980000000
P05:	EXP 600221360346,262546426550	;-.090879628821850000
P04:	EXP 177570707036,320713526601	; .111110916853003200
P03:	EXP 600133333333,171014105013	;-.142857142198848259
P02:	EXP 177663146314,314620200167	; .199999999998937080
P01:	EXP 600025252525,125252526621	;-.333333333333332690
ONE:	EXP 200140000000,000000000000	; 1.00000000000000000

EPS:	EXP 174372113547	;LARGEST X WITH DATAN(X) = X (HIGH WORD)
				;(IE, ALL DOUBLE PRECISION X WITH HIGH
				; WORD OF X <= EPS HAVE DATAN(X)=X).
MINX:	EXP 177562332734	;TAN(PI/32)	  (HIGH WORD, ROUNDED DOWN)
MAXX:	EXP 200450471543	;1/TAN(PI/32)	  (HIGH WORD, ROUNDED UP)


	SEGMENT	DATA

SGNFLG:	BLOCK	1		;SIGN TO BE ATTACHED TO RESULT
USAVE:	BLOCK	2		;TEMP STORAGE FOR U 

	PRGEND
TITLE	GCOSH	HYPERBOLIC COSINE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 7, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GCOSH
EXTERN	GCOSH.
GCOSH=GCOSH.
PRGEND
TITLE	GCOSH.	HYPERBOLIC COSINE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JUNE 7, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DCOSH(X) IS CALCULATED AS FOLLOWS

;  LET V BE APPROXIMATELY 2 SO THAT LNV AND ABS(X)+LN V
;  CAN BE EXACTLY REPRESENTED WHEN X IS EXACTLY REPRESENTABLE.
;  THEN, LETTING Y = ABS(X), FOR

;       Y <= 709.089565, DCOSH = (EXP(Y)+EXP(-Y))/2,
;       709.089565 <= Y < 1024 * LN(2)
;               DCOSH = (V/2)*EXP(Y - LN V),
;       Y >= 1024 * LN(2), DCOSH = +MACHINE INFINITY

;THE RANGE OF DEFINITION FOR DCOSH IS ABS(X) < 1024 * LN(2) AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE.  DCOSH WILL BE SET
;  TO + MACHINE INFINITY.

;REQUIRED (CALLED) ROUTINES:  GEXP

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GCOSH

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GCOSH,.)	;ENTRY FOR GCOSH ROUTINE
	DMOVE	T0,@(L)		;OBTAIN X
	JUMPGE	T0,DCSH		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;Y = -X
DCSH:	CAMGE	T0,TWENT2	;IF ABS(X) < 22, GO TO
	  JRST 	ALG2		;  FULL CALCULATION.
     	CAMG	T0,BIGX		;IF Y IS .LE. BIGX
	  JRST	EASY		;THEN GO TO EASY
	CAMGE	T0,YMAXHI	;IF HI OF Y < YMAXHI
	  JRST	GETW		;  GO TO GETW
	CAME	T0,YMAXHI	;IF HI OF Y > YMAXHI
	  JRST	MSTK		;GO TO MSTK
	CAMGE	T1,YMAXLO	;HI OF Y = YMAXHI
	  JRST	GETW		; IF LO OF Y .LE. YMAXLO, GO TO GETW
MSTK:	HRLOI	T0,377777	;SET RESULT TO
	SETO	T1,		;MACHINE INFINITY
	$LCALL	ROV
;LERR	(LIB,%,<DCOSH: result overflow>)
	GOODBY	(1)		;RETURN
GETW:	GFAD	T0,LNV		;W = Y-LNV
	DMOVEM	T0,TEMP		;MOVE W TO A TEMP REGISTER
	FUNCT	GEXP.,<TEMP>	;Z = EXP(W)
	DMOVEM	T0,TEMP		;SAVE A COPY OF Z
	GFMP	T0,CON1		;RESULT = Z*CON1
	GFAD	T0,TEMP		;+Z
	  JFCL	MSTK
	GOODBY	(1)		;RETURN
	
EASY:	DMOVEM	T0,TEMP		;MOVE Y TO TEMP
	FUNCT	GEXP.,<TEMP>	;GET ITS EXP
	EXTEND	T0,[GFSC -1]	;DIV BY 2
	GOODBY(1)		;RETURN
ALG2:	CAMG	T0,TM30		;IF ABS(X) .LT. 2**(-30) THE
	  JRST	TINY		;  RESULT IS 1.
     	PUSH	P,T2		;SAVE MORE ACCUMULATORS
	PUSH	P,T3		
	PUSH	P,T4
	PUSH	P,T5
	DMOVEM	T0,TEMP		;MOVE Y TO A TEMPORARY REGISTER
	FUNCT	GEXP.,<TEMP>	;Z = EXP(Y)
	DMOVE	T2,T0		;SAVE A COPY OF Z
	HRLZI	T4,200140	;SET T4 TO 1.0
	SETZ	T5,		;CLEAR SECOND WORD
	GFDV	T4,T2		;1/Z
	GFAD	T0,T4		;+Z
	EXTEND	T0,[GFSC -1]	;*1/2
	POP	P,T5
	POP	P,T4		;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN
TINY:	HRLZI	T0,200140	;RESULT IS 1
	SETZ	T1,
	GOODBY	(1)

ONE:	200140000000				;1.0
BIGX:	201254242673				;709.089565
YMAXHI:	201254271027				;709.782712893383998706
YMAXLO:	367643475715				;
LNV:	DOUBLE	577723506750,010134300000	;-.693147180559947174D0
CON1:	DOUBLE	172041452000,000000000000	;.186417721E-14
TWENT2:	200554000000				;22
TM30:	174340000000				;2**(-30)

	SEGMENT	DATA

TEMP:	DOUBLE	0,0
	PRGEND
TITLE	GEXP2.	DOUBLE ** INTEGER EXPONENTIATION
SUBTTL	MARY PAYNE/MHP/CKS

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1983, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

;	Mary Payne, July 30, 1982

	HELLO	(GEXP2,.)
	DMOVEM	T2,SAVE2	;Save T2-T4
	MOVEM	T4,SAVE4

	DMOVE	T0,ONE		;Floating 1 to T0-T1
	MOVM	T2,@1(L)	;|exponent| to T2
	JUMPE	T2,EXP0		;Exponent = 0 is special
	DMOVE	T3,@0(L)	;Base to T3-T4.
	JUMPN	T3,STEP1	;If base not 0 go to main flow
	JRST	BASE0		;Else to special code

LOOP:	GFMP	T3,T3		;Square current result
	  JOV	OVER2		;Over/underflow possible
STEP1:	TRNE	T2,1		;If exponent is odd
	  GFMP	T0,T3		;  update current result
	    JOV  OVER		;  Branch on over/underflow
	LSH	T2,-1		;Discard low bit of exponent
	JUMPN	T2,LOOP		;Iterate if not 0

	SKIPL	@1(L)		;If exponent > 0
	  JRST	RET		;  return
	DMOVE	T3,T0		;[3243] Copy result
	DMOVE	T0,ONE		;Get reciprocal of
	GFDV	T0,T3		;[3243] result; Underflow impossible
	  JOV	OVMSG		;  On overflow get message
	JRST	RET		;Else return

EXP0:	SKIPE	@0(L)		;Exponent 0. If base not
	  JRST	RET		;  0, result is 1. Return
	$LCALL	ZZZ		;zero**zero is indeterminate, result=zero
	SETZB	T0,T1		;Store 0
	JRST	RET

BASE0:	SKIPL	@1(L)		;If exponent > 0
	  JRST	ZERO		;  result is 0
	$LCALL	ROV		;Else overflow
	HRLOI	T0,377777	;Store +biggest
	HRLOI	T1,377777
	JRST	RET

ZERO:	SETZB	T0,T1		;Result is 0
	JRST	RET		;Return

;
;The following block of code deals with over/underflow in the
;square operation at LOOP:. Note that the "exponent" cannot be
;0 -- LOOP: is entered only if T2 is not 0. Moreover, if T2 is
;not 1 subsequent operations will aggravate the over/underflow
;condition in such a way that both the result of the iteration
;and its reciprocal will have the same exception as currently
;indicated. If, however, T2 = 1, and the square overflowed, it
;is possible that its reciprocal will be in range. We therefore
;complete the current pass through the loop, and if the LSH of T2
;makes it zero, we join the handling at OVER: for overflow/underflow
;on the MUL of T0 by T3. Note that no exception can occur on the
;MUL of T0 by a wrapped over/underflow of the square, so that the
;exception flags will still be valid after this step.
;

OVER2:	GFMP	T0,T3		;No over/underflow. Hence flags
				;  from square of T3 still valid
	LSH	T2,-1		;Discard low bit of exponent
	JUMPE	T2,OVER		;If T2 = 0, T0 has wrapped final
				;  result or its reciprocal
				;  which may be in range

				;Final product surely over/underflows.
	GETFLG	T2		;[3243] get exception flags into t2
	TLNE	T2,(PC%FUF)	;If underflow flag set, reciprocal
	  JRST	UNDER		;  overflows. Go test sign of exponent

	SKIPL	@1(L)		;For overflow, if exponent > 0
	  JRST	UNDMSG		;  final result underflows.
	JRST	OVMSGF		;Else reciprocal gives overflow

;
;The rest of the code handles over/underflow on the product of
;T0 by T3 and calculation of the reciprocal, if this is done.
;

OVER:	GETFLG	T2		;[3243] Get exception flags into t2
	TLNE	T2,(PC%FUF)	;If underflow flag set
	  JRST	UNDER		;  underflow on product
	SKIPL	@1(L)		;Else, overflow on result if
	  JRST	OVMSGF		;  exponent > 0. Get message
	DMOVE	T3,T0		;[3243] Copy result
	DMOVE	T0,ONE		;For exponent < 0, get
	GFDV	T0,T3		;[3243] reciprocal of wrapped overflow
	  JOV	RETF		;Underflow impossible; overflow
				;  compensates previous overflow
	JRST	UNDMSG		;Else, get underflow message

UNDER:	SKIPL	@1(L)		;Product underflowed. If exponent
	  JRST	UNDMSG		;  >/= 0, result underflows
				;Else reciprocal overflows

;[3243] entry point necessary because other branches join OVMSG
OVMSGF:	RESFLG	T2		;[3243] restore the flags, t2-t3 are dbl wd PC
OVMSG:	$LCALL	ROV		;Result overflow
	JUMPL	T0,NEGOV	;If result > 0
	HRLOI	T0,377777	;Store +BIGGEST
	HRLOI	T1,377777
	JRST	RET		;  and return

NEGOV:	MOVSI	T0,400000	;If result < 0, store -BIGGEST
	MOVEI	T1,1
	JRST	RET		;  and return

UNDMSG:	$LCALL	RUN		;Result underflow
	SETZ	T0,
RETF:	RESFLG	T2		;[3243] clear the PC flags
RET:	DMOVE	T2,SAVE2	;Restore T2-T4
	MOVE	T4,SAVE4
	POPJ	P,

ONE:	EXP	200140000000,0	;G-floating 1.0

	SEGMENT	DATA

SAVE2:	BLOCK	2		;Temp for T2-T3
SAVE4:	BLOCK	1		;Temp for T4

	PRGEND
TITLE	GEXP3.	POWER FUNCTION
;		(DOUBLE PRECISION)
SUBTTL	IMSL, INC.	JUNE 4, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.


;GEXP3 CALCULATES X**Y WHERE X AND Y ARE FLOATING POINT VALUES IN THE
;FOLLOWING RANGES
; 0.0 < X < + MACHINE INFINITY  (X MAY EQUAL 0.0 IF Y > 0.0 AND 
; X MAY BE LESS THAN 0.0 IF Y IS AN INTEGER. IF X IS < 0.0
; AND Y IS NOT AN INTEGER, A WARNING ERROR IS ISSUED AND
; DABS(X)**Y IS CALCULATED.)
;  -1024.9375 <= FLOAT(INT((Y*LOG2(X))*16))/16 < 1023.0
;  X**Y IS CALCULATED AS 2**W WHERE W = Y*LOG2(X). LOG2(X) IS 
;  CALCULATED AS FOLLOWS;
;      X = F*(2**M), 1/2 <= F < 1. PICK P SUCH THAT P IS AN ODD
;      INTEGER < 16 AND LET A = 2**(-P/16). NOW X = ((2**M)*A)*(F/A)
;      LOG2(X) = M+LOG2(A) + LOG2(F/A) OR
;      LOG2(X) = M-(P/16) + LOG2(F/A) .
;       LET U1 = M-(P/16) AND
;          U2 = LOG2(F/A) = LOG2((1+S)/(1-S)).
;       THEN LOG2(X) = U1 + U2.
;      AND S = (F-A)/(F+A). A RATIONAL
;      APPROXIMATION IS USED TO EVALUATE U2. U1 AND U2 ARE THEN
;      USED TO DETERMINE W1 AND W2 WHERE W=W1+W2
;      AND W1 = FLOAT(INT(W*16.0))/16.0. FINALLY Z=X**Y=2**W
;      IS RECONSTRUCTED AS Z = (2**W1) * (2**W2) WHERE
;      W1 = M1-P1/16 AND 2**W2 IS EVALUATED FROM ANOTHER RATIONAL
;      FUNCTION.
;
;THE RANGE OF DEFINITION FOR GEXP3 IS GIVEN ABOVE, AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE

;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, T5, P1, P2, P3, AND P4 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,EXP3

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GEXP3,.)		;ENTRY TO GEXP3. ROUTINE
	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3		
	DMOVE	T0,@(L)			;GET THE BASE
	CAMN	T0,[200140,,0]		;IS IT EXACTLY 1.0?
	 JUMPE	T1,POPRET		;YES. JUST RETURN EXACTLY 1.0
	DMOVE	T2,@1(L)		;GET THE EXPONENT
	SETZM	IFLAG			;SET INTEGER EXP FLAG TO 0
	JUMPG	T0,XOK			;IF X IS NOT GT 0
	JUMPE	T0,X0			;AND IF X IS NOT = 0
	PUSH	P,T4			;SAVE ACCUMULATORS
	PUSH	P,T5
	PUSH	P,P1
	PUSH	P,P2
	DMOVE	T4,T2			;GET A COPY OF Y
	JUMPGE	T4,YP			;IF Y IS NEGATIVE
	  DMOVN	T4,T4			;NEGATE IT
YP:	MOVE	P2,T4			;GET HIGH ORDER OF Y
	SETZ	P1,			;SET P1 TO ZERO
	LSHC	P1,14			;GET EXPONENT OF Y
	SUBI	P1,1765			;GET SHIFTING FACTOR
	LSHC	T4,(P1)			;SHIFT OFF EXP AND PART OF INTEGER
	TLNE	T4,400000		;IF Y IS ODD
	  SETOM	IFLAG			;SET IFLAG TO ONES
	LSHC	T4,1			;SHIFT OFF LAST OF INTEGER PART
	DMOVN	T0,T0			;NEGATE X
	JUMPN	T5,MSTK			;IF LO OF Y .NE. 0, GO TO MSTK
	JUMPE	T4,YINT			;IF HI OF Y = 0, GO TO YINT
MSTK:	$LCALL	NNA
;LERR	(LIB,%,<GEXP3: negative ** non-integer; ABS(base) used instead of base>)
	JRST	CONT			;GO TO CONT

X0:	JUMPG	T2,RET1			;IF Y IS NOT GT 0
	JUMPE	T2,ZERZER		;SEPARATE OUT 0**0 CASE
	$LCALL	ZNI
;LERR	(LIB,%,<GEXP3: base is 0 and exponent is le 0; result = infinity>)
	HRLOI	T0,377777		;RESULT = INFINITY
	HRLOI	T1,377777
POPRET:	POP	P,T3			;RESTORE ACCUMULATORS
	POP	P,T2
	GOODBY	(1)			;RETURN

ZERZER:	$LCALL	ZZZ			;"0**0 undefined, result = 0"
	JRST	POPRET			;GO RETURN

XOK:	PUSH	P,T4			;SAVE ACCUMULATORS
	PUSH	P,T5
	PUSH	P,P1
	PUSH	P,P2
YINT:	JUMPN	T3,CONT			;IF LO OF Y IS NOT 0, GO TO CONT
	JUMPN	T2,YONE			;IF Y .NE. 0 GO TO YONE, OTHERWISE
	DMOVE	T0,A1			;SET RESULT TO 1.0
	JRST	RET2			;GO TO RET2
YONE:	CAMN	T2,A1			;IF Y = 1.0 
	  JRST	RET2			;GO TO RET2
CONT:	PUSH	P,P3
	PUSH	P,P4
	MOVE	T4,T0			;OBTAIN THE EXPONENT
	ASH	T4,-30			;SHIFT MANTISSA OFF
	SUBI	T4,2000			;SUBTRACT 1024 FROM EXPONENT
	AND	T0,MASK1		;EXTRACT MANTISSA
	IOR	T0,MASK2		;SET EXPONENT TO 0
					;THE FOLLOWING TESTS DETERMINE P
	MOVEI	T5,2			;NP = 1
	CAMLE	T0,A1+20		;IF F GT A1(9)
	  JRST	NXT2			;THEN GO TO NXT2
	CAME	T0,A1+20		;
	  JRST	OK1			;
	CAMLE	T1,A1+21		;
	  JRST	NXT2
OK1:	ADDI	T5,20			;
NXT2:	CAMLE	T0,A1+6(T5)
	JRST	NXT4			;
	CAME	T0,A1+6(T5)		;
	  JRST	OK2			;
	CAMLE	T1,A1+7(T5)		;
	  JRST	NXT4			;
OK2:	ADDI	T5,10			;
NXT4:	CAMLE	T0,A1+2(T5)		;
	  JRST	NXT3			;
	CAME	T0,A1+2(T5)
	JRST	OK3			;
	CAMLE	T1,A1+3(T5)		;
	  JRST	NXT3			;
OK3:	ADDI	T5,4			;
NXT3:	MOVEM	T5,PTEMP		;SAVE A COPY OF P
	DMOVE	P1,T0			;SAVE A COPY OF F
	GFAD	P1,A1(T5)		;F+A1(2*P+2)
	GFSB	T0,A1(T5)		;F-A1(2*P+2)
	SUBI	T5,2			;FORM [(2*P+2)-
	ASH	T5,-1			;2/2]
	GFSB	T0,A2(T5)		;Z=(F-A1(P+1))-A2((P+1)/2)
	GFDV	T0,P1			;/(F+A1(P+1))
	EXTEND	T0,[GFSC 1]		;
	DMOVE	P1,T0			;SAVE A COPY OF Z
	GFMP	P1,P1			;FORM Z**2
	DMOVE	P3,P1			;SAVE A COPY OF Z**2
	GFMP	P3,RP4			;R(Z)=RP4*Z**2
	GFAD	P3,RP3			;+RP3
	GFMP	P3,P1			;*Z**2
	GFAD	P3,RP2			;+RP2
	GFMP	P3,P1			;*Z**2
	GFAD	P3,RP1			;+RP1
	GFMP	P3,P1			;*Z**2
	GFMP	P3,T0			;*Z
	GFAD	P3,T0			;+Z
	GFMP	P3,C			;U2 = R(Z) * C

	ASH	T4,4			;U1 = M*16
	MOVE	T5,PTEMP		;GET A COPY OF P
	ASH	T5,-1
	SUB	T4,T5			;-P
	FLTR	T4,T4			;FLOAT IT
	SETZ	T5,			;
	EXTEND	T4,[GDBLE T4]		;
	EXTEND	T4,[GFSC -4]		;
	DMOVE	T0,T2			;SAVE A COPY OF Y
	JUMPGE	T2,YPOS			;IF Y IS NEGATIVE
	  DMOVN	T0,T0			;NEGATE IT
YPOS:	MOVE	P1,T0			;GET COPY OF HIGH ORDER OF Y
	ASH	P1,-30			;SHIFT OFF MANTISSA
	CAIGE	P1,2024			;IF THE EXPONENT IS LESS THAN 2024
	  JRST	SMALLY			;GO TO SMALLY
	CAIE	P1,2024			;IF THE EXPONENT IS > 2024
	  JRST	LARGEY			;GO TO LARGEY
	SETZ	T1,			;ZERO SECOND WORD
	  JRST SGNCHK			;GO TO SGNCHK
SMALLY:	SUBI	P1,1775			;GET SHIFTING FACTOR
	SETZ	T1,			;SET Y1 TO 0
	JUMPGE	P1,GETY1			;IF P1 IS .GE. 0
	SETZ	T0,			;THEN GO TO GETY1
	JRST	GETY2			;
GETY1:	AND	T0,MSK3(P1)			;GET Y1
	JRST	SGNCHK			;GO TO CHECK SIGN
LARGEY:	CAIL	P1,2067			;IF EXPONENT IS .GE. 2067
	  JRST	SGNCHK			;GO TO SGNCHK
	SUBI	P1,2025			;GET SHIFTING FACTOR
	AND	T1,REDMSK(P1)			;GET LOW ORDER PART OF Y1
SGNCHK:	JUMPGE	T2,GETY2			;IF Y IS NEGATIVE
	  DMOVN	T0,T0			;NEGATE Y1
GETY2:	DMOVE	P1,T2			;SAVE A COPY OF Y
	GFSB	T2,T0			;Y2 = Y-Y1
	GFMP	T2,T4			;FORM Y2*U1
	JFCL
	GFMP	P1,P3			;FORM Y*U2
	JFCL
	GFAD	T2,P1			;W = Y*U2+Y2*U1
	GFMP	T4,T0			;FORM U1*Y1
	JFCL
	DMOVE	T0,T2			;RECONSTRUCT W
	GFAD	T0,T4
	JFCL
	CAMG	T0,BIGW			;IF W IS NOT TOO BIG
	  JRST	WOK			;GO TO WOK
OVFL:	$LCALL	ROV
;LERR	(LIB,%,<GEXP3: result overflow>)
	HRLOI	T0,377777		; RESULT = INFINITY
	HRLOI	T1,377777
	JRST	RET3			;GO TO RET3
WOK:	CAML	T0,SMALLW		;IF W IS NOT TOO SMALL
	  JRST	WOK2			;THEN PROCEED
UNFL:	$LCALL	RUN
;LERR	(LIB,%,<GEXP3: result underflow>)
	SETZ	T0,			; RESULT = 0
	SETZ	T1,
	POP	P,P4			;RESTORE ACCUMULATORS
	POP	P,P3
	POP	P,P2
	POP	P,P1
	POP	P,T5
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)			;RETURN
WOK2:	SKIPGE	P1,T2			;[4020] GET |W|
	DMOVN	P1,T2			;[4020] THE HARD WAY (IF NEGATIVE)
	SETZ	P2,			;ZERO P2
	MOVE	P3,P1			;
	ASH	P3,-30			;SHIFT OFF MANTISSA
	SUBI	P3,1775			;GET SHIFTING FACTOR
	JUMPGE	P3,GETW1		;IF P3 .GE. 0
					;THEN GO TO GETW1
	SETZ	P1,			;OTHERWISE, SET W1 = 0
	JRST	GETW2			;GO TO GETW2
GETW1:	CAIG	P1,MSKTOP		;[4002] BEYOND TABLE?
	AND	P1,MSK3(P3)		;GETW1
	JUMPG	T2,GETW2		;IF W IS NEGATIVE
	  DMOVN	P1,P1			;NEGATE W1
GETW2:	GFSB	T2,P1			;W2 = W-W1
	GFAD	T4,P1			;W = W1+U1*Y1
	SKIPGE	T0,T4			;[4020] SAVE A COPY OF ABS(W)
	DMOVN	T0,T4			;[4020] THE HARD WAY (IF NEGATIVE)
	MOVE	P1,T0			;
	SETZ	T1,			;ZERO T1
	ASH	P1,-30			;SHIFT OFF MANTISSA
	SUBI	P1,1775			;GET SHIFTING FACTOR
	JUMPGE	P1,GTW1			;IF P1 IS .GE. 0
					;THEN GO TO GTW1
	SETZ	T0,			;OTHERWISE, SET W1 TO 0
	JRST	GTW2			;GO TO GET W2
GTW1:	CAIG	P1,MSKTOP		;[4002] BEYOND TABLE?
	AND	T0,MSK3(P1)		;GET W1
	JUMPGE	T4,GTW2			;IF W IS NEGATIVE
	  DMOVN	T0,T0			;NEGATE W1
GTW2:	GFSB	T4,T0			;FORM W-W1
	GFAD	T2,T4			;W2 = W2+(W-W1)
	SKIPGE	T4,T2			;[4020] SAVE A COPY OF ABS(W2)
	DMOVN	T4,T2			;[4020] THE HARD WAY (IF NEGATIVE)
	MOVE	P1,T4			;
	SETZ	T5,			;ZERO	T5
	ASH	P1,-30			;SHIFT	OFF MANTISSA
	SUBI	P1,1775			;GET SHIFTING FACTOR
	JUMPGE	P1,GW			;IF P1 IS .GE. 0
					;THEN GO TO GW
	SETZ	T4,			;OTHERWISE, SET W=0
	JRST	GW2			;GO TO GW2
GW:	CAIG	P1,MSKTOP		;[4002] BEYOND TABLE?
	AND	T4,MSK3(P1)		;GET W
	JUMPGE	T2,GW2			;IF W2 IS NEGATIVE
	  DMOVN	T4,T4			;NEGATE W
GW2:	GFAD	T0,T4			;FORM W1 + W
	EXTEND T0,[GSNGL T0]		;
	FSC	T0,4			;*16
	FIX	T0,T0			;IW1
	GFSB	T2,T4			;W1 = W2 - W
	JUMPLE	T2,W2POS		;IF W2 .GT. 0
	GFSB	T2,SXTNTH		;W2 = W2-.0625
	ADDI	T0,1			;IW1 = IW1+1
W2POS:	MOVE	T5,T0			;SAVE A COPY OF IW1
	JUMPGE	T5,NPOS
	  ADDI	T5,17
NPOS:	ASH	T5,-4			;M1 = IW1/16
	JUMPL	T0,INEG			;IF IW1 .GE. 0
	ADDI	T5,1			;M1 = M1+1
INEG:	MOVE	T4,T5			;SAVE A COPY OF M1
	ASH	T4,4			;P1 = 16*M1
	SUB	T4,T0			;-IW1
	ASH	T4,1			;
	DMOVE	T0,T2			;SAVE A COPY OF W2
	GFMP	T2,Q7			;Z = Q7*W2
	GFAD	T2,Q6			;+Q6
	GFMP	T2,T0			;*W2
	GFAD	T2,Q5			;+Q5
	GFMP	T2,T0			;*W2
	GFAD	T2,Q4			;+Q4
	GFMP	T2,T0			;*W2
	GFAD	T2,Q3			;+Q3
	GFMP	T2,T0			;*W2
	GFAD	T2,Q2			;+Q2
	GFMP	T2,T0			;*W2
	GFAD	T2,Q1			;+Q1
	GFMP	T0,T2			;*W2
	GFMP	T0,A1(T4)		;*A1(P1+1)
	GFAD	T0,A1(T4)		;+A1(P1+1)
	EXTEND	T0,[GFSC (T5)]	;ADD M1 TO THE EXP OF Z
	JFCL	EXCPT
RET3:	POP	P,P4			;RESTORE ACCUMULATORS
	POP	P,P3
RET2:	SKIPE	IFLAG			;IF Y IS ODD NEGATIVE INTEGER
	  DMOVN	T0,T0			;NEGATE RESULT
	POP	P,P2
	POP	P,P1
	POP	P,T5
	POP	P,T4
RET1:	POP	P,T3
	POP	P,T2
	GOODBY	(1)			;RETURN
EXCPT:	JUMPG	T5,OVFL			;IF T5 >0 GO TO OVFL
	JRST	UNFL			;OTHERWISE, GO TO UNFL

SXTNTH:	DOUBLE	177540000000,000000000000	;.0625D0
BIGW:	201277740000				;UPPER BOUND FOR WW
SMALLW:	576440000000				;LOWER BOUND FOR W
RP1:    DOUBLE  177552525252,252525251443       ;.833333333333332114D-1
RP2:    DOUBLE  177263146314,314740401603       ;.125000000005037992D-01
RP3:    DOUBLE  177044444441,161170665342       ;.223214212859242590D-2
RP4:    DOUBLE  176570743756,332257605031       ;.434457756721631196D-3
C:      DOUBLE  200156125073,051270137606       ;1.442695040888963407D0
Q7:     DOUBLE  176076473357,034454617751       ;.149288526805956082D-4
Q6:     DOUBLE  176450275727,001604616375       ;.154002904409897646D-3
Q5:     DOUBLE  176753541760,306574663707       ;.133335413135857847D-02
Q4:     DOUBLE  177247312533,160461663013       ;.961812905951724170D-02
Q3:     DOUBLE  177470654106,270060115477       ;.555041086640855953D-1
Q2:     DOUBLE  177675376757,374054231615       ;.240226506959095371D0
Q1:     DOUBLE  200054271027,367643475706       ;.693147180559945296D0
A1:     DOUBLE  200140000000,000000000000       ;A1(I), I=1,17 =
A12:    DOUBLE  200075222575,025111033141       ;2**((1-I)/16). THIS
A13:    DOUBLE  200072540306,347672220712       ;TABLE IS SEARCHED
A14:    DOUBLE  200070146336,354125123411       ;TO DETERMINE P.
A15:    DOUBLE  200065642374,312655165530       ;
A16:    DOUBLE  200063422214,025077022007       ;
A17:    DOUBLE  200061263452,021252033327       ;
A18:    DOUBLE  200057204243,237260060666       ;
A19:    DOUBLE  200055202363,063763571444       ;
A110:   DOUBLE  200053254076,352205205126       ;
A111:   DOUBLE  200051377326,251542504707       ;
A112:   DOUBLE  200047572462,140443204215       ;
A113:   DOUBLE  200046033760,121433342514       ;
A114:   DOUBLE  200044341723,163526107032       ;
A115:   DOUBLE  200042712701,343725057267       ;
A116:   DOUBLE  200041325303,147630441731       ;
A117:   DOUBLE  200040000000,000000000000       ;
A2:	DOUBLE	170461734720,307535546011
A22:	DOUBLE	607304062611,120221564632
A23:	DOUBLE	170276106540,343712762662
A24:	DOUBLE	607625040217,333155037217
A25:	DOUBLE	170362230025,031075315002
A26:	DOUBLE	170466404421,360475356413
A27:	DOUBLE	607330176555,216025626545
A28:	DOUBLE	607323056225,270604125171
MASK1:	000077777777				;MASK FOR MANTISSA
MASK2:	200000000000				;MASK FOR EXPONENT
REDMSK:	600000000000
RDMSK2:	700000000000
RDMSK3:	740000000000
RDMSK4:	760000000000
RDMSK5:	770000000000
RDMSK6:	774000000000
RDMSK7:	776000000000
RDMSK8:	777000000000
MASK:	777400000000				;MASKS FOR FINDING W1
MSK1:	777600000000				;
MSK2:	777700000000				;
MSK3:	777740000000				;
MSK4:	777760000000				;
MSK5:	777770000000				;
MSK6:	777774000000				;
MSK7:	777776000000				;
MSK8:	777777000000				;
MSK9:	777777400000				;
MSK10:	777777600000				;
MSK11:	777777700000
MSK12:	777777740000
MSK13:	777777760000
MSK14:	777777770000
MSK15:	777777774000
MSK16:	777777776000
MSK17:	777777777000
MSK18:	777777777400
MSK19:	777777777600
MSK20:	777777777700
MSK21:	777777777740
MSK22:	777777777760
MSK23:	777777777770
MSK24:	777777777774
MSK25:	777777777776
MSK26:	777777777777
MSKTOP==MSK26-MSK3		;MAX INDEX TO USE FOR MSK3

	SEGMENT	DATA

PTEMP:	0
IFLAG:	0
	PRGEND
TITLE	GLOG10	LOG BASE 10 FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      APRIL 8, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GLOG10
EXTERN	GLG10.
GLOG10=GLG10.
PRGEND
TITLE	GLOG	NATURAL LOG FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	APRIL 8, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GLOG
EXTERN	GLOG.
GLOG=GLOG.
PRGEND
TITLE	GLOG.	NATURAL LOG FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	APRIL 8, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.


;DLOG10(X) AND DLOG(X) ARE CALCULATED AS FOLLOWS

;	FOR X = 0 AN ERROR MESSAGE IS ISSUED AND -MACHINE INFINITY
;	IS RETURNED AS THE RESULT. 
;	FOR X < 0 AN ERROR MESSAGE IS ISSUED, X IS SET T0 -X AND 
;	CALCULATION CONTINUES.
;	FOR X > 0, X = F*2**F(M), 1/2 < F < 1
;	DEFINE G AND N SO THAT F = G*2(-N), 1/SQRT(2) <= G < SQRT(2).
;	NOW
;		DLOG(X) = (K*M-N) * DLOG(2) + DLOG(G)
;	AND
;		DLOG10(X) = DLOG10(E) * DLOG(X) = DLOG(X)/DLOG(10)
;
;	DLOG(G) IS EVALUATED BY DEFINING S = (G-1)/(G+1) AND Z = 2*S
;	AND THEN CALCULATING DLOG(G) = DLOG((1+Z/2)/(1-Z/2)) USING
;	A MINIMAX RATIONAL APPROXIMATION.

;THE RANGE OF DEFINITION FOR DLOG/DLOG10 IS GIVEN ABOVE, AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE

;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GLOG  
;	OR
;  PUSHJ	P,GLOG10

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0.
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GLG10.,GLOG10)	;ENTRY TO LOG BASE 10 ROUTINE.
	SETZM	FLAG		;CLEAR FLAG FOR DLOG10 ENTRY
	JRST 	ALG		;GO TO ALGORITHM

	HELLO	(GLOG,.)	;ENTRY TO NATURAL LOG ROUTINE
	SETOM	FLAG		;SET FLAG FOR DLOG ENTRY
ALG:	DMOVE	T0,@(L)		;GET ARG
	JUMPG	T0,STRT		;IF ARG > ZERO GO TO STRT
	JUMPN	T0,ARGN		;OTHERWISE, IF ARG .NE. 0 GO TO ARGN
	$LCALL	AZM
;LERR	(LIB,%,<DLOG or DLOG10: zero arg; result = -infinity>)
	HRLZI	T0,400000	;SET RESULT TO
	HRRZI	T1,000001	;LARGE NEGATIVE NUMBER
	GOODBY	(1)
ARGN:	$LCALL	NAA
;LERR	(LIB,%,<DLOG or DLOG10: negative arg; result = LOG(ABS(arg))>)
	DMOVN	T0,T0		;ARG = -ARG

STRT:	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	HLRZ	T2,T0		;LEFT OF T0 TO RIGHT OF T2
	LSH	T2,-6		;ISOLATE EXPONENT
	SUBI	T2,2000		;SUBTRACT 1024 FROM THE EXP
	DMOVE	T4,T0		;GET COPY OF ARG
	AND	T4,MASK1	;EXTRACT MANTISSA
	IOR	T4,MASK2	;SET EXP TO 0
	CAMLE	T4,HI		;IS HIGH PART OF F > SQRT(.5)
	  JRST	FGT		;YES, GO TO FGT
	CAME	T4,HI		;NO, IS F = HIGH OF SQRT(.5)
	  JRST 	FLT		;NO, GO TO FLT
	CAMLE	T5,LOW		;YES, IS LOW PART > LOW PART OF SQRT(.5)
	  JRST 	FGT		;YES, GO TO FGT
FLT:	SUBI	T2,1		;N = N-1
	EXTEND	T2,[GFLTR T2]
	MOVEM	T2,N		; SAVE N
	GFAD	T4,MHALF	;ZNUM = F-.5
	DMOVE 	T2,T4		;GET COPY OF ZNUM
	EXTEND	T4,[GFSC -1]	;ZDEN = ZNUM * .5
	GFAD	T4,HALF		; + .5
	JRST	EVALRZ
FGT:	EXTEND	T2,[GFLTR T2]
	MOVEM	T2,N		; SAVE N
	DMOVE	T2,T4
	GFAD	T2,B3		;ZNUM = F - 1.0
	EXTEND	T4,[GFSC -1]	;ZDEN = F*.5
	GFAD	T4,HALF		; + .5
EVALRZ:	GFDV	T2,T4		;Z = ZNUM/ZDEN
	DMOVE 	T4,T2		;
	GFMP	T4,T4		;W = Z*Z
	DMOVE	T0,T4		;SAVE COPY OF W
	GFAD	T4,B2		; FORM B(W). B(W) = W + B2
	GFMP	T4,T0		; * W
	GFAD	T4,B1		; + B1
	GFMP	T4,T0		; * W
	GFAD	T4,B0		; + B0
	DMOVEM	T0,SAVEW	; SAVE A COPY OF W
	GFMP	T0,A2		;FORM A(W). A(W)= A2*W
	GFAD	T0,A1		; + A1
	GFMP	T0,SAVEW	; * W
	GFAD	T0,A0		; + A0
	GFDV	T0,T4		;R(Z) = A(W)/B(W)
	GFMP	T0,SAVEW	; * W
	GFMP	T0,T2		; *Z
	GFAD	T0,T2		; + Z
	MOVE	T2,N		;RETRIEVE N
	MOVEI	T3,0		;ZERO OUT SECOND WORD
	GFMP	T2,C1		;FORM N*C1
	MOVE	T4,N		;RETRIEVE N
	MOVEI	T5,0		;ZERO OUT SECOND WORD
	GFMP	T4,C2		;FORM N*C2
	GFAD	T0,T4		;RESULT = N*C2 + R(Z)
	GFAD	T0,T2		; + N*C1
	SKIPN	FLAG
	  GFMP	T0,C3		;IF DLOG10 ROUTINE, RESULT = RESULT*C3
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)

C1:     DOUBLE  200054300000,000000000000       ;.693359375D0                   
C2:     DOUBLE  601310277575,034757152744       ;-2.12194440054690583D-4        
C3:     DOUBLE  177767455730,251156241615       ;.434294481903251828D0          
A0:     DOUBLE  577037740007,152304514557       ;-.641249434237455811D2         
A1:     DOUBLE  200540611121,000552775450       ;.163839435630215342D2          
A2:     DOUBLE  577715357522,145224132710       ;-.789561128874912573D0         
B0:     DOUBLE  576517720013,037446761043       ;-.769499321084948798D3         
B1:     DOUBLE  201147002037,320522317573       ;.312032220919245328D3          
B2:     DOUBLE  577134251775,244603076112       ;-.356679777390346462D2         
B3:     DOUBLE  577640000000,000000000000       ;-.1D1                          
HALF:   DOUBLE  200040000000,000000000000       ;.5D0                           
MHALF:  DOUBLE  577740000000,000000000000       ;-.5D0                          
HI:		200055202363
LOW:		063763571444
MASK1:		000077777777
MASK2:		200000000000

	SEGMENT	DATA

N:		0
FLAG:		0
SAVEW:  DOUBLE  000000000000,000000000000                                       
	PRGEND
TITLE	GMOD	DOUBLE PRECISION REMAINDER FUNCTION
SUBTTL	MARY PAYNE /MHP/CKS	25-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GMOD
EXTERN	GMOD.
GMOD=GMOD.
PRGEND
TITLE	GMOD.	G-FLOATING REMAINDER

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1983, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

; This routine was rewritten during edit 3233.

; G-floating MOD routine with no limit on the magnitudes
; of the arguments. This routine calculates
;
;	|ARG1| - [|ARG1/ARG2|]*|ARG2|
;
; and returns the result with the sign of ARG1. Special code is
; provided for avoiding intermediate exceptions. The final result
; will not overflow. It may underflow, in which case a result of 0
; is returned and an error message is written.
; 
; For A and B positive floating point numbers, let A = 2^n*F and
; B = 2^m*G. We would like to compute R = A - B*int[A/B]. We introduce
; the following sequences of numbers to assist in computing R. For
; p and d positive integers, let J = 2^p*G, I = 2^p*F and define
; I(0) = I - q*J, where q = 0 if I < J and 1 otherwise. For i >= 0 let
;
;                       L(i+1) = 2^d*I(i) 
;
;               I(i+1) = L(i+1) - J*int[L(i+1)/J].
;
; Using an induction argument it is not difficult to show that
;
;            I(k) = 2^(kd)*I(0) - J*int[2^(kd)*I(0)/J].
;
; Now let 0 = < n - m = kd + r, where 0 =< r < d, then 
;
;                 R = A - B*int[A/B]
;                   = 2^n*F - 2^m*G*int[2^(n-m)*F/G]
;     ==>     2^p*R = 2^n*[I(0) + q*J] - 2^m*J*int[2^(n-m)*(I(0) + q*J)/J]
;                   = 2^n*I(0) - 2^m*J*int[2^(n-m)*I(0)]
;     ==> 2^(p-m)*R = 2^(kd+r)*I(0) - J*int[2^(kd+r)*I(0)/J]
;                   = 2^r*I(k) - J*int[2^r*I(k)/J]
;                   = L - J*int[L/J],
;
; where L = 2^r*I(k).
;
; Based on the above, we have the following algorithm for computing R:
;
;             Step 1: m <--- the exponent value of B
;                     c <--- the exponent value of A - m
;                     if c < 0, end with R = A
;             Step 2: I <--- 2^p*(fraction field of A)
;                     J <--- 2^p*(fraction field of B)
;                     If I >= J, I <--- I - J
;                     go to step 5
;             Step 3: L <--- 2^d*I         
;             Step 4: L <--- L - J*int[L/J]
;             Step 5: c <--- c - d
;             Step 6: if c > 0 go to step 3
;	      Step 7: if c = -d, exit with R = L
;             Step 8: L <--- 2^(d+c)*I = 2^r*I
;             Step 9: L <--- L - J*int[L/J]
;             Step 10: R <--- 2^(m-p)*L

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	T6=6			;Additional AC names
	T7=7
	T8=8

DEFINE DMOVM (AC,E) <		;Double absolute value
	DMOVE AC,E
	TLNE AC,(1B0)
	 DMOVN AC,AC
> ;DMOVM

DEFINE DCAML (AC,E) <		;Double integer compare
	CAMN	AC, E
	CAMGE	AC+1, E+1
	CAMLE	AC, E
> ;DCAML

	HELLO (GMOD,.)

	DMOVEM	T2, SAV23	;Save registers 2, 3,
	DMOVEM	T4, SAV45	;  4, 5
	DMOVEM	T6, SAV67	;  6, 7
	MOVEM	T8, SAV8	;  and 8

	DMOVM	T0, @0(L)	; T0 = |A|
	DMOVM	T4, @1(L)	; T4 = |B|
	JUMPE	T4,GMRETZ	;IF ARG2=0, GO RETURN ZERO
;
; Step 1
;
	MOVE	T6, T4		; T6 = |B|hi
	AND	T6, [777700000000]
				; T6 = Exponent field of B (including bias)
	MOVE	T7, T0		; T7 = |A|hi
	SUB	T7, T6		; High 12 bits of T7 = c
	JUMPL	T7, TSTSGN	; Done if c < 0

;
; Step 2
;
STEP2:	ASHC	T6, -30		; Get c in the low bits of T7
				;  and m+2000 in low bits of T6
	TLZ	T0, 777700	; T0 = I
	TLZ	T4, 777700	; T4 = J
	DCAML	T0, T4		; Compare I to J
	DSUB	T0, T4		; If I >= J, I <--- I - J
	JRST	STEP5
;
; Steps 3 through 6
;
STEP3:	SETZB	T2, T3		; T0/T3 = L = 2^35*I -- d = 35
	DDIV	T0, T4		; T0 = int(L/J), T2 = L - J*int(L/J)
	DMOVE	T0, T2		; T0 = L - J*int(L/J)
STEP5:	SUBI	T7, 106		; c <--- c - d
	JUMPG	T7, STEP3	; If c > 0 go to Step 3
;
; Steps 7 and 8 9: At this point c = r - d or r = c + d;
;
	SETZB	T2, T3		; T0/T3 = 2^35*I

				; Shift T0/T3 right T7 places
	CAMG	T7, [-43]	; More than 1 word to shift
	 JRST	[EXCH T1, T2	; Yes, do first 35 bits with word operations
		 EXCH T0, T1
		 MOVN T8, T7	; Get negative shift count
		 ASHC T2, 43(T7)  ; Shift T2-T3
		 ASH  T2, -43(T8) ; Reposition T2
		 ASHC T1, 43(T7)  ; Shift T1-T2
		 JRST STEP8]
	MOVN	T8, T7		; Get negative shift count
	ASHC	T1, (T7)	; Shift T1-T2
	ASH	T1, (T8)	; Reposition T1
	ASHC	T0, (T7)	; Shift T0-T1

STEP8:	DDIV	T0, T4		; T2 = 2^(p-m)*R
;
; Step 9 - Obtain R in floating point format and check for underflow
;
	DMOVE	T0, T2		; Copy fraction into T0
	EXTEND	T0, [GFSC (T6)]	; Insert biased exponent and normalize
	  JFCL	UNDER		; Can underflow

TSTSGN:	SKIPGE	@0(L)		;Remainder in T0. If A < 0
	  DMOVN	T0,T0		;  negate it
RESTOR:	DMOVE	T2,SAV23	;Restore registers 2, 3
	DMOVE	T4,SAV45	;  4, 5
	DMOVE	T6,SAV67	;  6, 7
	MOVE	T8,SAV8		;  and 8
	POPJ	P,		; Return
;
; If processing continues here, the remainder has underflowed.
;
UNDER:	$LCALL	RUN		;Result underflow
	SETZB	T0, T1		;Store 0 for result
	JRST	RESTOR		;Restore registers and return

;
;IF ARG2 IS ZERO, RETURN 0 WITH A WARNING MESSAGE
;
GMRETZ:	SETZB	T0,T1		;RETURN ZERO
	$LCALL	MZZ		;WITH MESSAGE
	JRST	RESTOR

	SEGMENT	DATA
SAV23:	BLOCK	2
SAV45:	BLOCK	2
SAV67:	BLOCK	2
SAV8:	BLOCK	1

	PRGEND
TITLE	GCOS	SINE AND COSINE FUNCTIONS
;	        (DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 1, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GCOS
EXTERN	GCOS.
GCOS=GCOS.
PRGEND
TITLE	GSIN	SINE AND COSINE FUNCTIONS
;	        (DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 1, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GSIN
EXTERN	GSIN.
GSIN=GSIN.
PRGEND
TITLE	GSIN.	SINE AND COSINE FUNCTIONS
;	        (DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 1, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY

;COS(X) AND SIN(X) ARE CALCULATED AS FOLLOWS

;  NOTE THAT COS(X) = COS(-X) AND SIN(X) = -SIN(-X)
;  LET ABS(X) = N*PI + F WHERE ABS(F) <= PI/2
;  WHEN ABS(F) < EPS (SEE CONSTANTS, BELOW), SIN(F) = F
;  OTHERWISE, THE ARGUMENT REDUCTION IS:
;       LET X1 = ABS(X)
;	N = THE NEAREST INTEGER TO ABS(X), FOR SIN, OR
;				TO ABS(X) + PI/2, FOR COS
;       THEN XN = [N/PI], THE GREATEST INTEGER IN N/PI.

;       THEN THE REDUCED ARGUMENT F = (((X1 - XN*C1) -XN*C2) -XN*C3)
;               WHERE C1+C2+C3 = PI TO EXTRA PRECISION AND ARE GIVEN BELOW.

;	LET G = F**2
;		THEN R(G) = (G*XNUM/XDEN+RP1)*G 
;	WHERE XNUM = ((RP5*G+RP4)*G+RP3)*G+RP2
;	      XDEN = ((G*Q2)*G+Q1)*G+Q0
;	AND RP5,RP4,RP3,RP2,RP1,Q2,Q1, AND Q0 ARE GIVEN BELOW
;	AND SIN(F) = F + F*R(G).  THE R(I) APPEAR BELOW

;	FINALLY,
;		SIN(X) = SIGN(X)*SIN(F)*(-1)**N, AND
;		COS(X) = SIN(X + PI/2)

;THE RANGE OF DEFINITION FOR SIN IS ABS(X) <= 1686629713.  ABS(X)+PI/2, IN
;  COS(X), MUST BE LESS THAN 1686629713.  SIN(X) = COS(X) = 0.0 AND AN
;  ERROR MESSAGE WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
 
;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED
 
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
 
;  MOVEI	L,ARG
;  PUSHJ	P,GSIN
;             OR
;  PUSHJ	P,GCOS
 
;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0.
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GCOS,.)	;ENTRY TO COSINE ROUTINE
	DMOVE	T0,@(L)		;OBTAIN ARGUMENT
	PUSH	P,T2		;SAVE AN ACCUMULATOR
	HRLZI	T2,201400	;SET SGN TO 1.
	JUMPGE	T0,XPOS		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;Y = -X
XPOS:	DMOVEM	T0,XDEN		;SAVE A COPY OF ABS(X)
	SETOM	FLAG		;SET FLAG FOR COSINE ENTRY
	GFAD	T0,PID2		;Y = Y+PID2
	JRST	ALG		;PROCEED TO MAIN ALGORITHM

	HELLO	(GSIN,.)	;ENTRY TO SIN ROUTINE
	DMOVE	T0,@(L)		;GET ARG
	PUSH	P,T2		;SAVE AN ACCUMULATOR
	HRLZI	T2,201400	;SET SIGN TO 1.0
	SETZM	FLAG		;CLEAR FLAG FOR SINE ROUTINE
	JUMPGE	T0,ALG		;IF ARG IS NEGATIVE
	DMOVN	T0,T0		;THEN, SET Y=-X AND
	HRLZI	T2,576400	;SET VALUE OF SIGN TO -1.
ALG:	MOVEM	T2,SGN		;MOVE SIGN TO MEMORY
	CAMGE	T0,YMAX1	;IF HI OF Y<HI OF YMAX
	  JRST YOK		;PROCEED TO YOK
	CAME	T0,YMAX1	;IF HI OF Y > HI OF YMAX
	  JRST	ERR1		;GO TO ERR1
	CAMGE	T1,YMAX2	;HI OF Y=HI OF YMAX, IS LO OF Y < LO OF YMAX?
	  JRST	YOK		;YES, PROCEED TO YOK
ERR1:	$LCALL	ATZ
;LERR	(LIB,%,<DSIN or DCOS: ABS(arg) too large; result = zero>)
	MOVEI	T0,0		;SET RESULT
	MOVEI	T1,0		;TO ZERO
	POP	P,T2		;RESTORE ACCUMULATOR
	GOODBY	(1)		;
YOK:	PUSH	P,T3		;SAVE ACCUMULATORS
	PUSH	P,T4		;
	PUSH 	P,T5
	DMOVE	T2,T0		;SAVE A COPY OF ABS(X)
	SKIPE	FLAG		;IF THIS IS THE COS ROUTINE
	  DMOVE	T2,XDEN		;RESTORE ABS(X) TO T2
	GFMP	T0,ODPI		;RN = Y/PI
	JFCL
	EXTEND	T0,[DGFIXR T0]	;FIX RN
	MOVE	T5,T1		;SAVE N
	EXTEND	T0,[DGFLTR T0]	;XN=FLOAT(N)
	TRNE	T5,1		;IS N ODD?
	  MOVNS	SGN		;YES,NEGATE SIGN
	SKIPE	FLAG		;IF THE COSINE IS WANTED
	  GFAD	T0,PT5		;THEN XN=XN-.5
	DMOVE	T4,T0		;SAVE A COPY OF XN
	GFMP	T0,C1		;FORM XN*C1
	GFSB	T2,T0		;F=ABS(X)-(XN*C1)
	DMOVE	T0,T4		;SAVE A COPY OF XN
	GFMP	T0,C3		;FORM XN*C3
	GFMP	T4,C2		;FORM XN*C2
	GFSB	T2,T4		;-XN*C2
	GFSB	T2,T0		;F=F-XN*C3
	DMOVE	T0,T2		;SAVE A COPY OF F
	JUMPGE	T2,FPOS		;IF F IS NEGATIVE
	  DMOVN	T2,T2		;F = -F
FPOS:	CAML	T2,EPS		;IF F IS .GE. EPS
	  JRST	GEEPS		;GO TO GEEPS
	SKIPGE	SGN		;IF SGN IS NEGATIVE
	  DMOVN	T0,T0		;F = -F
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN
GEEPS:	GFMP	T2,T2		;G = F*F
	DMOVE	T4,T2		;SAVE A COPY OF G
	GFAD	T2,Q2		;XDEN = G+Q2
	GFMP	T2,T4		;*G
	GFAD	T2,Q1		;+Q1
	GFMP	T2,T4		;*G
	GFAD	T2,Q0		;+Q0
	DMOVEM	T2,XDEN		;MOVE XDEN TO MEMORY
	DMOVE	T2,T4		;GET A COPY OF G
	GFMP	T2,RP5		;XNUM = G*RP5
	GFAD	T2,RP4		;+RP4
	GFMP	T2,T4		;*G
	GFAD	T2,RP3		;+RP3
	GFMP	T2,T4		;*G
	GFAD	T2,RP2		;+RP2
	GFMP	T2,T4		;*G
	GFDV 	T2,XDEN		;R(G) = XNUM/XDEN
	GFAD	T2,RP1		;+RP1
	GFMP	T2,T4		;*G
	GFMP	T2,T0		;*F
	GFAD	T0,T2		;+F
	SKIPGE	SGN		;IF SGN IS NEGATIVE
	  DMOVN	T0,T0		;THEN NEGATE T0
RET:	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4		;
	POP	P,T3
	POP	P,T2
	GOODBY	(1)

EPS:	174340000000				;2**(-30)
PID2:   DOUBLE 	200162207732,242102643022       ;PI/2
ODPI:   DOUBLE  177750574603,156234420251       ;1/PI
C1:     DOUBLE  200262207732,240000000000       ;HIGH 30 BITS OF PI
C2:     DOUBLE  174442055060,200000000000	;NEXT 28 BITS OF PI
C3:	DOUBLE	171064611431,212134015604	;C1+C2+C3=PI TO 120 BITS
Q2:     DOUBLE  201161273135,127076131616       ;0.394924723520450141 D+3
Q1:     DOUBLE  202142232235,112027153730       ;0.702492288221842518D+5
Q0:     DOUBLE  202751252025,266402115312       ;0.541748285645351853D+7
RP5:    DOUBLE  600516152672,335644427111       ;-0.121560740596710190D-1
RP4:    DOUBLE  200342202301,360464200740       ;0.428183075897778265D+01
RP3:    DOUBLE  576602640645,001056761116       ;-0.489487151969463797D+03
RP2:    DOUBLE  202054054660,302527505074       ;0.451456904704461990D+05
RP1:    DOUBLE  600125252525,125252525242       ;-.166666666666666667D0
PT5:    DOUBLE  577740000000,000000000000       ;-.5
YMAX1:	203762207732				;YMAX = 
YMAX2:	242102643021				;  (PI*2**29)

	SEGMENT	DATA

FLAG:	0
SGN:	DOUBLE	0,0
XDEN:	DOUBLE	0,0
	PRGEND
TITLE	GSINH	HYPERBOLIC SINE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 7, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GSINH
EXTERN	GSINH.
GSINH=GSINH.
PRGEND
TITLE	GSINH.	HYPERBOLIC SINE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JUNE 7, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;GSINH(X) IS CALCULATED AS FOLLOWS

;  LET V BE APPROXIMATELY 2 SO THAT LNV AND ABS(X)+LN V
;  CAN BE EXACTLY REPRESENTED WHEN X IS EXACTLY REPRESENTABLE.
;  THEN, LETTING Y = ABS(X), AND NOTING THAT -SINH(-X)=SINH(X),
;  FOR
;       0 <= Y < EPS, GSINH = Y*SIGN(X)
;       EPS <= Y <= 1, GSINH = X-X*(Z*R(Z)), WHERE Z = Y**2 AND
;       R(Z) IS GIVEN BELOW.
;       1 < Y <= 709.089565, GSINH = SIGN(X)*(EXP(Y)-EXP(-Y))/2,
;       709.089565 <= Y < 1024 * LN(2)
;               GSINH = SIGN(X)*((V/2)*EXP(Y - LN V)),
;       Y >= 1024 * LN(2), SINH = +MACHINE INFINITY * SIGN(X)

;       LET Z = Y**2. THEN
;       R(Z) = (RP0 + Z*(RP1 + Z*(RP2 + Z*RP3)))/(Q0 + Z*(Q1 + Z*(Q2 + Z)))
;       WHERE RP0, RP1, RP2, Q0, Q1, AND Q2 ARE GIVEN BELOW.

;THE RANGE OF DEFINITION FOR GSINH IS ABS(X) < 1024 * LN(2) AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE.  GSINH WILL BE SET
;  TO + MACHINE INFINITY * SIGN(X).

;REQUIRED (CALLED) ROUTINES:  DEXP

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GSINH

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GSINH,.)	;ENTRY TO GSINH ROUTINE
	DMOVE	T0,@(L)		;OBTAIN X
	PUSH	P,T5		;SAVE REGISTER T5
	HRLZI	T5,200140	;SET FLAG TO 1
	JUMPGE	T0,XPOS		;IF X IS NEGATIVE
	  MOVN	T5,T5		;NEGATE FLAG
	  DMOVN	T0,T0		;Y = -X
XPOS:	CAMLE	T0,ONE		;IF Y IS .GT. 1.0
	  JRST	DCSH		;GO TO DCSH
LE:	CAMG	T0,TWOM30	;IF Y IS .LE. 2**-30
	  JRST	RET2		;GO TO RET1
	JUMPGE	T5,NXT		;IF X IS NEGATIVE
	DMOVN	T0,T0		;Y = X
NXT:	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	DMOVE	T2,T0		;SAVE A COPY OF X
	DMOVNM	T0,TEMP		;MOVE A COPY OF -X TO MEMORY
	GFMP	T2,T2		;Z = X*X
	DMOVE	T4,T2		;SAVE A COPY OF Z
	GFAD	T4,Q2		;XDEN = Z + Q2
	GFMP	T4,T2		;*Z
	GFAD	T4,Q1		;+Q1
	GFMP	T4,T2		;*Z
	GFAD	T4,Q0		;+Q0
	DMOVEM	T2,Z		;MOVE A COPY OF Z TO MEMORY
	GFMP	T2,RP3		;XNUM = Z*RP3
	GFAD	T2,RP2		;+RP2
	GFMP	T2,Z		;*Z
	GFAD	T2,RP1		;+RP1
	GFMP	T2,Z		;*Z
	GFAD	T2,RP0		; +RP0
	GFDV	T2,T4		;R(Z) = XNUM/XDEN
	GFMP	T2,Z		;*Z
	GFMP	T2,TEMP		;*(-X)
	GFAD	T0,T2		;+X
	POP	P,T4		;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	POP	P,T5
	GOODBY	(1)		;RETURN
 
DCSH:	CAMGE	T0,TWENT2		;IF ABS(X) < 22, GO TO
	  JRST	ALG2			;  FULL CALCULATION
     	CAMG	T0,BIGX		;IF Y IS .LE. BIGX
	  JRST	EASY		;THEN GO TO EASY
	CAMGE	T0,YMAXHI	;IF HI OF Y < YMAXHI
	  JRST	GETW		;  GO TO GETW
	CAME	T0,YMAXHI	;IF HI OF Y > YMAXHI
	  JRST	MSTK		;GO TO MSTK
	CAMGE	T1,YMAXLO	;HI OF Y = YMAXHI
	  JRST	GETW		;IF LO OF Y .LE. YMAXLO, GO TO GETW
MSTK:	HRLOI	T0,377777	;SET RESULT TO
	SETO	T1,		;MACHINE INFINITY
	  JUMPG	T5,DSNH		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;RESULT = -RESULT
DSNH:	$LCALL	ROV
;LERR	(LIB,%,<DSINH: result overflow>)
	POP	P,T5		;RESTORE ACCUMULATOR
	GOODBY	(1)		;RETURN

GETW:	GFAD	T0,LNV		;W = Y-LNV
	DMOVEM	T0,TEMP		;MOVE W TO A TEMP REGISTER
	FUNCT	GEXP.,<TEMP>	;Z = EXP(W)
	DMOVEM	T0,TEMP		;SAVE A COPY OF Z
	GFMP	T0,CON1		;RESULT = Z*CON1
	GFAD	T0,TEMP		;+Z
	  JFCL	MSTK		;OVERFLOW POSSIBLE
RET2:	JUMPGE	T5,RET3		;IF SGNFLG IS NEGATIVE
	DMOVN	T0,T0		;RESULT = -RESULT
RET3:	POP	P,T5		;RESTORE ACCUMULATOR
	GOODBY	(1)		;RETURN
EASY:	DMOVEM	T0,TEMP		;MOVE Y TO TEMP
	FUNCT	GEXP.,<TEMP>	;GET ITS EXP
	EXTEND	T0,[GFSC -1] 	;DIV BY 2
	JRST	RET2		;GET RIGHT SIGN FOR RESULT.
ALG2:	PUSH	P,T2		;SAVE MORE ACCUMULATORS
	PUSH	P,T3		
	PUSH	P,T4
	DMOVEM	T0,TEMP		;MOVE Y TO A TEMPORARY REGISTER
	FUNCT	GEXP.,<TEMP>	;Z = DEXP(Y)
	DMOVN	T2,T0		;SAVE A COPY OF Z
	MOVEM	T5,SGNFLG	;MOVE FLAG TO MEMORY
	HRLZI	T4,200140	;SET T4 TO 1.0
	SETZ	T5,		;CLEAR SECOND WORD
	GFDV	T4,T2		;1/Z
	GFAD	T0,T4		;+Z
	EXTEND	T0,[GFSC -1]	;*1/2
	POP	P,T4		;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
RET1:	SKIPGE	SGNFLG		;IF SGNFLG IS NEGATIVE
	  DMOVN	T0,T0		;RESULT = -RESULT
RET:	POP	P,T5		;RESTORE ACCUMULATOR
	GOODBY	(1)		;RETURN

RP0:    DOUBLE  202352744232,262463203065       ;.35181283430177117881D+6       
RP1:    DOUBLE  201655127025,264501221757       ;.11563521196851768270D+5       
RP2:    DOUBLE  201050741013,034133711232       ;.16375798202630751372D+3       
RP3:    DOUBLE  200062423475,303374403264       ;.78966127417357099479D+0       
Q0:     DOUBLE  575137624613,372031435521       ;-.21108770058106271242D+7      
Q1:     DOUBLE  202043241271,035545730675       ;.36162723109421836460D+5       
Q2:     DOUBLE  576635220743,361550001577       ;-.27773523119650701667D+3      
LNV:    DOUBLE  577723506750,010134300000       ;-.693147180559947174D0         
ONE:	200140000000				;1.0
BIGX:		201254242673			;709.089565
YMAXHI:		201254271027			;709.782712893383998706
YMAXLO:		367643475715			;
TWOM30:		174340000000			;2**-30
CON1:		172041452000			;.186417721E-14
TWENT2:		200554000000			;

	SEGMENT	DATA

SGNFLG:	0
TEMP:	DOUBLE	0,0
Z:	DOUBLE	0,0
	PRGEND
TITLE	GSQRT	SQUARE ROOT FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	MARCH 19, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GSQRT
EXTERN	GSQRT.
GSQRT=GSQRT.
PRGEND
TITLE	GSQRT.	SQUARE ROOT FUNCTION  
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	MARCH 19, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE  WITH  THE  TERMS OF  SUCH  LICENSE  AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;GSQRT(X) IS CALCULATED AS FOLLOWS

;THE ROUTINE CALCULATES THE SQUARE ROOT OF THE ABSOLUTE VALUE OF A
;DOUBLE PRECISION ARGUMENT BY DOING A LINEAR SINGLE PRECISION 
;APPROXIMATION ON THE HIGH ORDER WORD, FOLLOWED BY TWO SINGLE PRECISION
;NEWTON ITERATIONS AND TWO DOUBLE PRECISION NEWTON ITERATIONS. AN ERROR
;MESSAGE RESULTS FOR NEGATIVE ARGUMENTS.

;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GSQRT

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE

	HELLO	(GSQRT,.)		;ENTRY TO GSQRT ROUTINE

	DMOVE	T0,@(L)			;GET DP ARGUMENT
	JUMPG	T0,GSQRTP		;ARGUMENT IS GREATER THAN 0
	JUMPE	T0,GSQRT4		;ARGUMENT IS ZERO
	$LCALL	NAA
;LERR	(LIB,%,<GSQRT: negative arg; result = GSQRT(DABS(arg))>)

	DMOVN	T0,T0			;ARG = -ARG
GSQRTP:	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	DMOVE	T4,T0			;COPY ARG
	LSH	T4,-1			;COMPUTE LINEAR APPROXIMATION
	TLZE	T4,40
	  JRST	GSQRT2			;YES, GO TO GSQRT2
	ADD	T4,[XWD 26,760700]
	GFMP	T4,[EXP 300145400000,0]
	JRST	GSQRT3
GSQRT2:	ADD	T4,[XWD 26,760700]
	GFMP	T4,[EXP 300165000000,0]
GSQRT3:	DMOVE	T2,T0			;COPY ORIGINAL ARG
	GFDV	T2,T4			;DO NEWTON ITERATIONS
	GFAD	T2,T4
	EXTEND	T2,[GFSC -1]
	DMOVE	T4,T0
	GFDV	T4,T2
	GFAD	T4,T2
	EXTEND	T4,[GFSC -1]
	DMOVE	T2,T0
	GFDV	T2,T4
	GFAD	T2,T4
	EXTEND	T2,[GFSC -1]
	GFDV	T0,T2
	GFAD	T0,T2
	EXTEND	T0,[GFSC -1]
	POP	P,T5			;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
GSQRT4:	GOODBY	(1)			;RETURN
	PRGEND
TITLE   GCOTAN	TANGENT AND COTANGENT FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.      APRIL 16, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GCOTAN
EXTERN	GCOTN.
GCOTAN=GCOTN.
PRGEND
TITLE	GTAN	TANGENT AND COTANGENT FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC       APRIL 16, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GTAN
EXTERN	GTAN.
GTAN=GTAN.
PRGEND
TITLE	GTAN.	TANGENT AND COTANGENT FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC       APRIL 16, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;TAN(X) AND COTAN(X) ARE CALCULATED AS FOLLOWS;

;	THE IDENTITIES
;	 TAN(PI/2.0-G)=1.0/TAN(G)
;	 TAN(N*PI+H)=TAN(H), -PI/2.0 < H <= PI/2.0
;	 TAN(-X)=-TAN(X)
;	 COTAN(X)=1.0/TAN(X)
;	 COTAN(-X)=-COTAN(X)
;	ARE USED TO REDUCE THE ORIGINAL PROBLEM TO A PROBLEM WITH
;	AN ARGUMENT GREATER THAN -PI/2.0 AND LESS THAN OR EQUAL TO
;	PI/2.0, PI=3.14159265358979324.
;	THEN DEFINE N AND F SO THAT
;	 X=N*PI/4.0+F, 0.0 <= F <= PI/4.0, WITH PI=3.14159265358979324.
;	THEN
;	 TAN(F) = F, IF F<2**(-30)
;		= R(F), OTHERWISE
;	 WHERE
;	    R(F)=(((XP3*G+XP2)*G+XP1)*G+XP0)/((((Q4*G+Q3)*G+Q2)*G+Q1)*G+Q0)
;	 AND
;	    G = F*F
;	    XPi AND Qi ARE GIVEN BELOW
;	THE APPROXIMATION IS DERIVED FROM ONE GIVEN IN CODY AND WAITE,
;	"SOFTWARE MANUAL FOR THE ELEMENTARY FUNCTIONS"

;	THE RESULT IS THEN RECIPROCATED, IF NECESSARY, AND GIVEN
;	THE APPROPRIATE SIGN.

;THE RANGE OF DEFINITION FOR GTAN IS 0 < ABS(X) < ((2**29) * (PI/2)) = 843314856.D0

;  AND FOR GCOTAN(X), (2**(-1023))*(1+(2**(-58))) < ABS(X) < ((2**29)*(PI/2)) IS NECESSARY.
;  OTHERWISE, ERROR MESSAGES WILL RESULT.

;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GTAN
;	OR
;  PUSHJ	P,GCOTAN

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0,
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GCOTN.,GCOTAN)	;ENTRY TO GCOTAN ROUTINE
	PUSH	P,T4		;SAVE ACCUMULATORS
	PUSH	P,T5
	SETZM	FLAG		;CLEAR FLAG FOR GCOTAN
	DMOVE	T0,@(L)		;GET ARG
	DMOVE	T4,T0		;Y = X
	JUMPGE	T0,XPOS		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;SET Y = -X
XPOS:	CAMGE	T0,EPS1		;IF THE HI ORDER PART OF Y IS TOO SMALL
	  JRST	ERR1		;THEN GO TO ERR1
	CAME	T0,EPS1		;IF HI OF Y  >  HI OF EPS
	  JRST	YOK		;THEN GO TO YOK
	JUMPN	T1,YOK		;HI OF Y = HI OF EPS, IS LO OF Y  OK?
ERR1:	$LCALL	ROV
;LERR	(LIB,%,DCOTAN: result overflow)
	HRLOI	T0,377777	;SET RESULT TO MACHINE
	HRLOI	T1,377777	;INFINITY
	JUMPGE	T4,RET		;IF ARG IS NEGATIVE
	  DMOVN	T0,T0		;RESULT = - RESULT
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	GOODBY	(1)		;RETURN

	HELLO	(GTAN,.)	;ENTRY TO TAN ROUTINE
	PUSH 	P,T4		;SAVE ACCUMULATORS
	PUSH	P,T5
	SETOM	FLAG		;SET FLAG FOR GTAN ENTRY
	DMOVE	T0,@(L)		;GET COPY OF ARG
	DMOVE	T4,T0		;Y = X
	JUMPGE	T0,YOK		;IF X IS NEGATIVE
	DMOVN	T0,T0		;Y = -X
YOK:	CAMGE	T0,YMAX1	;IF HI OF Y < HI OF YMAX
	  JRST	ALG		;PROCEED TO ALG
	CAME	T0,YMAX1	;IF HI OF Y > HI OF YMAX
	  JRST ERR2		;GO TO ERR2
	CAMG	T1,YMAX2	;HI OF Y=HI OF YMAX, IS LO OF Y .LE. LO OF YMAX?
	  JRST ALG		;YES, PROCEED TO ALG
ERR2:	$LCALL	ATZ
;LERR	(LIB,%,<DTAN or DCOTAN: ABS(arg) too large; result = zero>)
	MOVEI	T0,0		;SET RESULT TO ZERO
	MOVEI	T1,0		;SET SECOND WORD TO ZERO
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	GOODBY	(1)		;RETURN
ALG:	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	CAML	T0,PIO4HI	;COMPARE |ARG| WITH PI/4
	  JRST	REDUCE		; REDUCTION IS NECESSARY
	DMOVE	T2,T0		;|F| TO T2
	DMOVE	T0,T4		;AND SIGNED F TO T0.
	SETZM	N		;STORE 0 FOR N
	JRST	FPOS		;BYPASS REDUCTION

REDUCE:	GFMP	T0,TWODPI	;Y*(2/PI)
	EXTEND	T0,[GFIXR T0]	;FIX IT
	MOVEM	T0,N		;MOVE N TO MEMORY
	EXTEND	T0,[GFLTR T0]	;FLOAT IT
	JUMPLE	T4,NEXT		;IF X IS POSITIVE
	  DMOVN	T0,T0		;NEGATE XN
NEXT:	DMOVE	T2,T0		;GET A COPY OF -XN
	GFMP	T2,C1		;FORM -XN*C1
	GFAD	T4,T2		;F=X - XN*C1
	DMOVE	T2,T0		;GET A COPY OF -XN
	GFMP	T2,C3		;FORM -XN*C3
	GFMP	T0,C2		;FORM -XN*C2
	GFAD	T4,T0		;F=F - XN*C2
	GFAD	T2,T4		;F=F - XN*C3
	DMOVE	T0,T2		;MOVE COPY OF F INTO T0
	JUMPGE	T2,FPOS		;IF F IS NEGATIVE
	  DMOVN	T2,T2		;F = -F
FPOS:	CAML	T2,HIEPS	;IS F < EPS?
	  JRST	NO		;NO, GO TO NO
	HRLZI	T4,200140	;YES , F< EPS,
	MOVEI	T5,0		;XDEN = 1.0
	JRST	FLGCHK		;GO TO CHECK FLAG
	
NO:	GFMP	T2,T2		;NO, F .GE. EPS, SET G=F*F
	DMOVE	T4,T2		;SAVE A COPY OF G
	GFMP	T2,XP3		;XNUM = XP3*G +
	GFAD	T2,XP2		; P2 *
	GFMP	T2,T4		; G +
	GFAD	T2,XP1		; P1 *
	GFMP	T2,T4		; G *
	GFMP	T2,T0		; F +
	GFAD	T0,T2		; F
	DMOVE	T2,T4		;SAVE A COPY OF G
	GFMP	T4,Q4		;XDEN = Q4 * G +
	GFAD	T4,Q3		; Q3 *
	GFMP	T4,T2		; G +
	GFAD	T4,Q2		; Q2 *
	GFMP	T4,T2		; G +
	GFAD	T4,Q1		; Q1 *
	GFMP	T4,T2		; G +
	GFAD	T4,ONE		; 1.0
	
FLGCHK:	MOVE	T2,N		;GET COPY OF N
	SKIPE	FLAG		;IF FLAG IS NOT ZERO
	  JRST	DA		;THEN GO TO DA
	TRNN	T2,1		;OTHERWISE, IS N ODD?
	  JRST	DOVN		;NO, GO GET RESULT
	DMOVN	T0,T0		;XNUM = -XNUM
	JRST	NOVD		;GO TO  NOVD
DA:	TRNN	T2,1		;FLAG MUST BE ONE, IS N ODD?
		
	  JRST	NOVD		;NO, GO TO NOVD
	DMOVN	T0,T0		;XNUM = -XNUM
DOVN:	GFDV	T4,T0		;RESULT = XDEN/XNUM
	DMOVE	T0,T4
	JRST RET		;GO TO RET
NOVD:	GFDV	T0,T4		;RESULT = XNUM/XDEN
RET:	POP	P,T3		;RESTORE ACCUMULATORS
	POP	P,T2
	POP	P,T5
	POP	P,T4
	GOODBY	(1)

XP1:    DOUBLE  600135665120,325457340031       ;-.133383500064219607D0
XP2:    DOUBLE  177070072025,061672533663       ;.342488782358905900D-2
XP3:    DOUBLE  601632425106,212723433423       ;-.178617073422544267D-4
Q1:     DOUBLE  600004205175,300102305265       ;-.466716833397552942D0
Q2:     DOUBLE  177364436365,013742053124       ;.256638322894401129D-1
Q3:     DOUBLE  601227102333,067553367667       ;-.311815319070100273D-3
Q4:     DOUBLE  175441335647,203547435105       ;.498194339937865123D-6
C1:     DOUBLE  200162207732,240000000000       ;HIGH 30 BITS OF PI/2
C2:     DOUBLE  174342055060,200000000000       ;C1+C2+C3=PI/2 TO EXTRA PREC
C3:	DOUBLE	170764611431,212134015604
TWODPI: DOUBLE  200050574603,156234420251
ONE:    DOUBLE  200140000000,000000000000
PIO4HI:		200062207733			;HIGH ORDER PART OF PI/4
EPS1:		000240000000			;HIGH ORDER PART OF EPS1
YMAX1:		203662207732			;HIGH ORDER PART OF YMAX
YMAX2:		242102643021			;LOW ORDER PART OF YMAX
HIEPS:	 	174340000000			; 2**(-31)

	SEGMENT	DATA

N:      0
FLAG:	0
	PRGEND
TITLE	GTANH	HYPERBOLIC TANGENT FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979


;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GTANH
EXTERN	GTANH.
GTANH=GTANH.
PRGEND
TITLE	GTANH.	HYPERBOLIC TANGENT FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

;TANH(X) IS CALCULATED AS FOLLOWS

;  TANH(X) = -TANH(-X).  LET F = ABS(X)
;	IF F > 21.1409890, TANH = 1.0*SIGN(X)
;	IF F > LN(3)/2 AND F <= 21.1409890, TANH = RESULT 1 =
;		SIGN(X)*(1 - 2/(EXP(2*F) + 1)))  
;	IF F < 2**(-29), TANH = F*SIGN(X)

;	OTHERWISE, LET G = F**2, AND OBTAIN RESULT 2 AS
;		TANH = SIGN(X)*(F + F*R(G)), WHERE
;               R(G) = G*(RP0+G*(RP1+RP2*G))/(Q0+G*(Q1+G*(Q2+G))).
;               RP0, RP1, RP2, Q0, Q1, AND Q2 APPEAR BELOW.

;THE RANGE OF DEFINITION FOR TANH IS THE REPRESENTABLE REAL NUMBERS

;REQUIRED (CALLED) ROUTINES:  DEXP

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,TANH

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GTANH,.)		;ENTRY TO GTANH ROUTINE
	DMOVE	T0,@(L)			;OBTAIN X
	PUSH	P,T5			;SAVE AN ACCUMULATOR
	MOVE	T5,T0			;SAVE A COPY OF X
	JUMPGE	T0,XPOS			;IF X IS NEGATIVE
	  DMOVN	T0,T0			;F = -X
XPOS:	CAMG	T0,LN2TC		;IF F IS .LE. LN2TC
	  JRST	CALCF			;GO TO CALCF
	HRLZI	T0,200140		;SET RESULT TO 1.0
	SETZ	T1,			;ZERO SECOND WORD
	JUMPGE	T5,RET1			;IF X IS NEGATIVE
	DMOVN	T0,T0			;RESULT = -RESULT
RET1:	POP	P,T5			;RESTORE ACCUMULATORS
	GOODBY	(1)			;RETURN

CALCF:	PUSH	P,T2			;SAVE MORE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	CAMG	T0,LN3D2		;IF F IS .LE. LN3D2
	  JRST	ALG1			;GO TO ALG1
	EXTEND	T0,[GFSC 1]		;F = F+F
	DMOVEM	T0,TEMP			;SAVE F IN A TEMP REGISTER
	FUNCT	GEXP.,<TEMP>		;EXP(2*F)
	GFAD	T0,ONE			;+1.0
	HRLZI	T2,577540		;SET T2 TO -2.0
	SETZ	T3,			;ZERO SECOND WORD
	GFDV	T2,T0			;-2/(EXP(2*F)+1)
	DMOVE	T0,T2
	GFAD	T0,ONE			;1. - 2/(EXP(2*F)+1)
	JUMPGE	T5,RET			;IF X IS NEGATIVE
	  DMOVN	T0,T0			;RESULT = -RESULT
	POP	P,T4			;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2		
	POP	P,T5
	GOODBY	(1)			;RETURN
ALG1:	CAML	T0,CON1			;IF F IS .GE. 2**(-29)
	JRST	ALG2			;GO TO ALG2
	JUMPGE	T5,RET			;IF X IS NEGATIVE
	  DMOVN	T0,T0			;RESULT = -RESULT
	POP	P,T4			;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	POP	P,T5
	GOODBY	(1)			;RETURN
ALG2:	JUMPGE	T5,POSARG		;IF ARGUMENT < 0
	  DMOVN	T0,T0			;  NEGATE |ARGUMENT|
POSARG:	DMOVEM	T0,TEMP			;SAVE SIGNED ARGUMENT.
	GFMP	T0,T0			;G = F*F
	DMOVE	T2,T0			;SAVE A COPY OF G
	GFAD	T2,Q2			;XDEN = G+Q2
      	GFMP	T2,T0			;*G
	GFAD	T2,Q1			;+Q1
	GFMP	T2,T0			;*G
	GFAD	T2,Q0			; + Q0
	DMOVE	T4,T0			;SAVE A COPY OF G
	GFMP	T4,RP2			;XNUM = G*RP2
	GFAD	T4,RP1			; + RP1
	GFMP	T4,T0			; * G
	GFAD	T4,RP0			; + RP0
	GFMP	T0,T4			; *G
	GFDV	T0,T2			;R(G) = XNUM/XDEN
	GFMP	T0,TEMP			;RESULT = F*R
	GFAD	T0,TEMP			; + F
RET:	POP	P,T4			;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	POP	P,T5
	GOODBY	(1)			;RETURN

RP0:    DOUBLE  576415451321,262036074743               ;-.161341190239962281D+4
RP1:    DOUBLE  577016306122,362132146345               ;-.992259296722360833D+2
RP2:    DOUBLE  577702217271,210105420140               ;-.964374927772254698D0 
Q0:     DOUBLE  201545640742,272351322265               ;.484023570719886887D+4 
Q1:     DOUBLE  201442716132,150037036260               ;.22337720718962312926D+
Q2:     DOUBLE  200770276517,017277377240               ;.112744743805349493D+3 
LN2TC:		200552220277				;21.1409890E0
CON1:		174340000000				;2**(-30)
LN3D2:	      	200043117523             		;.549306144E0
ONE:	DOUBLE	200140000000,000000000000		;1.0

	SEGMENT	DATA

TEMP:	DOUBLE	0,0
	PRGEND
TITLE	GEXP	EXPONENTIAL FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL.INC	APRIL 3, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GEXP
EXTERN	GEXP.
GEXP=GEXP.
PRGEND
TITLE	GEXP.	EXPONENTIAL FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.APRIL 3, 1979

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORDANCE  WITH  THE  TERMS OF  SUCH  LICENSE  AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DEXP(X) IS CALCULATED AS FOLLOWS

;  IF X <= -710.475860073943942, EXP = 0
;  IF X >  709.089565712824051, EXP = +MACHINE INFINITY
;  IF X = 0.0, EXP = 1
;  OTHERWISE,
;       THE ARGUMENT REDUCTION IS:
;		LET X1 = [X], THE GREATEST INTEGER IN X
;			X2 = X - X1
;		    N = THE NEAREST INTEGER TO X/LN(2)
;		THE REDUCED ARGUMENT IS G = ((X1 - N*C1)+X2)+N*C2
;		    WHERE C1 = .543 (OCTAL),
;		    AND C2 IS GIVEN BELOW
;		  THE CALCULATION IS:
;		EXP = R(G)*2**(N+1)
;		    WHERE R(G) = 0.5 + G*P/(Q - G*P)
;			P = ((((P2*G**2)+P1)*G**2)+P0)*G**2
;			Q = (((((Q3*G**2)+Q2)*G**2)+Q1)*G**2)+Q0
;		P0, P1, P2, Q0, Q1, Q2, AND Q3 ARE GIVEN BELOW AS
;		XP0, XP1, XP2, XQ0, XQ1, XQ2, AND XQ3 .

;THE RANGE OF DEFINITION FOR DEXP IS GIVEN ABOVE, AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE

;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, AND T5 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  MOVEI	L,ARG
;  PUSHJ	P,GEXP

;THE HIGH ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T0,
;THE LOW ORDER PART OF THE ANSWER IS RETURNED IN ACCUMULATOR T1.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GEXP,.)		;ENTRY TO GEXP ROUTINE
	DMOVE	T0,@(L)			;GET DP ARGUMENT
	CAMLE	T0,SMLXHI		;IF HI OF X > SMLXHI
	  JRST	NXTCHK			;  GO TO NXTCHK
	CAME	T0,SMLXHI		;IF HI OF X < SMLXHI
	  JRST OUT2			;  GO TO OUT2
	CAMGE	T1,SMLXLO		;HI PART = SMLXHI,
	  JRST 	OUT2			;  IF LO OF X < SMLXLO, OUT2
NXTCHK:	CAMGE	T0,BIGXHI		;IF HI OF X < BIGXHI
	  JRST	EXP1			;  GO TO EXP1
	CAME	T0,BIGXHI		;IF HI OF X > BIGXHI
	  JRST	MSTK			;  GO TO MSTK
	CAMG	T1,BIGXLO		;IF LO OF X .LE. BIGXLO,
	  JRST	EXP1			;  GO TO EXP1

MSTK:	$LCALL	ROV
;LERR	(LIB,%,<DEXP: result overflow>)
	HRLOI	T0,377777		;DEXP = +MACHINE INFINITY
	HRLOI	T1,377777
	GOODBY	(1)			;RETURN

OUT2:	$LCALL	RUN
;LERR	(LIB,%,<DEXP: result underflow>)
	MOVEI 	T0,0			;EXP = 0
	MOVEI	T1,0
	GOODBY	(1)			;RETURN

EXP1:	JUMPE	T0,[MOVSI T0,200140	;RETURN 1.0 FOR ARG OF ZERO
		GOODBY (1)]		;EXIT
	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3	
	PUSH	P,T4
	PUSH	P,T5
	MOVM	T2,T0			;GET ABS(ARGHI)
	CAML	T2,LN2OV2		;IF .GE. (LN(2))/2
	  JRST	REDUCE			;  MUST REDUCE ARGUMENT.
	SETZ	T4,
	MOVEM	T4,SAVEN		;SET N TO ZERO.
	DMOVE	T4,T0			;GET COPY OF ARG.
	JRST	MERGE			;MERGE WITH MAIN FLOW.

REDUCE:	DMOVE	T2,T0			;GET COPY OF ARG
	GFMP	T2,RNDLN2		;ARG/LN2
	EXTEND	T2,[GFIXR T2]		;NEAREST INTEGER = N
	EXTEND	T2,[GFLTR T2]		;FLOAT N
	MOVEM	T2,SAVEN		;SAVEN
	GFMP	T2,C1			;N*C1
	GFAD	T0,T2			;X + N*C1
	MOVE	T2,SAVEN		;RETRIEVE N
	MOVEI	T3,0			;ZERO SECOND WORD
	GFMP	T2,C2			;FORM N*C2
	GFAD	T0,T2			;N*C2 + (N*C1 + X)
	DMOVE	T4,T0			;SAVE G
	MOVM	T2,T4			;GET ABS(G)
MERGE:	CAML	T2,TWOM30		;IF REDUCED ARG IS>= 2**-30
	  JRST	APPRX			;  GO TO APPRX
	GFAD	T0,ONE			;R(G) = 1. + G
	EXTEND	T0,[GFSC -1]		;*.5
	JRST	BRNCH			;GO TO BRNCH
APPRX:	GFMP	T4,T4			;Z = G*G
	DMOVE	T2,T4			;SAVE Z
	GFMP	T4,XP2			;Z*XP2
	GFAD	T4,XP1			;+XP1
	GFMP	T4,T2			;* Z
	GFAD	T4,XP0			;+ XP0
	GFMP	T0,T4			;* G
	DMOVE	T4,T2			;SAVE Z
	GFMP	T4,XQ3			;XQ3*Z
	GFAD	T4,XQ2			;+XQ2
	GFMP	T4,T2			;*Z
	GFAD	T4,XQ1			;+ XQ1
	GFMP	T4,T2			;* Z
	GFAD	T4,XQ0			; + XQ0
	GFSB	T4,T0			; XQ - G*XP
	GFDV	T0,T4			;(G*XP)/(XQ-G*XP)
	GFAD	T0,XQ0			; + .5
BRNCH:	MOVE	T4,SAVEN		;RETRIEVE N
	SETZ	T5,			;T5=0
	EXTEND	T4,[GFIX T4]
	ADDI	T4,1			;N = N+1
	EXTEND	T0,[GFSC 0(T4)]		;ADD N TO THE EXPONENT

	POP	P,T5			;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2

RET:	GOODBY	(1)			; RETURN

SMLXHI:	576523460613				;-710.475860073943942
SMLXLO:	202140360224				;
BIGXHI:	201254242673				;709.089565712824051
BIGXLO:	161647554056				;
RNDLN2: DOUBLE  200156125073,051270137606       ;1.44269504088896341 = 	1/LN2   
ONE:    DOUBLE  200140000000,000000000000       ;1.0D0                          
TWOM30: DOUBLE  174340000000,000000000000       ;2**-30                         
C1:     DOUBLE  577723500000,000000000000       ;-0.693359375D0                 
C2:     DOUBLE  176467500202,343020625033       ;2.12194440054690583D-4         
XP0:    DOUBLE  177740000000,000000000000       ;0.250D0                        
XP1:    DOUBLE  177176035137,221241124545       ;0.757531801594227767D-2        
XP2:    DOUBLE  176241055041,127151103610       ;0.315551927656846464D-4        
XQ0:    DOUBLE  200040000000,000000000000       ;0.5D0                          
XQ1:    DOUBLE  177472134502,216775477655       ;0.568173026985512218D-1        
XQ2:    DOUBLE  176651274142,341215233732       ;0.631218943743985036D-3        
XQ3:    DOUBLE  175462315430,147120212512       ;0.751040283998700461D-6        
LN2OV2:	DOUBLE	177754271027,367643475715

	SEGMENT	DATA

SAVEN:	0
	PRGEND
TITLE	GINT	DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GINT
EXTERN	GINT.
GINT=GINT.
PRGEND
TITLE	GABS	DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GABS
EXTERN	DABS.
GABS=DABS.
PRGEND
TITLE	GABS.	DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GABS.
EXTERN	DABS.
GABS.=<DABS.+0>			;[4015]
PRGEND
TITLE	GMAX1	DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GMAX1
EXTERN	DMAX1.
GMAX1=DMAX1.
PRGEND
TITLE	GMAX1.	DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GMAX1.
EXTERN	DMAX1.
GMAX1.=<DMAX1.+0>			;[4015]
PRGEND
TITLE	GMIN1	DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GMIN1
EXTERN	DMIN1.
GMIN1=DMIN1.
PRGEND
TITLE	GMIN1.	DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GMIN1.
EXTERN	DMIN1.
GMIN1.=<DMIN1.+0>			;[4015]
PRGEND
TITLE	GSIGN	DOUBLE PRECISION TRANSFER OF SIGN FUNCTION
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GSIGN
EXTERN	DSIGN.
GSIGN=DSIGN.
PRGEND
TITLE	GSIGN.	DOUBLE PRECISION TRANSFER OF SIGN FUNCTION
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1980, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GSIGN.
EXTERN	DSIGN.
GSIGN.=<DSIGN.+0>			;[4015]
PRGEND
TITLE	DTOGA	
SUBTTL	M. R. BOUCHER/MRB		11-JUN-84

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	DTOGA
EXTERN	DTOGA.
DTOGA=DTOGA.
PRGEND
TITLE	DTOG	DOUBLE PRECISION TO GFLOAT DOUBLE PRECISION CONVERSION

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

					;ARG BLOCK OFFSETS
	SRC==0				;SOURCE ARRAY ADDRESS
	DST==1				;DESTINATION ARRAY ADDRESS
	CNT==2				;NUMBER OF ITEMS TO CONVERT

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(DTOG)
					;CONVERT ONE NUMBER
	DMOVE	T0,@(L)			;GET THE D NUMBER
	PUSHJ	P,DTOGS			;CONVERT IT
	GOODBY

	HELLO	(DTOGA.)		;[4012]
					;CONVERT AN ARRAY OF NUMBERS
	PUSH	P,T3			;SAVE SOME AC'S
	PUSH	P,T4
	PUSH	P,T5
	XMOVEI	T3,@SRC(L)		;GET D PNTR
	XMOVEI	T4,@DST(L)		;GET G PNTR
	MOVE	T5,@CNT(L)		;GET COUNT
DTOGLP:	DMOVE	T0,(T3)			;GET A D NUMBER
	PUSHJ	P,DTOGS			;CONVERT IT
	DMOVEM	T0,(T4)			;SAVE IT
	ADDI	T3,2			;INCR PNTRS
	ADDI	T4,2
	SOJG	T5,DTOGLP		;LOOP FOR ARRAY
	POP	P,T5
	POP	P,T4
	POP	P,T3
	GOODBY

DTOGS:	JUMPL	T0,NEGDTG		;DO SEPARATELY IF NEGATIVE
	JUMPE	T0,DTGZER		;AND NOTHING IF ZERO
	PUSHJ	P,EXPDTG		;NOW DO THE MOST STUFF
	POPJ	P,

NEGDTG:	DMOVN	T0,T0			;MAKE IT POSITIVE
	PUSHJ	P,EXPDTG		;DO MOST STUFF
	DMOVN	T0,T0			;MAKE IT NEGATIVE AGAIN
	POPJ	P,

DTGZER:	DMOVE	T0,[EXP 0,0]		;LOAD ALL ZEROES
	POPJ	P,

EXPDTG:	PUSH	P,T2			;SAVE AN AC
	CAMN	T0,[377777,,-1]		;'OVERFLOW' AMOUNT?
	 CAME	T1,[377777,,-1]
	  JRST	DTGOK			;NO
	JRST	DTGDON			;YES. LEAVE AS IS
DTGOK:	LDB	T2,[POINT 8,T0,8]	;GET THE EXPONENT
	ADDI	T2,1600			;CONVERT TO G-TYPE EXP
	TLZ	T0,777000		;CLEAR THE EXPONENT
	DADD	T0,[EXP 0,4]		;ROUND UP
	ASHC	T0,-3			;SHIFT FRACTION
	TLNN	T0,100			;DID WE OVERFLOW?
	 JRST	EXPOK			;NO
	ASHC	T0,-1			;YES. SHIFT 1 MORE
	ADDI	T2,1			;AND INCR THE EXPONENT
EXPOK:	DPB	T2,[POINT 12,T0,11]	;DROP IN THE EXP
DTGDON:	POP	P,T2			;RESTORE THE AC
	POPJ	P,

	PRGEND
TITLE	GTODA
SUBTTL	M. R. BOUCHER/MRB		11-JUN-84

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1973, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

SEARCH	MTHPRM
NOSYM
ENTRY	GTODA
EXTERN	GTODA.
GTODA=GTODA.
PRGEND
TITLE	GTOD	GFLOAT DOUBLE PRECISION TO DOUBLE PRECISION CONVERSION

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

					;ARG BLOCK OFFSETS
	SRC==0				;SOURCE ARRAY ADDRESS
	DST==1				;DESTINATION ARRAY ADDRESS
	CNT==2				;NUMBER OF ITEMS TO CONVERT

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GTOD)
					;CONVERT ONE NUMBER
	DMOVE	T0,@(L)			;GET THE G NUMBER
	PUSHJ	P,GTODS			;CONVERT IT
	GOODBY

	HELLO	(GTODA.)		;[4012]
					;CONVERT AN ARRAY OF NUMBERS
	PUSH	P,T3			;SAVE SOME AC'S
	PUSH	P,T4
	PUSH	P,T5
	XMOVEI	T3,@SRC(L)		;GET D PNTR
	XMOVEI	T4,@DST(L)		;GET G PNTR
	MOVE	T5,@CNT(L)		;GET COUNT
GTODLP:	DMOVE	T0,(T3)			;GET A G NUMBER
	PUSHJ	P,GTODS			;CONVERT IT
	DMOVEM	T0,(T4)			;SAVE IT
	ADDI	T3,2			;INCR PNTRS
	ADDI	T4,2
	SOJG	T5,GTODLP		;LOOP FOR ARRAY
	POP	P,T5
	POP	P,T4
	POP	P,T3
	GOODBY

GTODS:	JUMPL	T0,NEGGTD		;DO SEPARATELY IF NEGATIVE
	JUMPE	T0,GTDZER		;AND NOTHING IF ZERO
	PUSHJ	P,EXPGTD		;NOW DO THE MOST STUFF
	POPJ	P,

NEGGTD:	DMOVN	T0,T0			;MAKE IT POSITIVE
	PUSHJ	P,EXPGTD		;DO MOST STUFF
	DMOVN	T0,T0			;MAKE IT NEGATIVE AGAIN
	POPJ	P,

GTDZER:	DMOVE	T0,[EXP 0,0]		;LOAD ALL ZEROES
	POPJ	P,

EXPGTD:	PUSH	P,T2			;SAVE AN AC
	CAMN	T0,[377777,,-1]		;'OVERFLOW' AMOUNT?
	 CAME	T1,[377777,,-1]
	  JRST	GTDOK			;NO
	JRST	GTDDON			;YES. LEAVE AS IS
GTDOK:	LDB	T2,[POINT 11,T0,11]	;GET THE EXPONENT
	SUBI	T2,1600			;CONVERT TO D-TYPE EXP
	JUMPL	T2,GTDLOW		;EXPONENT TO SMALL
	CAIL	T2,400			;TEST IF TOO BIG
	 JRST	GTDHGH			;YUP. IT'S TO BIG
	TLZ	T0,777700		;CLEAR THE EXPONENT
	ASHC	T0,3			;SHIFT THE FRACTION
	DPB	T2,[POINT 9,T0,8]	;DROP IN THE EXP
	JRST	GTDDON
GTDLOW:	$LCALL	RUN
;LERR	(LIB,%,<GTOD: result underflow>)
	DMOVE	T0,[EXP 0,0]		;LOAD ZEROES
	JRST	GTDDON			;AND LEAVE
GTDHGH:	$LCALL	ROV
;LERR	(LIB,%,<GTOD: result overflow>)
	DMOVE	T0,[EXP 377777777777,377777777777] ;LOAD AN OVERFLOW
GTDDON:	POP	P,T2			;RESTORE THE AC
	POPJ	P,

	PRGEND

TITLE	GSN.0

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.0
GSN.0:	GSNGL	0
	PRGEND


TITLE	GSN.2
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.2
GSN.2:	GSNGL	2
	PRGEND


TITLE	GSN.4
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.4
GSN.4:	GSNGL	4
	PRGEND


TITLE	GSN.6
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.6
GSN.6:	GSNGL	6
	PRGEND


TITLE	GSN.10
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.10
GSN.10:	GSNGL	10
	PRGEND


TITLE	GSN.12
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.12
GSN.12:	GSNGL	12
	PRGEND


TITLE	GSN.14
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GSN.14
GSN.14:	GSNGL	14
	PRGEND

TITLE	GDB.0

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.0
GDB.0:	GDBLE	0
	PRGEND


TITLE	GDB.2
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.2
GDB.2:	GDBLE	2
	PRGEND


TITLE	GDB.4
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.4
GDB.4:	GDBLE	4
	PRGEND


TITLE	GDB.6
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.6
GDB.6:	GDBLE	6
	PRGEND


TITLE	GDB.10
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.10
GDB.10:	GDBLE	10
	PRGEND


TITLE	GDB.12
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.12
GDB.12:	GDBLE	12
	PRGEND


TITLE	GDB.14
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDB.14
GDB.14:	GDBLE	14
	PRGEND

TITLE	GFX.0

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.0
GFX.0:	GFIX	0
	PRGEND


TITLE	GFX.2
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.2
GFX.2:	GFIX	2
	PRGEND


TITLE	GFX.4
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.4
GFX.4:	GFIX	4
	PRGEND


TITLE	GFX.6
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.6
GFX.6:	GFIX	6
	PRGEND


TITLE	GFX.10
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.10
GFX.10:	GFIX	10
	PRGEND


TITLE	GFX.12
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.12
GFX.12:	GFIX	12
	PRGEND


TITLE	GFX.14
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFX.14
GFX.14:	GFIX	14
	PRGEND

	TITLE	GFL.0

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1979, 1987
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.0
GFL.0:	GFLTR	0
	PRGEND


TITLE	GFL.2
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.2
GFL.2:	GFLTR	2
	PRGEND


TITLE	GFL.4
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.4
GFL.4:	GFLTR	4
	PRGEND


TITLE	GFL.6
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.6
GFL.6:	GFLTR	6
	PRGEND


TITLE	GFL.10
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.10
GFL.10:	GFLTR	10
	PRGEND


TITLE	GFL.12
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.12
GFL.12:	GFLTR	12
	PRGEND


TITLE	GFL.14
	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GFL.14
GFL.14:	GFLTR	14

	PRGEND
	TITLE GDIM G-floating positive difference

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	ENTRY	GDIM
	EXTERN	GDIM.
	GDIM=GDIM.
	PRGEND

	TITLE	GDIM. G-Floating Positive Difference.
	SUBTTL	C. McCutcheon	6/26/81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

; Fortran Intrinsic Function for Positive Difference.
; Passed and returned arguments are G-floating double precision.

; Passed: A1,A2
; Returns:	(from page 15-23 of standard)
; 	A1-A2	if A1 .GT. A2
; 	0	if A1 .LE. A2

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL
	HELLO	(GDIM,.)

	PUSH	P,T2
	PUSH	P,T3

	DMOVE	T0,@0(L)		;A1
	DMOVE	T2,@1(L)		;A2

	CAMN	T0,T2			;If A1 .LT. A2,
	CAML	T1,T3
	CAMGE	T0,T2
	JRST	RET0			; return 0.

	GFSB	T0,@1(L)		;Subtract A1-A2
	JFCL	EXCEP			;Can underflow and overflow

RET:	POP	P,T3
	POP	P,T2
	GOODBYE				;Return

RET0:	SETZB	T0,T1			;Return zero
	JRST	RET

EXCEP:	JUMPE	T0,UNDER		;Underflow?
	$LCALL	ROV,RET			;No, result overflow
UNDER:	$LCALL	RUN,RET			;Result underflow

	PRGEND
	TITLE	GINT.  G-Floating truncation.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

; Fortran Instrinsic Function to return a G-floating double
; precision truncation of the double precision number passed.

; If Magnitude(A) .LT. 1, then 0 is returned,
; otherwise return the largest integer that does not exceed the
; magnitude of A, and whose sign is the same as that of A.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GINT,.)	;Enter routine

	DMOVE	T0,@0(L)	;Get argument to truncate.
	JUMPN	T0,NZERO	;If zero passed, return.
	GOODBYE			;Return

NZERO:	PUSH	P,T2		;Save ac.

	CAIGE	T0,0		;If original number .LT. zero,
	 DMOVN	T0,T0		; negate it.

	HLRZ	T2,T0		;Get exponent
	LSH	T2,-6		;Put rightmost
	CAIG	T2,^O2000	;[3215] If exponent .LE. 2000 then
	 JRST	ZERO		; return zero (number .LT. 1)
	CAIL	T2,^O2073	;If exponent .GE. 2073 then 
	 JRST	DONE		; return the number passed.

; Now shift out the fractional part of the number.

	ASHC	T0,-2073(T2)	;Shift into integer position.
	MOVN	T2,T2		;Negate exponent.
	ASHC	T0,2073(T2)	;Shift back to where found.

	SKIPGE	@0(L)		;If original number .LT. zero,
	 DMOVN	T0,T0		; negate it again.

BYE:	POP	P,T2		;Pop saved ac's
	GOODBYE			;Return

DONE:	DMOVE	T0,@0(L)	;No fractional part, return
	JRST	BYE		; the number passed.

ZERO:	SETZB	T0,T1		;Return 0
	JRST	BYE

	PRGEND
	TITLE	GNINT.  G-Floating nearest whole number.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

; Fortran Instrinsic Function to return a G-floating double
; precision nearest whole number of the double precision number
; passed.

; Number returned is:	Integer(A+.5)	if A .GE. 0
; 			Integer(A-.5)	if A .LT. 0.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GNINT,.)	;Enter routine

	DMOVE	T0,@0(L)		;Get argument to truncate.
	JUMPN	T0,NZERO	;If zero passed,
	GOODBYE			; return.

NZERO:	PUSH	P,T2		;Save ac's
	CAIGE	T0,0		;If original number .LT. zero,
	 DMOVN	T0,T0		; negate it.

; Now Shift out the fractional part of the number.

	HLRZ	T2,T0		;Get exponent
	LSH	T2,-6		;Put rightmost
	CAIL	T2,^O2073	;If exponent .GE. 2073 then 
	 JRST	DONE		; return the number passed.
	CAIGE	T2,2000		;number much .LT. 1.
	 JRST	ZERO		; return 0.

	GFAD	T0,[200040,,0
			0,,0]	;Add .5 before truncation.

	HLRZ	T2,T0		;Get exponent
	LSH	T2,-6		;Put rightmost
	ASHC	T0,-2073(T2)	;Shift into integer position.
	MOVN	T2,T2		;Negate exponent.
	ASHC	T0,2073(T2)	;Shift back to where found.

	SKIPGE	@0(L)		;If original number .LT. zero,
	 DMOVN	T0,T0		; negate it again.

BYE:	POP	P,T2		;Pop saved ac's
	GOODBYE			;Return

DONE:	DMOVE	T0,@0(L)		;No fractional part, return
	JRST	BYE		; the number passed.

ZERO:	SETZB	T0,T1		;Set to 0.
	JRST	BYE

	PRGEND
	TITLE	GPROD.  G-Floating product for single prec. factors.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

; Fortran Instrinsic Function to return a double precision
; product for two real arguments passed it.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(GPROD,.)	;Enter routine
	EXTEND	T0,[GDBLE @1(L)] ;Convert 2nd arg to G-floating.
	DMOVEM	T0,MULT		;Save 2nd argument

	EXTEND	T0,[GDBLE @0(L)] ;Convert 1st arg to G-floating.

	GFMP	T0,MULT		;Multiply and leave result in AC0.
	JFCL	EXCEP		;Result can under/overflow

RET:	POPJ	P,		;Return

EXCEP:	JUMPE	T0,UNDER	;Underflow?
	$LCALL	ROV,RET
UNDER:	$LCALL	RUN,RET

	SEGMENT	DATA

MULT:	BLOCK	2		;Multiplier
	PRGEND
	TITLE	IGNIN.  Integer nearest whole number for G-Gloating.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (c) DIGITAL EQUIPMENT CORPORATION 1981, 1987
;ALL RIGHTS RESERVED.

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE

;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

;THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT  NOTICE
;AND  SHOULD  NOT  BE  CONSTRUED  AS A COMMITMENT BY DIGITAL EQUIPMENT
;CORPORATION.

;DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY  OF  ITS
;SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL.

; Fortran Instrinsic Function to return the nearest
; integer to the G-floating double precision number passed.

; Number returned is:	Integer(A+.5)	if A .GE. 0
; 			Integer(A-.5)	if A .LT. 0.

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(IGNIN.,IGNINT)	;Enter routine

	DMOVE	T0,@0(L)	;Get argument to round.
	JUMPN	T0,NZERO	;If number passed = 0,
	GOODBYE			;Return

NZERO:	CAIGE	T0,0		;If original number .LT. zero,
	 DMOVN	T0,T0		; negate it.
	TLNN	T0,200000	;Is |number| .LT. 1/2 ?
	 JRST	ZERO		;Yes, go return zero
	PUSH	P,T2		;Save an ac

	GFAD	T0,[200040,,0
			 0,,0]	;Add 0.5G0 before truncation.

; Now shift out the fractional part of the number.

	HLRZ	T2,T0		;Get exponent
	LSH	T2,-6		;Put rightmost
	CAILE	T2,2043		;If exponent .GE. 2043 then 
	 JRST	DONE		; return largest integer

	TLZ	T0,777700	;Eliminate exponent.
	ASHC	T0,-2030(T2)	;Shift into integer position,
				; integer result in T0.

	SKIPGE	@0(L)		;If original number .LT. zero,
	 MOVN	T0,T0		; negate it again.

BYE:	POP	P,T2		;Pop saved ac's
	GOODBYE			;Return

; Number is .LT. 1/2, return zero quickly (don't use the normal flow
; because GFAD 0.5G0 will round)

ZERO:	SETZ	T0,		;Return . . .
	GOODBYE			; . . . 0

; Number too large to represent as integer.  Return largest
; integer.

DONE:	$LCALL	ROV
;LERR	(LIB,%,<IDNINT: Result overflow>)
	HRLOI	T0,377777	;Largest positive integer.
	SKIPGE	@0(L)		;If original number .LT. zero,
	 HRLZI	T0,400000	; largest negative integer.
	JRST	BYE		;

	END