Google
 

Trailing-Edge - PDP-10 Archives - BB-D480F-SB_FORTRAN10_V10 - 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(4015)

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;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.
\
	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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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:	SETZ	P2,			;ZERO P2
	MOVM	P1,T2			;GET HIGH ORDER OF W
	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
	MOVM	T0,T4			;SAVE A COPY OF ABS(W)
	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)
	MOVM	T4,T2			;SAVE A COPY OF ABS(W2)
	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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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, 1985
;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