Google
 

Trailing-Edge - PDP-10 Archives - BB-4157E-BM - fortran-ots-debugger/fordbl.mac
There are 3 other files named fordbl.mac in the archive. Click here to see a list.
	SEARCH	FORPRM
	TV	FORDBL	DOUBLE PRECISION ROUTINES ,6(2031)

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 *****

1101	SWG	16-Aug-79
	Cleanup FORDAR for V6 - take out F40 and KA conditionals
	remove following macros: FLADD, DFA, DAM, DFS, DSM
	FLMUL, DFM, DMM, FLDIV, DFD, DDM, all of which are KA
	only library routines.  Replace DFN in IDF macro
	with DMOVN.  remove KI conditionals because no
	longer necessary

1351	EDS	16-Mar-81	Q10-04786
	Fix TWOSEG and RELOC problems.

1405	DAW	6-Apr-81
	Extended addressing support.

1464	DAW	12-May-81
	New error message format.

1673	CDM	9-Sept-81
	Added new double precision routines.
	(DDIM, DINT, DNINT, DPROD, IDNINT)

1713	BL	15-Sep-81
	Added 'DFL.3'.

1721	CDM	15-Sep-81
	Edited DDIM to eliminate DDIM(-INF,+INF) overflow message.

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

1746	CDM	25-Sep-81
	Simple change to DDIM.

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

\
	PRGEND
TITLE	DACOS	ARC SINE AND ARC COSINE FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 19, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DACOS
EXTERN	DACOS.
DACOS=DACOS.
PRGEND
TITLE	DASIN	ARC SINE AND ARC COSINE FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 19, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DASIN
EXTERN	DASIN.
DASIN=DASIN.
PRGEND
TITLE	DASIN.	ARC SINE AND ARC COSINE FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JUNE 19, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 + 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:

;  XMOVEI	L,ARG
;  PUSHJ	P,DASIN
;	OR
;  PUSHJ	P,DACOS

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DACOS,.)	;ENTRY TO DACOS ROUTINE
	HRRZI	T0,2		;SET FLAG TO TWO
	MOVEM	T0,FLAG		;MOVE FLAG TO MEMORY
	JRST	GETX		;GO TO GETX
	
	HELLO	(DASIN,.)	;ENTRY TO DASIN 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:	CAMGE	T0,ONE		;IF Y IS .LT. ONE
	  JRST	ALG		;GO TO ALG
	CAMN	T0,ONE		;IF Y IS .GT. ONE, GO TO ERR
	  JUMPE	T1,ALG
ERR0:	LERR	(LIB,%,DASIN or DACOS; ABS(arg) > one; 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,201400	;SET T2,T3 TO
	SETZ	T3,		;ONE
	DFSB	T2,T0		;G = 1-Y
	FSC	T2,-1		;* 1/2
	DMOVEM	T2,TEMP		;SAVE A COPY OF G
	FUNCT	DSQRT.,<TEMP> ;Y = DSQRT(G)
	FSC	T0,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,TNM11H	;IF Y < TNM11H
	JRST	GOTRES		;GO TO GOTRES
     	DMOVE	T2,T0		;SAVE A COPY OF Y
	DFMP	T2,T2		;G = Y*Y
GOTY:	DMOVE	T4,T2		;SAVE A COPY OF G
	DFAD	T4,Q4		;XDEN = G + Q4
	DFMP	T4,T2		;*G
	DFAD	T4,Q3		;+Q3
	DFMP	T4,T2		;*G
	DFAD	T4,Q2		;+Q2
	DFMP	T4,T2		;*G
	DFAD	T4,Q1		;+Q1
	DFMP	T4,T2		;*G
	DFAD	T4,Q0		;+Q0
	DMOVEM	T2,TEMP		;SAVE A COPY OF G
	DFMP	T2,RP5		;XNUM = G*RP5
	DFAD	T2,RP4		;+RP4
	DFMP	T2,TEMP		;*G
	DFAD	T2,RP3		;+RP3
	DFMP	T2,TEMP		;*G
	DFAD	T2,RP2		;+RP2
	DFMP	T2,TEMP		;*G
	DFAD	T2,RP1		;+RP1
	DFMP	T2,TEMP		;*G
	DFDV	T2,T4		;RESULT = XNUM/XDEN
	DFMP	T2,T0		;*Y
	DFAD	T0,T2		;+Y
GOTRES:	MOVE	T5,I		;GET A COPY OF I
	SKIPE	FLAG		;IF FLAG IS .NE. 0
	  JRST	NE		;GO TO NE
	DFAD	T0,A(T5)	;RESULT = RESULT + A(I)
	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
	DFAD	T0,A(T5)	;+A(I)
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN

NEGX:	DFAD	T0,B(T5)	;RESULT = RESULT + B(I)
	POP	P,T5		;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN

RP1:	DOUBLE	572112065225,361372635617	;-.27368494524164255994D+2
RP2:	DOUBLE	206711524715,202552212153	;.57208227877891731407D+2
RP3:	DOUBLE	571302372325,226222540767	;-.39688862997504877339D+2
RP4:	DOUBLE	204504702731,073034447126	;.10152522233806463645D+2
RP5:	DOUBLE	577233210222,206360633365	;-.69674573447350646411D+0
Q0:	DOUBLE	567267447760,165074063632	;-.16421096714498560795D+3
Q1:	DOUBLE	211641111704,007534102703	;.41714430248260412556D+3
Q2:	DOUBLE	566202106100,352253670303	;-.38186303361750149284D+3
Q3:	DOUBLE	210455717445,226216157457	;.15095270841030604719D+3
Q4:	DOUBLE	572202642744,101474740606	;-.23823859153670238830D+2
TNM11H:	134537657770				;HIGH ORDER PART OF 1.D-11
A:	DOUBLE	0,0
A1:	DOUBLE	201622077325,021026430215	;PI/2
B:	DOUBLE	202622077325,021026430215	;PI
B1:	DOUBLE	201622077325,021026430215	;PI/2
ONE:		201400000000			;1.0
HALF:		200400000000			;1/2

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DATAN2
EXTERN	DATN2.
DATAN2=DATN2.
PRGEND
TITLE	DATAN	ARC TAN FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 9, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DATAN
EXTERN	DATAN.
DATAN=DATAN.
PRGEND
TITLE	DATAN.	ARC TAN FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	SALL

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



	HELLO (DATAN,.)		;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,9		;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-2000+24(T2) ;GET OFFSET INTO XHI TABLES

	DMOVE	T3,T0		;GET A COPY OF X
	DFSB	T0,XHI(T2)	;GET X-XHI
	DFMP	T3,XHI(T2)	;    X*XHI
	DFAD	T3,ONE		;    1 + X*XHI
	DFDV	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
	DFMP	T3,T3		;GET Z**2
	DMOVE	T5,P06		;GET P(Z**2)
	DFMP	T5,T3
	DFAD	T5,P05
	DFMP	T5,T3
	DFAD	T5,P04
	DFMP	T5,T3
	DFAD	T5,P03
	DFMP	T5,T3
	DFAD	T5,P02
	DFMP	T5,T3
	DFAD	T5,P01
	DFMP	T3,T5
	DFMP	T3,T0		; * Z
	DFAD	T0,T3		; + Z = ATAN(Z)

SMALLX:	DFAD	T0,ATANLO(T2)	;  + ATAN(XHI) LOW
	DFAD	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
	DFDV	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 (DATN2.,DATAN2)	;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
	DFDV	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,9		;GET INDEX INTO OFFSET TABLES
	ASHC	T2,3
	HLRZ	T2,OFFSET-2000+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
	DFSB	T7,T5		;GET LOW 35 BITS OF V = VLO
	DFMP	T5,XHI(T2)	;GET V*XHI = VHI * XHI
	DFMP	T7,XHI(T2)	;	    + VLO * XHI
	DFSB	T0,T5		;GET (U - VHI*XHI)
	DFSB	T0,T7		;		   - VLO*XHI
	DMOVE	T5,USAVE	;GET U BACK
	DFMP	T5,XHI(T2)	;GET U * XHI
	DFAD	T3,T5		;GET V + U*XHI
	DFDV	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
	DFDV	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
	  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
	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    175657740000,0		; .1054534912109375000    X1
 	EXP    176407740000,0		; .1288757324218750000    X2
 	EXP    176477740000,0		; .1562194824218750000    
 	EXP    176617640000,0		; .1952209472656250000    
 	EXP    176777400000,0		; .2497558593750000000    
 	EXP    177477540000,0		; .3121948242187500000    
 	EXP    177617200000,0		; .3898925781250000000    
 	EXP    177776300000,0		; .4984130859375000000    
 	EXP    200515740000,0		; .6522216796875000000    
 	EXP    200674040000,0		; .8673095703125000000    
 	EXP    201453440000,0		; 1.170166015625000000    
 	EXP    201645100000,0		; 1.645019531250000000    
 	EXP    202511040000,0		; 2.570800781250000000    
 	EXP    203527000000,0		; 5.359375000000000000    X14

ATANHI:	EXP  000000000000,000000000000		; .0000000000000000000    ATAN(0)
 	EXP  175656261520,251717717115		; .1050651824695016827    ATAN(X1)
 	EXP  176406373154,305103547431		; .1281692623534198424    
 	EXP  176475276500,226075230141		; .1549669515088296698    
 	EXP  176612661305,353341503017		; .1927961221331547523    
 	EXP  176765175625,123431412160		; .2447488705187413410    
 	EXP  177465675077,207755453362		; .3026068193134274483    
 	EXP  177574536624,346575570716		; .3717628303965550497    
 	EXP  177731362665,337061624571		; .4623772720674932799    
 	EXP  200447716236,341667321523		; .5779354489017924548    
 	EXP  200555632634,072362412013		; .7144577222199199908    
 	EXP  200672140433,153152445005		; .8636495725719767213    
 	EXP  201406227204,252435004161		; 1.024591515455200379    
 	EXP  201463117110,112456461235		; 1.199822549393342831    
 	EXP  201542714675,142761043345		; 1.386328658261967261    ATAN(X14)
PI2: 	EXP  201622077325,021026430215		; 1.570796326794896619    PI/2
MPI: 	EXP  575155700452,356751347563		;-3.141592653589793238    -PI
 	EXP  575173246105,164207746145		;-3.036527471120291556    -PI + ATAN(X1)
 	EXP  575176220221,273215536144		;-3.013423391236373396    -PI + ATAN(X2)
 	EXP  575201554376,370255221171		;-2.986625702080963569    
 	EXP  575206433527,115527433724		;-2.948796531456638486    
 	EXP  575215150344,104133030272		;-2.896843783071051897    
 	EXP  575224470162,337747115121		;-2.838985834276365790    
 	EXP  575235354335,213631126655		;-2.769829823193238188    
 	EXP  575251036741,252657532242		;-2.679215381522299958    
 	EXP  575267664122,247327234107		;-2.563657204688000784    
 	EXP  575311247221,375446052165		;-2.427134931369873248    
 	EXP  575334330561,311604060764		;-2.277943081017816517    
 	EXP  575361014155,104167751653		;-2.117001138134592860    
 	EXP  576016720236,050401400603		;-1.941770104196450407    
 	EXP  576076516023,100703762713		;-1.755263995327825977    -PI + ATAN(X14)
 	EXP  576155700452,356751347563		;-1.570796326794896619    -PI/2

ATANLO:	EXP  000000000000,000000000000		; .0000000000000000000    
 	EXP  701136265552,021570267746		;-.1105497383317563340E-19
 	EXP  075632563733,177450755317		; .5435918950106886920E-20
 	EXP  077603757522,264676361123		; .2053885961335437547E-19
 	EXP  077723405617,212207745036		; .2474984160195000159E-19
 	EXP  700044202740,142243000710		;-.2518569148307390217E-19
 	EXP  100702445617,347031410060		; .4770635578227287401E-19
 	EXP  100414524241,267545144713		; .2844597940229199763E-19
 	EXP  100625025353,262726507370		; .4288548085099372680E-19
 	EXP  676255576133,135473517361		;-.7162797697818810113E-19
 	EXP  676323642342,201773750626		;-.6356616556133364840E-19
 	EXP  077413353720,052424757433		; .1415925447541254247E-19
 	EXP  100744743210,367241225634		; .5134543068800780929E-19
 	EXP  101423305250,056513174673		; .5831512827058072309E-19
 	EXP  102617464410,310763545272		; .1692382723892392993E-18
 	EXP  101611431424,134015603464		; .8333742918520878328E-19
 	EXP  675166346353,243762174314		;-.1666748583704175666E-18
 	EXP  102634261642,105051607712		; .1746358738541957411E-18
 	EXP  103601511025,077765605721		; .3266520381981663157E-18
 	EXP  676115710653,365124065054		;-.9192589013278796940E-19
 	EXP  675060707135,225203171017		;-.1961351253927427867E-18
 	EXP  675072766647,260206474405		;-.1918605498534914687E-18
 	EXP  101716137637,073361174657		; .9787193190895619423E-19
 	EXP  674134635612,010745612637		;-.3550693134652264558E-18
 	EXP  700335536464,205477162116		;-.1536916027087339636E-19
 	EXP  103746522614,251310022042		; .4122184681426969927E-18
 	EXP  103760133656,162370070313		; .4202802795595514454E-18
 	EXP  101457607713,122451564536		; .6432483060209586270E-19
 	EXP  103567657506,360715220731		; .3183514413117920163E-18
 	EXP  676000222177,166457565523		;-.1083597300998368435E-18
 	EXP  074603276433,074574160521		; .2563414018821732668E-20
 	EXP  676166346353,243762174314		;-.8333742918520878328E-19
 
;COEFFICIENTS OF APPROXIMATION POLYNOMIAL ATAN(X) = X*P(X**2)

P06:	EXP 175461762413,330264551705	; .0747006049800000000
P05:	EXP 602213603465,225464265500	;-.0908796288218500000
P04:	EXP 175707070366,207135266005	; .1111109168530032000
P03:	EXP 601333333333,310141050127	;-.1428571421988482590
P02:	EXP 176631463146,146202001672	; .1999999999989370798
P01:	EXP 600252525252,252525266207	;-.3333333333333326895
ONE:	EXP 201400000000,000000000000   ; 1.000000000000000000

EPS:	EXP 142561156421	;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 175623327343	;TAN(PI/32)	  (HIGH WORD, ROUNDED DOWN)
MAXX:	EXP 204504715427	;1/TAN(PI/32)	  (HIGH WORD, ROUNDED UP)


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

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DCOSH
EXTERN	DCOSH.
DCOSH=DCOSH.
PRGEND
TITLE	DCOSH.	HYPERBOLIC COSINE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JUNE 7, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 <= 88.029678, DCOSH = (EXP(Y)+EXP(-Y))/2,
;       88.029678 <= Y < 128 * LN(2)
;               DCOSH = (V/2)*EXP(Y - LN V),
;       Y >= 128 * LN(2), DCOSH = +MACHINE INFINITY

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

;REQUIRED (CALLED) ROUTINES:  DEXP

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,DCOSH

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DCOSH,.)	;ENTRY FOR DCOSH 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	ERRD		;GO TO ERRD
	CAMGE	T1,YMAXLO	;HI OF Y = YMAXHI
	  JRST	GETW		; IF LO OF Y .LE. YMAXLO, GO TO GETW
ERRD:	HRLOI	T0,377777	;SET RESULT TO
	SETO	T1,		;MACHINE INFINITY
	LERR	(LIB,%,DCOSH: result overflow)
	GOODBY	(1)		;RETURN

GETW:	DFAD	T0,LNV		;W = Y-LNV
	DMOVEM	T0,TEMP		;MOVE W TO A TEMP REGISTER
	FUNCT	DEXP.,<TEMP>	;Z = EXP(W)
	DMOVEM	T0,TEMP		;SAVE A COPY OF Z
	FMP	T0,CON1		;RESULT = Z*CON1
	SETZ	T1,		;ZERO SECOND WORD
	DFAD	T0,TEMP		;+Z
	  JFCL	ERRD
	GOODBY	(1)		;RETURN
EASY:	DMOVEM	T0,TEMP		;MOVE Y TO TEMP
	FUNCT	DEXP.,<TEMP>	;GET ITS EXPONENT
	FSC	T0,-1		;DIV BY 2
	GOODBY	(1)		;RETURN
ALG2:	CAMG	T0,TM32		;IF ABS(X) < 2**(-32)
	JRST	TINY		; THE RESULT = 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	DEXP.,<TEMP>	;Z = DEXP(Y)
	DMOVE	T2,T0		;SAVE A COPY OF Z
	HRLZI	T4,201400	;SET T4 TO 1.0
	SETZ	T5,		;CLEAR SECOND WORD
	DFDV	T4,T2		;1/Z
	DFAD	T0,T4		;+Z
	FSC	T0,-1		;*1/2
	POP	P,T5
	POP	P,T4		;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	GOODBY	(1)		;RETURN
TINY:	HRLZI	T0,201400	;RESULT IS 1.
	SETZ	T1,
	GOODBY	(1)		;RETURN

ONE:	201400000000				;1.0
BIGX:		207540074635			;88.029691
YMAXHI:	207542710277				;88.722839111672999604
YMAXLO:	276434757153				;
LNV:	DOUBLE 	577235067500,101343000000	;-.693147180559947174D0
CON1:		120414520000			;.186417721E-14
TWENT2:		205540000000	;22
TM32:		141400000000	;2**(-32)

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DEXP
EXTERN	DEXP.
DEXP=DEXP.
PRGEND
TITLE	DEXP.	EXPONENTIAL FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	APRIL 3, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 <= -89.415986292232944914, EXP = 0
;  IF X >  88.029691931113054295, 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:

;  XMOVEI	L,ARG
;  PUSHJ	P,DEXP

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DEXP,.)		;ENTRY TO DEXP 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	ERRD			;  GO TO ERR
	CAMG	T1,BIGXLO		;IF LO OF X .LE. BIGXLO,
	  JRST	EXP1			;  GO TO EXP1

ERRD:	LERR	(LIB,%,DEXP: result overflow)
	HRLOI	T0,377777		;DEXP = +MACHINE INFINITY
	HRLOI	T1,377777
	GOODBY	(1)			;RETURN

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

EXP1:	JUMPE	T0,[MOVSI T0,(1.0)	;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 |ARGUMENT| IN T2
	CAML	T2,LN2OV2		;IS IT .LT. (1/2)*LN(2)?
	  JRST	REDUCE			;IF SO, IT MUST BE REDUCED.
	SETZ	T4,			;MAKE SCALE INDEX 0.
	MOVEM	T4,SAVEN		;  AND SAVE IT.
	DMOVE	T4,T0			;GET COPY OF ARGUMENT.
	JRST	MERGE			;MERGE WITH MAIN FLOW.

REDUCE:	DMOVE	T2,T0			;GET COPY OF ARG
	DFMP	T2,RNDLN2		;ARG/LN2
	FIXR	T2,T2			;NEAREST INTEGER = N
	FLTR	T2,T2			;FLOAT N
	MOVEI	T3,0			;
	FIX	T4,T0			;[X]
	FLTR	T4,T4			;X1 = FLOAT[X]
	MOVEI	T5,0
	DFSB	T0,T4			;X2 = X - X1
	MOVEM	T2,SAVEN		;SAVEN
	DFMP	T2,C1			;N*C1
	DFAD	T4,T2			;X1 + (N*C1)
	DFAD	T4,T0			;+X2
	MOVE	T2,SAVEN		;RETRIEVE N
	MOVEI	T3,0			;ZERO SECOND WORD
	DFMP	T2,C2			;FORM N*C2
	DFAD	T4,T2			;(N*C2)+X1+(N*C1)+X2
	DMOVE	T0,T4			;SAVE G
	MOVM	T2,T4			;GET ABS(G)

MERGE:	CAML	T2,TWOM32		;IF REDUCED ARG IS>= 2**-32
	  JRST	APPRX			;  GO TO APPRX
	DFAD	T0,ONE			;R(G) = 1. + G
	FSC	T0,-1			;*.5
	JRST	BRNCH			;GO TO BRNCH
APPRX:	DFMP	T4,T4			;Z = G*G
	DMOVE	T2,T4			;SAVE Z
	DFMP	T4,XP2			;Z*XP2
	DFAD	T4,XP1			;+XP1
	DFMP	T4,T2			;* Z
	DFAD	T4,XP0			;+ XP0
	DFMP	T0,T4			;* G
	DMOVE	T4,T2			;SAVE Z
	DFMP	T4,XQ3			;XQ3*Z
	DFAD	T4,XQ2			;+XQ2
	DFMP	T4,T2			;*Z
	DFAD	T4,XQ1			;+ XQ1
	DFMP	T4,T2			;* Z
	DFAD	T4,XQ0			; + XQ0
	DFSB	T4,T0			; XQ - G*XP
	DFDV	T0,T4			;(G*XP)/(XQ-G*XP)
	DFAD	T0,XQ0			; + .5
BRNCH:	MOVE	T4,SAVEN		;RETRIEVE N
	FIX	T4,T4
	ADDI	T4,1			;N = N+1
	FSC	T0,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:	570232254036				;-89.415986292232944914
SMLXLO:	301750634730				;
BIGXHI:	207540074636				;88.029691931113054295
BIGXLO:	077042573256				;
LN2OV2:	177542710300				;(1/2) * LN(2)
RNDLN2:	DOUBLE 201561250731,112701376057	;1.44269504088896341 = 	1/LN2
ONE:	DOUBLE	201400000000,000000000000	;1.0D0
TWOM32:	DOUBLE	141400000000,000000000000	;2**-32
C1:	DOUBLE 577235000000,000000000000	;-0.693359375D0
C2:	DOUBLE 164675002027,030206250331	;2.12194440054690583D-4
XP0:	DOUBLE 177400000000,000000000000	;0.250D0
XP1:	DOUBLE 171760351374,212411245446	;0.757531801594227767D-2
XP2:	DOUBLE 162410550412,271511036101	;0.315551927656846464D-4
XQ0:	DOUBLE 200400000000,000000000000	;0.5D0
XQ1:	DOUBLE 174721345024,167754776545	;0.568173026985512218D-1
XQ2:	DOUBLE 166512741427,012152337316	;0.631218943743985036D-3
XQ3:	DOUBLE 154623154303,071202125117	;0.751040283998700461D-6

	RELOC					;DATA
SAVEN:	0
	RELOC
	PRGEND
TITLE	DEXP2.	DOUBLE ** INTEGER EXPONENENTIATION
SUBTTL	CHRIS SMITH/CKS		28-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DEXP2,.)
	PUSH	P,T2		;SAVE TEMP REGISTERS
	PUSH	P,T3
	PUSH	P,T4

	DMOVE	T0,[EXP 201400000000,0]	;SET RESULT TO 1
	MOVM	T4,@1(L)	;GET ABS(EXPONENT)
	JUMPE	T4,D2RET	;EXPONENT ZERO, GO RETURN 1

	DMOVE	T2,@0(L)	;GET BASE
	JUMPN	T2,D2POS	;BASE NOT ZERO, GO GET BASE**ABS(EXP)
	JRST	D2ZERO		;BASE ZERO, GO HANDLE

D2LP:	DFMP	T2,T2		;SQUARE BASE
	  JFCL	D2OVUN		;CATCH OVERFLOW
D2POS:	TRNE	T4,1		;CHECK LOW BIT OF EXPONENT
	  DFMP	T0,T2		;MULTIPLY ANSWER BY BASE
	   JFCL	D2OVUN		;CATCH OVERFLOW
	LSH	T4,-1		;DIVIDE EXPONENT BY 2
	JUMPN	T4,D2LP		;HANDLE ALL BITS OF EXPONENT

	SKIPL	@1(L)		;NEGATIVE EXPONENT?
	  JRST	D2RET		;NO, WE HAVE RESULT
	DMOVE	T2,[EXP 201400000000,0]	;YES, GET RECIPROCAL
	DFDV	T2,T0
	DMOVE	T0,T2

D2RET:	POP	P,T4		;RESTORE TEMP REGISTERS
	POP	P,T3
	POP	P,T2
	POPJ	P,		;DONE

D2ZERO:	SKIPL	T4,@1(L)	;BASE 0, RESULT IS 0 IF EXPONENT POSITIVE
	  JRST	D2RTZ		;POSITIVE EXPONENT, GO RETURN 0
	JRST	D2OVFL		;ELSE RESULT OVERFLOWS

D2OVUN:	SKIPG	T4,@1(L)	;NEGATIVE EXPONENT?
	  JRST	D2UNFL		;YES, RESULT UNDERFLOWS
D2OVFL:	LERR	(LIB,%,DEXP2: result overflow)
	HRLOI	T0,377777	;OVERFLOW, GUESS POSITIVE RESULT
	HRLOI	T1,377777
	SKIPL	@0(L)		;CHECK SIGN OF BASE
	  JRST	D2RET		;BASE NONNEGATIVE, GO RETURN +INFINITY
	TRNE	T4,1		;ODD EXPONENT?
	  DMOVN	T0,T0		;NEGATIVE**ODD, RETURN -INFINITY
	JRST	D2RET

D2UNFL:	LERR	(LIB,%,DEXP2: result underflow)
D2RTZ:	SETZB	T0,T1		;RETURN 0
	JRST	D2RET

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.


;DEXP3 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.)
;  -128.9375 < FLOAT(INT((Y*LOG2(X))*16))/16 < 127.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 DEXP3 IS GIVEN ABOVE, AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE

;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2, T3, T4, T5, G1, G2, G3, AND G4 WERE SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DEXP3,.)		;ENTRY TO DEXP3. ROUTINE
	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3		
	DMOVE	T0,@(L)			;GET THE BASE
	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,G1
	PUSH	P,G2
	DMOVE	T4,T2			;GET A COPY OF Y
	JUMPGE	T4,YP			;IF Y IS NEGATIVE
	  DMOVN	T4,T4			;NEGATE IT
YP:	MOVE	G2,T4			;GET HIGH ORDER OF Y
	SETZ	G1,			;SET G1 TO ZERO
	LSHC	G1,11			;GET EXPONENT OF Y
	SUBI	G1,170			;GET SHIFTING FACTOR
	LSHC	T4,(G1)			;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,ERR0			;IF LO OF Y .NE. 0, GO TO ERR
	JUMPE	T4,YINT			;IF HI OF Y = 0, GO TO YINT
ERR0:	LERR	(LIB,%,<DEXP3: negative ** non-integer;  ABS(base) used instead of base>)
	JRST	CONT			;GO TO CONT

X0:	JUMPG	T2,RET1			;IF Y IS NOT GT 0
	LERR	(LIB,%,<DEXP3: base is 0 and exponent is le 0; result = infinity>)
	HRLOI	T0,377777		;RESULT = INFINITY
	HRLOI	T1,377777
	POP	P,T3			;RESTORE ACCUMULATORS
	POP	P,T2
	GOODBY	(1)			;RETURN

XOK:	PUSH	P,T4			;SAVE ACCUMULATORS
	PUSH	P,T5
	PUSH	P,G1
	PUSH	P,G2
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,G3
	PUSH	P,G4
	MOVE	T4,T0			;OBTAIN THE EXPONENT
	ASH	T4,-33			;SHIFT MANTISSA OFF
	SUBI	T4,200			;SUBTRACT 128 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	G1,T0			;SAVE A COPY OF F
	DFAD	G1,A1(T5)		;F+A1(P+1)
	DFSB	T0,A1(T5)		;F-A1(P+1)
	SUBI	T5,2			;FORM P + 1
	ASH	T5,-1			;(P+1)/2
	DFSB	T0,A2(T5)		;Z=(F-A1(P+1))-A2((P+1)/2)
	DFDV	T0,G1			;/(F+A1(P+1))
	FSC	T0,1			;
	DMOVE	G1,T0			;SAVE A COPY OF Z
	DFMP	G1,G1			;FORM Z**2
	DMOVE	G3,G1			;SAVE A COPY OF Z**2
	DFMP	G3,RP4			;R(Z)=RP4*Z**2
	DFAD	G3,RP3			;+RP3
	DFMP	G3,G1			;*Z**2
	DFAD	G3,RP2			;+RP2
	DFMP	G3,G1			;*Z**2
	DFAD	G3,RP1			;+RP1
	DFMP	G3,G1			;*Z**2
	DFMP	G3,T0			;*Z
	DFAD	G3,T0			;+Z
	DMOVE	G1,G3			;SAVE A COPY OF R(Z)
	DFMP	G3,C			;U2 = R(Z) * C
	DFAD	G3,G1			;+ R(Z)

	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,			;
	FSC	T4,-4			;
	DMOVE	T0,T2			;SAVE A COPY OF Y
	JUMPGE	T2,YPOS			;IF Y IS NEGATIVE
	  DMOVN	T0,T0			;NEGATE IT
YPOS:	MOVE	G1,T0			;GET COPY OF HIGH ORDER OF Y
	ASH	G1,-33			;SHIFT OFF MANTISSA
	CAIGE	G1,227			;IF THE EXPONENT IS LESS THAN 227
	  JRST	SMALLY			;GO TO SMALLY
	CAIE	G1,227			;IF THE EXPONENT IS > 227
	  JRST	LARGEY			;GO TO LARGEY
	SETZ	T1,			;ZERO SECOND WORD
	  JRST SGNCHK			;GO TO SGNCHK
SMALLY:	SUBI	G1,175			;GET SHIFTING FACTOR
	SETZ	T1,			;SET Y1 TO 0
	JUMPGE	G1,GETY1			;IF G1 IS .GE. 0
	SETZ	T0,			;THEN GO TO GETY1
	JRST	GETY2			;
GETY1:	AND	T0,MASK(G1)			;GET Y1
	JRST	SGNCHK			;GO TO CHECK SIGN
LARGEY:	CAIL	G1,272			;IF EXPONENT IS .GE. 272
	  JRST	SGNCHK			;GO TO SGNCHK
	SUBI	G1,230			;GET SHIFTING FACTOR
	AND	T1,REDMSK(G1)			;GET LOW ORDER PART OF Y1
SGNCHK:	JUMPGE	T2,GETY2			;IF Y IS NEGATIVE
	  DMOVN	T0,T0			;NEGATE Y1
GETY2:	DMOVE	G1,T2			;SAVE A COPY OF Y
	DFSB	T2,T0			;Y2 = Y-Y1
	DFMP	T2,T4			;FORM Y2*U1
	JFCL
	DFMP	G1,G3			;FORM Y*U2
	JFCL
	DFAD	T2,G1			;W = Y*U2+Y2*U1
	DFMP	T4,T0			;FORM U1*Y1
	JFCL
	DMOVE	T0,T2			;RECONSTRUCT W
	DFAD	T0,T4
	JFCL
	CAMG	T0,BIGW			;IF W IS NOT TOO BIG
	  JRST	WOK			;GO TO WOK
	LERR	(LIB,%,DEXP3: 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
	LERR	(LIB,%,DEXP3: result underflow)
	SETZ	T0,			; RESULT = 0
	SETZ	T1,
	POP	P,G4			;RESTORE ACCUMULATORS
	POP	P,G3
	POP	P,G2
	POP	P,G1
	POP	P,T5
	POP	P,T4
	POP	P,T3
	POP	P,T2
	GOODBY	(1)			;RETURN
WOK2:	SETZ	G2,			;ZERO G2
	MOVM	G1,T2			;GET HIGH ORDER OF W
	MOVE	G3,G1			;
	ASH	G3,-33			;SHIFT OFF MANTISSA
	SUBI	G3,175			;GET SHIFTING FACTOR
	JUMPGE	G3,GETW1		;IF G3 .GE. 0
					;THEN GO TO GETW1
	SETZ	G1,			;OTHERWISE, SET W1 = 0
	JRST	GETW2			;GO TO GETW2
GETW1:	AND	G1,MASK(G3)		;GETW1
	JUMPG	T2,GETW2		;IF W IS NEGATIVE
	  DMOVN	G1,G1			;NEGATE W1
GETW2:	DFSB	T2,G1			;W2 = W-W1
	DFAD	T4,G1			;W = W1+U1*Y1
	MOVM	T0,T4			;SAVE A COPY OF ABS(W)
	MOVE	G1,T0			;
	SETZ	T1,			;ZERO T1
	ASH	G1,-33			;SHIFT OFF MANTISSA
	SUBI	G1,175			;GET SHIFTING FACTOR
	JUMPGE	G1,GTW1			;IF G1 IS .GE. 0
					;THEN GO TO GTW1
	SETZ	T0,			;OTHERWISE, SET W1 TO 0
	JRST	GTW2			;GO TO GET W2
GTW1:	AND	T0,MASK(G1)		;GET W1
	JUMPGE	T4,GTW2			;IF W IS NEGATIVE
	  DMOVN	T0,T0			;NEGATE W1
GTW2:	DFSB	T4,T0			;FORM W-W1
	DFAD	T2,T4			;W2 = W2+(W-W1)
	MOVM	T4,T2			;SAVE A COPY OF ABS(W2)
	MOVE	G1,T4			;
	SETZ	T5,			;ZERO	T5
	ASH	G1,-33			;SHIFT	OFF MANTISSA
	SUBI	G1,175			;GET SHIFTING FACTOR
	JUMPGE	G1,GW			;IF G1 IS .GE. 0
					;THEN GO TO GW
	SETZ	T4,			;OTHERWISE, SET W=0
	JRST	GW2			;GO TO GW2
GW:	AND	T4,MASK(G1)		;GET W
	JUMPGE	T2,GW2			;IF W2 IS NEGATIVE
	  DMOVN	T4,T4			;NEGATE W
GW2:	DFAD	T0,T4			;FORM W1 + W
	FSC	T0,4			;*16
	FIX	T0,T0			;IW1
	DFSB	T2,T4			;W1 = W2 - W
	JUMPLE	T2,W2POS		;IF W2 .GT. 0
	DFSB	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
	DFMP	T2,Q7			;Z = Q7*W2
	DFAD	T2,Q6			;+Q6
	DFMP	T2,T0			;*W2
	DFAD	T2,Q5			;+Q5
	DFMP	T2,T0			;*W2
	DFAD	T2,Q4			;+Q4
	DFMP	T2,T0			;*W2
	DFAD	T2,Q3			;+Q3
	DFMP	T2,T0			;*W2
	DFAD	T2,Q2			;+Q2
	DFMP	T2,T0			;*W2
	DFAD	T2,Q1			;+Q1
	DFMP	T0,T2			;*W2
	DFMP	T0,A1(T4)		;*A1(P1+1)
	DFAD	T0,A1(T4)		;+A1(P1+1)
	FSC	T0,(T5)			;ADD M1 TO THE EXP OF Z
RET3:	POP	P,G4			;RESTORE ACCUMULATORS
	POP	P,G3
RET2:	SKIPE	IFLAG			;IF Y IS ODD NEGATIVE INTEGER
	  DMOVN	T0,T0			;NEGATE RESULT
	POP	P,G2
	POP	P,G1
	POP	P,T5
	POP	P,T4
RET1:	POP	P,T3
	POP	P,T2
	GOODBY	(1)			;RETURN

SXTNTH:	DOUBLE	175400000000,000000000000	;.0625D0
BIGW:		127.0E0				;UPPER BOUND FOR WW
SMALLW:	-128.9375E0				;LOWER BOUND FOR W
RP1:	DOUBLE	175525252525,125252514430	;.833333333333332114D-1
RP2:	DOUBLE	172631463146,147404016032	;.125000000005037992D-01
RP3:	DOUBLE	170444444413,211706653423	;.223214212859242590D-2
RP4:	DOUBLE	165707437566,322576050305	;.434457756721631196D-3
C:	DOUBLE	177705243545,053405770274	;.442695040888963407D0
Q7:	DOUBLE	160764733570,344546177513	;.149288526805956082D-4
Q6:	DOUBLE	164502757270,016046163745	;.154002904409897646D-3
Q5:	DOUBLE	167535417606,065746637065	;.133335413135857847D-02
Q4:	DOUBLE	172473125333,204616630130	;.961812905951724170D-02
Q3:	DOUBLE	174706541065,300601154766	;.555041086640855953D-1
Q2:	DOUBLE	176753767577,340542316147	;.240226506959095371D0
Q1:	DOUBLE	200542710277,276434757056	;.693147180559945296D0
A1:	DOUBLE	201400000000,000000000000	;A1(I), I=1,17 =
A12:	DOUBLE	200752225750,251110331413	;2**((1-I)/16). THIS
A13:	DOUBLE	200725403067,076722207114	;TABLE IS SEARCHED 
A14:	DOUBLE	200701463367,141251234104	;TO DETERMINE P.
A15:	DOUBLE	200656423746,126551655275	;
A16:	DOUBLE	200634222140,250770220071	;
A17:	DOUBLE	200612634520,212520333267	;
A18:	DOUBLE	200572042434,372600606660	;
A19:	DOUBLE	200552023631,237635714441	;
A110:	DOUBLE	200532540767,122052051262	;
A111:	DOUBLE	200513773265,115425047073	;
A112:	DOUBLE	200475724623,004432042153	;
A113:	DOUBLE	200460337602,214333425134	;
A114:	DOUBLE	200443417233,235261070316	;
A115:	DOUBLE	200427127017,037250572672	;
A116:	DOUBLE	200413253033,076304417305	;
A117:	DOUBLE	200400000000,000000000000	;
A2:	DOUBLE	077756350307,256663004531
A22:	DOUBLE	101406261124,022156463124
A23:	DOUBLE	702043260343,345371331424
A24:	DOUBLE	676250402176,331550372175
A25:	DOUBLE	676111401243,043654640102
A26:	DOUBLE	101640442174,047535641316
A27:	DOUBLE	676017655543,202562654434
A28:	DOUBLE	100613445334,141025236276
MASK1:	000777777777				;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

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DLOG10
EXTERN	DLG10.
DLOG10=DLG10.
PRGEND
TITLE	DLOG	NATURAL LOG FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	APRIL 8, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DLOG
EXTERN	DLOG.
DLOG=DLOG.
PRGEND
TITLE	DLOG.	NATURAL LOG FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	APRIL 8, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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:

;  XMOVEI	L,ARG
;  PUSHJ	P,DLOG  
;	OR
;  PUSHJ	P,DLOG10

;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	FORPRM
	TWOSEG	400000
	SALL

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

	HELLO	(DLOG,.)	;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
	LERR	(LIB,%,DLOG or DLOG10: zero arg; result = -infinity)
	HRLZI	T0,400000	;SET RESULT TO
	HRRZI	T1,000001	;LARGE NEGATIVE NUMBER
	GOODBY	(1)
ARGN:	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
	MOVE	T2,T0		;GET COPY OF HIGH ORDER PART OF ARG
	ASH	T2,-33		;SHIFT OFF THE MANTISSA
	SUBI	T2,200		;SUBTRACT 128 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
	FLTR	T2,T2
	MOVEM	T2,N		; SAVE N
	DFAD	T4,MHALF	;ZNUM = F-.5
	DMOVE 	T2,T4		;GET COPY OF ZNUM
	FSC	T4,-1		;ZDEN = ZNUM * .5
	DFAD	T4,HALF		; + .5
	JRST	EVALRZ
FGT:	FLTR	T2,T2
	MOVEM	T2,N		; SAVE N
	DMOVE	T2,T4
	DFAD	T2,B3		;ZNUM = F - 1.0
	FSC	T4,-1		;ZDEN = F*.5
	DFAD	T4,HALF		; + .5
EVALRZ:	DFDV	T2,T4		;Z = ZNUM/ZDEN
	DMOVE 	T4,T2		;
	DFMP	T4,T4		;W = Z*Z
	DMOVE	T0,T4		;SAVE COPY OF W
	DFAD	T4,B2		; FORM B(W). B(W) = W + B2
	DFMP	T4,T0		; * W
	DFAD	T4,B1		; + B1
	DFMP	T4,T0		; * W
	DFAD	T4,B0		; + B0
	DMOVEM	T0,SAVEW	; SAVE A COPY OF W
	DFMP	T0,A2		;FORM A(W). A(W)= A2*W
	DFAD	T0,A1		; + A1
	DFMP	T0,SAVEW	; * W
	DFAD	T0,A0		; + A0
	DFDV	T0,T4		;R(Z) = A(W)/B(W)
	DFMP	T0,SAVEW	; * W
	DFMP	T0,T2		; *Z
	JFCL
	DFAD	T0,T2		; + Z
	MOVE	T2,N		;RETRIEVE N
	MOVEI	T3,0		;ZERO OUT SECOND WORD
	DFMP	T2,C1		;FORM N*C1
	MOVE	T4,N		;RETRIEVE N
	MOVEI	T5,0		;ZERO OUT SECOND WORD
	DFMP	T4,C2		;FORM N*C2
	DFAD	T0,T4		;RESULT = N*C2 + R(Z)
	DFAD	T0,T2		; + N*C1
	SKIPN	FLAG
	  DFMP	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	200543000000,000000000000	;.693359375D0
C2:	DOUBLE	613102775750,347571527443	;-2.12194440054690583D-4
C3:	DOUBLE	177674557305,111562416145	;.434294481903251828D0
A0:	DOUBLE	570377400073,123045145570	;-.641249434237455811D2
A1:	DOUBLE	205406111210,005527754477	;.163839435630215342D2
A2:	DOUBLE	577153575223,052241327104	;-.789561128874912573D0
B0:	DOUBLE	565177200130,374467610432	;-.769499321084948798D3
B1:	DOUBLE	211470020376,205223175724	;.312032220919245328D3
B2:	DOUBLE	571342517755,046030761117	;-.356679777390346462D2
B3:	DOUBLE	576400000000,000000000000	;-.1D1
HALF:	DOUBLE	200400000000,000000000000	;.5D0
MHALF:	DOUBLE	577400000000,000000000000	;-.5D0
HI:		200552023631
LOW:		237635714441
MASK1:		000777777777
MASK2:		200000000000

	RELOC					;DATA
N:		0
FLAG:		0
SAVEW:	DOUBLE	000000000000,000000000000
	RELOC
	PRGEND
TITLE	DMOD	DOUBLE PRECISION
SUBTTL	MARY PAYNE /MHP/CKS	25-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

NOSYM
ENTRY	DMOD
EXTERN	DMOD.
DMOD=DMOD.
PRGEND
TITLE	DMOD.	DOUBLE PRECISION
SUBTTL	MARY PAYNE /MHP/CKS	25-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	SALL

	T6=6			;EXTRA TEMP REGISTERS
	T7=7
	T8=10
	T9=11


	HELLO	(DMOD,.)
	DMOVEM	T2,SAVE2	;SAVE TEMP REGISTERS
	DMOVEM	T4,SAVE4
	DMOVEM	T6,SAVE6
	DMOVEM	T8,SAVE8

	DMOVE	T0,@0(L)	;GET ABS(A) IN T0-T1
	CAIGE	T0,0
	  DMOVN	T0,T0

REPEAT:	DMOVE	T2,@1(L)	;GET ABS(B) IN T2-T3
	CAIGE	T2,0
	  DMOVN	T2,T2

	CAMG	T0,T2		;CHECK IF A .LE. B
	  JRST	FAST		;YES, GO BE FAST

	DMOVEM	T2,B		;SAVE ABS(B)
	DMOVE	T8,T0		;SAVE ABS(A) IN T8-T9
				;BREAK B INTO HIGH, MIDDLE, AND LOW PARTS 
				;B IN T2/ XXXHHHHHHHHH MMMMMMMMMLLL
	MOVE	T6,T2		;GET  T6/ XXXHHHHHHHHH 000000000000
	SETZ	T7,		;
	DFSB	T2,T6		;     T2/ XXXMMMMMMMMM LLL000000000

	MOVE	T4,T2		;     T4/ XXXMMMMMMMMM 000000000000
	SETZ	T5,		;
	DFSB	T2,T4		;     T2/ XXXLLL000000 000000000000

	DFDV	T8,B		;GET A/B
	  JFCL	OVFL		;CATCH OVERFLOW

	CAML	T8,[244400000000] ;IS QUOTIENT OVER 2**35?
	  JRST	[AND T9,[777000000000] ;YES, TRUNCATE IT TO 35 BITS
		 JRST BIG]	;	TO KEEP PRODUCTS EXACT
	FIX	T9,T8		;UNDER 2**35, TRUNCATE TO NEAREST INTEGER
	MOVSI	T8,276000	;DFLOAT THE INTEGER
	DFAD	T8,[EXP 0,0]

BIG:	DFMP	T6,T8		;GET A-B*[A/B] IN T0-T1
	DFSB	T0,T6
	DFMP	T4,T8
	DFSB	T0,T4
	DFMP	T2,T8
	DFSB	T0,T2

OVCONT:	CAIGE	T0,0		;IF RESULT NEGATIVE, [A/B] WAS 1 TOO HIGH
	  DFAD	T0,B		;CORRECT THE RESULT

	CAMGE	T0,B		;IF RESULT .LT. B, OK
	  JRST	MRET
	CAMG	T0,B		;ELSE [A/B] WAS TOO LOW
	CAML	T1,B+1		;  AND MUST SUBTRACT SOME MORE MULTIPLES OF B
	  JRST	REPEAT		;GO GET T0 MOD B IN T0

MRET:	SKIPGE	@0(L)		;GIVE RESULT SIGN OF A
	  DMOVN	T0,T0

MRET1:	DMOVE	T8,SAVE8	;RESTORE REGISTERS
	DMOVE	T6,SAVE6
	DMOVE	T4,SAVE4
	DMOVE	T2,SAVE2
	POPJ	P,		;DONE
FAST:	CAMN	T0,T2		;HERE IF A HIGH .LE. B HIGH
	CAMGE	T1,T3		;HIGH WORDS EQUAL, COMPARE LOW WORDS
	  JRST	MRET		;A .LT. B, GO RETURN DMOD = A
	DFSB	T0,T2		;ELSE RETURN DMOD = A - B
	JRST	MRET


OVFL:	JUMPE	T6,RTZ		;IF B IS ZERO, DMOD IS 0

	FSC	T0,-177		;SCALE A TO AVOID OVERFLOW IN DIVIDE
	DMOVE	T8,T0		;GET A
	DFDV	T8,B		;GET A/B
	AND	T9,[777000000000] ;TRUNCATE TO 35 BITS SO PRODUCTS ARE EXACT

	DFMP	T6,T8		;GET A-B*[A/B]
	DFSB	T0,T6
	DFMP	T4,T8
	DFSB	T0,T4
	DFMP	T2,T8
	DFSB	T0,T2

	FSC	T0,177		;UNDO SCALING
	JRST	OVCONT		;CONTINUE

RTZ:	SETZB	T0,T1		;A .EQ. B, RETURN DMOD = 0
	JRST	MRET1


	RELOC			;DATA
B:	BLOCK	2		;ABS(B)
SAVE2:	BLOCK	2		;TEMP REGISTERS
SAVE4:	BLOCK	2
SAVE6:	BLOCK	2
SAVE8:	BLOCK	2
	RELOC
	PRGEND
TITLE	DCOS	SINE AND COSINE FUNCTIONS
;	        (DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DCOS
EXTERN	DCOS.
DCOS=DCOS.
PRGEND
TITLE	DSIN	SINE AND COSINE FUNCTIONS
;	        (DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DSIN
EXTERN	DSIN.
DSIN=DSIN.
PRGEND
TITLE	DSIN.	SINE AND COSINE FUNCTIONS
;	        (DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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)], THE GREATEST INTEGER IN ABS(X)
;	X2 = ABS(X) - X1
;	N = THE NEAREST INTEGER TO ABS(X), FOR SIN, OR
;				TO ABS(X) + PI/2, FOR COS

;	THEN THE REDUCED ARGUMENT F = ((X1 - N*C1) + X2) - N*C2
;		WHERE C1 = 3.1104 (OCTAL) AND C2 IS 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) <= 6746518852.  ABS(X)+PI/2, IN
;  COS(X), MUST BE LESS THAN 6746518852.  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:
 
;  XMOVEI	L,ARG
;  PUSHJ	P,DSIN
;             OR
;  PUSHJ	P,DCOS
 
;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DCOS,.)	;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
	DFAD	T0,PID2		;Y = Y+PID2
	JRST	ALG		;PROCEED TO MAIN ALGORITHM

	HELLO	(DSIN,.)	;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:	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
	DFMP	T0,ODPI		;RN = Y/PI
	JFCL
	DMOVE	T4,T0		;SAVE A COPY OF RN
	CAMG	T0,TWO26	;IF RN IS .LE. 2**26
	  JRST	LESS		;THEN GO TO LESS
	ASH	T4,-33		;SHIFT OFF MANTISSA
	SUBI	T4,275		;OBTAIN SHIFTING FACTOR
	ASH	T1,(T4)		;SHIFT SECOND WORD APPROPRIATELY
	MOVE	T5,T1		;SAVE BITS OF SECOND WORD
	ASH	T1,-1		;SHIFT OFF REMAINDER
	MOVN	T4,T4		;PREPARE TO SHIFT BACK
	ADDI	T4,1		;ADD 1 TO SHIFTING FACTOR
	ASH	T1,(T4)		;SHIFT BACK SECOND WORD
	TRNN	T5,1		;IF RIGHTMOST OF SAVED BITS IS OFF
	  JRST 	NORND		;THEN GO TO NORND
	DFAD	T0,ONE		;OTHERWISE ADD ONE TO XN
	ASH	T5,-1		;SHIFT REMAINDER OFF OF N
	ADDI	T5,1		;RND N UP
	JRST	GOTN		;GO TO GOTN
NORND:	ASH	T5,-1		;SHIFT OFF FRACTIONAL PART
	JRST	GOTN		;GO TO GOTN
LESS:	FIXR	T0,T0		;FIX RN
	MOVE	T5,T0		;SAVE N
	FLTR	T0,T0		;XN=FLOAT(N)
	MOVEI	T1,0		;ZERO SECOND WORD
GOTN:	TRNE	T5,1		;IS N ODD?
	  MOVNS	SGN		;YES,NEGATE SIGN
	SKIPE	FLAG		;IF THE COSINE IS WANTED
	  DFAD	T0,PT5		;THEN XN=XN-.5
	DMOVE	T4,T0		;SAVE A COPY OF XN
	DFMP	T0,C1		;FORM XN*C1
	DFSB	T2,T0		;GET |X| - XN * C1
	DMOVE	T0,T4		;GET A COPY OF XN
	DFMP	T4,C2		;GET XN * C2
	DFMP	T0,C3		;GET XN * C3
	DFSB	T2,T4		;F = (|X| - XN*C1) - XN*C2
	DFSB	T2,T0		;	- 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:	DFMP	T2,T2		;G = F*F
	DMOVE	T4,T2		;SAVE A COPY OF G
	DFAD	T2,Q2		;XDEN = G+Q2
	DFMP	T2,T4		;*G
	DFAD	T2,Q1		;+Q1
	DFMP	T2,T4		;*G
	DFAD	T2,Q0		;+Q0
	DMOVEM	T2,XDEN		;MOVE XDEN TO MEMORY
	DMOVE	T2,T4		;GET A COPY OF G
	DFMP	T2,RP5		;XNUM = G*RP5
	DFAD	T2,RP4		;+RP4
	DFMP	T2,T4		;*G
	DFAD	T2,RP3		;+RP3
	DFMP	T2,T4		;*G
	DFAD	T2,RP2		;+RP2
	DFMP	T2,T4		;*G
	DFDV 	T2,XDEN		;R(G) = XNUM/XDEN
	DFAD	T2,RP1		;+RP1
	DFMP	T2,T4		;*G
	DFMP	T2,T0		;*F
	DFAD	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)

PID2:	DOUBLE 	201622077325,021026430215	;PI/2
ODPI:	DOUBLE	177505746033,162344202512	;1/PI
EPS:	DOUBLE	142400000000,000000000000	;2**(-31)
C1:	DOUBLE	202622077325,020000000000	;HIGH 34 BITS OF PI
C2:	DOUBLE	140413214106,220000000000	;NEXT 31 BITS OF PI
C3:	DOUBLE	101423063050,270033407151	;C1+C2+C3=PI TO XTRA PREC
Q2: 	DOUBLE	211612731352,270761316154	;0.394924723520450141D+3
Q1:	DOUBLE	221422322352,120271537275	;0.702492288221842518D+5
Q0:	DOUBLE	227512520255,264021153121	;0.541748285645351853D+7
RP5:	DOUBLE	605161526726,356444271113	;-0.121560740596710190D-1
RP4:	DOUBLE	203422023017,204642007402	;0.428183075897778265D+01
RP3:	DOUBLE	566026406450,010567611157	;-0.489487151969463797D+03
RP2:	DOUBLE	220540546606,025275050734	;0.451456904704461990D+05
RP1:	DOUBLE	601252525252,252525252421	;-.166666666666666667D0
TWO26:	DOUBLE	233400000000,000000000000	;2**26
ONE:	DOUBLE	201400000000,000000000000	;1.0
PT5:	DOUBLE	577400000000,000000000000	;-.5
YMAX1:	241622077325				;YMAX = (PI*2**31)
YMAX2:	021026430215

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DSINH
EXTERN	DSINH.
DSINH=DSINH.
PRGEND
TITLE	DSINH.	HYPERBOLIC SINE FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JUNE 7, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

;DSINH(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, DSINH = Y*SIGN(X)
;       EPS <= Y <= 1, DSINH = X-X*(Z*R(Z)), WHERE Z = Y**2 AND
;       R(Z) IS GIVEN BELOW.
;       1 < Y <= 88.029678, DSINH = SIGN(X)*(EXP(Y)-EXP(-Y))/2,
;       88.029678 <= Y < 128 * LN(2)
;               DSINH = SIGN(X)*((V/2)*EXP(Y - LN V)),
;       Y >= 128 * 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 DSINH IS ABS(X) < 128 * LN(2) AND ERROR MESSAGES
;  WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE.  DSINH 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:

;  XMOVEI	L,ARG
;  PUSHJ	P,DSINH

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DSINH,.)	;ENTRY TO DSINH ROUTINE
	DMOVE	T0,@(L)		;OBTAIN X
	PUSH	P,T5		;SAVE REGISTER T5
	HRLZI	T5,201400	;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,TWOM31	;IF Y IS .LE. 2**-31
	  JRST	RET2		;GO RETURN
	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
	DFMP	T2,T2		;Z = X*X
	DMOVE	T4,T2		;SAVE A COPY OF Z
	DFAD	T4,Q2		;XDEN = Z + Q2
	DFMP	T4,T2		;*Z
	DFAD	T4,Q1		;+Q1
	DFMP	T4,T2		;*Z
	DFAD	T4,Q0		;+Q0
	DMOVEM	T2,Z		;MOVE A COPY OF Z TO MEMORY
	DFMP	T2,RP3		;XNUM = Z*RP3
	DFAD	T2,RP2		;+RP2
	DFMP	T2,Z		;*Z
	DFAD	T2,RP1		;+RP1
	DFMP	T2,Z		;*Z
	DFAD	T2,RP0		; +RP0
	DFDV	T2,T4		;R(Z) = XNUM/XDEN
	DFMP	T2,Z		;*Z
	DFMP	T2,TEMP		;*(-X)
	DFAD	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	ERRD		;GO TO ERRD
	CAMGE	T1,YMAXLO	;HI OF Y = YMAXHI
	  JRST	GETW		;IF LO OF Y .LE. YMAXLO, GO TO GETW
ERRD:	HRLOI	T0,377777	;SET RESULT TO
	SETO	T1,		;MACHINE INFINITY
	  JUMPG	T5,DSNH		;IF X IS NEGATIVE
	  DMOVN	T0,T0		;RESULT = -RESULT
DSNH:	LERR	(LIB,%,DSINH: result overflow)
	POP	P,T5		;RESTORE ACCUMULATOR
	GOODBY	(1)		;RETURN
GETW:	DFAD	T0,LNV		;W = Y-LNV
	DMOVEM	T0,TEMP		;MOVE W TO A TEMP REGISTER
	FUNCT	DEXP.,<TEMP>	;Z = EXP(W)
	DMOVEM	T0,TEMP		;SAVE A COPY OF Z
	FMP	T0,CON1		;RESULT = Z*CON1
	SETZ	T1,		;ZERO SECOND WORD
	DFAD	T0,TEMP		;+Z
	  JFCL	ERRD		;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	DEXP.,<TEMP>	;GET ITS EXP
	FSC	T0,-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	DEXP.,<TEMP>	;Z = DEXP(Y)
	DMOVN	T2,T0		;SAVE A COPY OF Z
	MOVEM	T5,SGNFLG	;MOVE FLAG TO MEMORY
	HRLZI	T4,201400	;SET T4 TO 1.0
	SETZ	T5,		;CLEAR SECOND WORD
	DFDV	T4,T2		;1/Z
	DFAD	T0,T4		;+Z
	FSC	T0,-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	223527442325,224632030652	;.35181283430177117881D+6
RP1:	DOUBLE	216551270255,245012217565	;.11563521196851768270D+5
RP2:	DOUBLE	210507410130,341337112321	;.16375798202630751372D+3
RP3:	DOUBLE	200624234756,033744032643	;.78966127417357099479D+0
Q0:	DOUBLE	551376246137,720314355210	;-.21108770058106271242D+7
Q1:	DOUBLE	220432412710,355457306744	;.36162723109421836460D+5
Q2:	DOUBLE	566352207437,615500015770	;-.27773523119650701667D+3
ONE:	201400000000				;1.0
BIGX:		207540074635			;88.029691
YMAXHI:		207542710277			;88.722839111672999604
YMAXLO:		276434757153			;
TWOM31:		142400000000			;2**-31
LNV:	DOUBLE	577235067500,101343000000	;-.693147180559947174D0
CON1:		120414520000			;.186417721E-14
TWENT2:		205540000000			;22

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

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DSQRT
EXTERN	DSQRT.
DSQRT=DSQRT.
PRGEND
TITLE	DSQRT.	SQUARE ROOT FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	MARCH 19, 1979 

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

;DSQRT(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:

;  XMOVEI	L,ARG
;  PUSHJ	P,DSQRT

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DSQRT,.)		;ENTRY TO DSQRT ROUTINE
	DMOVE	T0,@(L)			;GET DP ARGUMENT
	JUMPG	T0,DSQRTP		;ARGUMENT IS GREATER THAN 0
	JUMPE	T0,DSQRT4		;ARGUMENT IS ZERO
	LERR	(LIB,%,DSQRT: negative arg; result = SQRT(ABS(arg)))

	DMOVN	T0,T0			;ARG = -ARG
DSQRTP:	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	MOVE	T2,T0			;GET SPARE COPY OF HIGH ORDER BITS
	MOVE	T3,T0			;GET SPARE COPY OF HIGH ORDER BITS
	MOVE	T4,T0			;GET SPARE COPY OF HIGH ORDER BITS
	LSH	T3,-1
	TLZE	T3,400			;WAS EXPONENT ODD?
	  JRST	DSQRT2			;YES, GO TO DSQRT2
	ADD	T3,[XWD 267,607000]	;NO,COMPUTE LINEAR APPROXIMATION
	FMPRI	T3,301454		;RESCALE
	JRST	DSQRT3
DSQRT2:	ADD	T3,[XWD 267,607000]	;ODD EXPONENT, COMPUTE LINEAR APPROXIMATION
	FMPRI	T3,301650		;RESCALE
DSQRT3:	FDV	T4,T3			;ORIGINAL ARG / INITIAL GUESS
	FAD	T4,T3			;AVERAGE THEM
	FSC	T4,-1			;
	FDV	T2,T4			;ORIGINAL ARG / SECOND GUESS
	FAD	T2,T4			;AVERAGE THEM
	FSC	T2,-1
	MOVEI	T3,0			;LOW HALF OF 1ST APPROX. IS 0
	DMOVE	T4,T0			;GET ARG BACK
	DFDV	T4,T2			;ORIGINAL ARG / THIRD GUESS
	DFAD	T2,T4			;AVERAGE THEM
	FSC	T2,-1			;
	DFDV	T0,T2			;ORIGINAL ARG / FOURTH GUESS			
	DFAD	T0,T2			;AVERAGE THEM
	FSC	T0,-1
	POP	P,T5			;RESTORE ACCUMULATORS
	POP	P,T4
	POP	P,T3
	POP	P,T2
DSQRT4:	GOODBY	(1)			;RETURN
	PRGEND
TITLE   DCOTAN	TANGENT AND COTANGENT FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.      APRIL 16, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DCOTAN
EXTERN	DCOTN.
DCOTAN=DCOTN.
PRGEND
TITLE	DTAN	TANGENT AND COTANGENT FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC       APRIL 16, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DTAN
EXTERN	DTAN.
DTAN=DTAN.
PRGEND
TITLE	DTAN.	TANGENT AND COTANGENT FUNCTIONS
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC       APRIL 16, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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**(-31)
;		= 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, Qi ARE GIVEN BELOW
;	THE RESULT IS THEN RECIPROCATED, IF NECESSARY, AND GIVEN
;	THE APPROPRIATE SIGN.
;	THE APPROXIMATION USED IS A MODIFICATION OF HART 4287.

;THE RANGE OF DEFINITION FOR DTAN IS 0 < ABS(X) < ((2**31) * (PI/2)) = 3373259426.D0

;  AND FOR DCOTAN(X), (2**(-127))*(1+(2**(-61))) < ABS(X) < ((2**31)*(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:

;  XMOVEI	L,ARG
;  PUSHJ	P,DTAN
;	OR
;  PUSHJ	P,DCOTAN

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DCOTN.,DCOTAN)	;ENTRY TO DCOTAN ROUTINE
	PUSH	P,T4		;SAVE ACCUMULATORS
	PUSH	P,T5
	SETZM	FLAG		;CLEAR FLAG FOR DCOTAN
	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:	LERR	(LIB,%,DCOTAN: ABS(arg) too large; result = zero)
	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	(DTAN,.)	;ENTRY TO TAN ROUTINE
	PUSH 	P,T4		;SAVE ACCUMULATORS
	PUSH	P,T5
	SETOM	FLAG		;SET FLAG FOR DTAN 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:	LERR	(LIB,%,DTAN or DCOTAN: result underflow)
	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 T0
	DMOVE	T0,T4		;AND SIGNED F TO T0.
	SETZM	N		;STORE 0 FOR N.
	JRST	FPOS		;BYPASS REDUCTION

REDUCE:	DFMP	T0,TWODPI	;Y*(2/PI)
	MOVE	T2,T0		;SAVE A COPY OF Y*(2/PI)
	CAMG	T2,TWO26	;IF Y*(2/PI) IS .LE. 2**26
	  JRST	LES		;THEN GO TO LES
	ASH	T2,-33		;SHIFT OFF MANTISSA
	SUBI	T2,275		;
	ASH	T1,(T2)		;SHIFT SECOND WORD APPROPRIATELY
	MOVE	T3,T1		;SAVE BITS OF SECOND WORD
	ASH	T1,-1		;SHIFT REMAINDER OFF
	MOVN	T2,T2		;PREPARE TO SHIFT BACK
	ADDI	T2,1		;
	ASH	T1,(T2)		;SHIFT BACK
	TRNN	T3,1		;IF RIGHTMOST OF SAVED BITS IS OFF
	  JRST	NORND		;THEN GO TO NORND
	  DFAD	T0,ONE		;OTHERWISE, ADD 1 TO XN
	ASH	T3,-1		;PREPARE TO LOOK AT NEXT BIT
	ADDI	T3,1		;ADD 1 TO N
	JRST	RND		;GO TO RND
NORND:	ASH	T3,-1		;SHIFT OFF FRACTIONAL PART
RND:	MOVEM	T3,N		;MOVE N TO MEMORY
	  JRST	GOTIT		;GO TO GOTIT
LES:	FIXR	T0,T0		;FIX IT
	MOVEM	T0,N		;MOVE N TO MEMORY
	FLTR	T0,T0		;FLOAT IT
	MOVEI	T1,0		;ZERO SECOND WORD
GOTIT:	JUMPLE	T4,NEXT		;IF X IS POSITIVE
	  DMOVN	T0,T0		;NEGATE XN
NEXT:	DMOVE	T2,T0		;GET A COPY OF -XN
	DFMP	T2,C1		;FORM -XN*C1
	DFAD	T4,T2 		;GET |X| - N*C1
	DMOVE	T2,T0		;GET A COPY OF -XN
	DFMP	T0,C2		;GET -XN*C2
	DFMP	T2,C3		;GET -XN*C3
	DFAD	T2,T4		;F = (|X| - XN*C1) - XN*C2
	DFAD	T2,T0		;  - 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,201400	;YES , F< EPS,
	MOVEI	T5,0		;XDEN = 1.0
	JRST	FLGCHK		;GO TO CHECK FLAG
	
NO:	DFMP	T2,T2		;NO, F .GE. EPS, SET G=F*F
	DMOVE	T4,T2		;SAVE A COPY OF G
	DFMP	T2,XP4		;XNUM = XP4 * G +
	DFAD	T2,XP3		; P3 *
	DFMP	T2,T4		; G +
	DFAD	T2,XP2		; P2 *
	DFMP	T2,T4		; G +
	DFAD	T2,XP1		; P1 *
	DFMP	T2,T4		; G *
	DFMP	T2,T0		; F +
	DFAD	T0,T2		; F
	DMOVE	T2,T4		;SAVE A COPY OF G
	DFMP	T4,Q4		;XDEN = Q4 * G +
	DFAD	T4,Q3		; Q3 *
	DFMP	T4,T2		; G +
	DFAD	T4,Q2		; Q2 *
	DFMP	T4,T2		; G +
	DFAD	T4,Q1		; Q1 *
	DFMP	T4,T2		; G +
	DFAD	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:	DFDV	T4,T0		;RESULT = XDEN/XNUM
	DMOVE	T0,T4
	JRST RET		;GO TO RET
NOVD:	DFDV	T0,T4		;RESULT = XNUM/XDEN
RET:	POP	P,T3		;RESTORE ACCUMULATORS
	POP	P,T2
	POP	P,T5
	POP	P,T4
	GOODBY	(1)

XP1:	DOUBLE	 601346652066,115432245623	;-.1372889460941120802    
XP2:	DOUBLE	 171401224404,126140152543	; .3925934686364577602E-02
XP3:	DOUBLE	 616034314474,026356056117	;-.2882482747560198194E-04
XP4:	DOUBLE	 147766720604,360442362056	; .2927308283322907641E-07
Q1:	DOUBLE	 600036052305,321342375437	;-.4706222794274454135    
Q2:	DOUBLE	 173702007252,226306364361	; .2746669449551304872E-01
Q3:	DOUBLE	 612131325464,137044675153	;-.4030063705745304384E-03
Q4:	DOUBLE	 155540343710,045637644377	; .1312960309685759549E-05
C1:	DOUBLE	201622077325,020000000000	;FIRST 34 BITS OF PI/2
C2:	DOUBLE	137413214106,220000000000	;NEXT 31 BITS OF PI/2
C3:	DOUBLE	100423063050,270033407151	;NEXT 62 BITS OF PI/2.
EPS1:		002400000000			;HIGH ORDER PART OF EPS1
YMAX1:		240622077325			;HIGH ORDER PART OF YMAX
YMAX2:		021026430215			;LOW ORDER PART OF YMAX
HIEPS:	 	142400000000			; 2**(-31)
PIO4HI:		200622077325			;HIGH ORDER PART OF PI/4
TWODPI:	DOUBLE	200505746033,162344202512
ONE:	DOUBLE	201400000000,000000000000
TWO26:	      	233400000000             	;2**26

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


;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DTANH
EXTERN	DTANH.
DTANH=DTANH.
PRGEND
TITLE	DTANH.	HYPERBOLIC TANGENT FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 > 22.1807100, TANH = 1.0*SIGN(X)
;	IF F > LN(3)/2 AND F <= 22.1807100, TANH = RESULT 1 =
;		SIGN(X)*(1 - 2/(EXP(2*F) + 1)))  
;	IF F < 2**(-31)*SQRT(3), 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:

;  XMOVEI	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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DTANH,.)		;ENTRY TO DTANH 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,LN2T32		;IF F IS .LE. LN2T32
	  JRST	CALCF			;GO TO CALCF
	HRLZI	T0,201400		;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
	FSC	T0,1			;F = F+F
	DMOVEM	T0,TEMP			;SAVE F IN A TEMP REGISTER
	FUNCT	DEXP.,<TEMP>		;EXP(2*F)
	DFAD	T0,ONE			;+1.0
	HRLZI	T2,575400		;SET T2 TO -2.0
	SETZ	T3,			;ZERO SECOND WORD
	DFDV	T2,T0			;-2/(EXP(2*F)+1)
	DMOVE	T0,T2
	DFAD	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. .806549008E-9
	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
	DFMP	T0,T0			;G = F*F
	DMOVE	T2,T0			;SAVE A COPY OF G
	DFAD	T2,Q2			;XDEN = G+Q2
      	DFMP	T2,T0			;*G
	DFAD	T2,Q1			;+Q1
	DFMP	T2,T0			;*G
	DFAD	T2,Q0			; + Q0
	DMOVE	T4,T0			;SAVE A COPY OF G
	DFMP	T4,RP2			;XNUM = G*RP2
	DFAD	T4,RP1			; + RP1
	DFMP	T4,T0			; * G
	DFAD	T4,RP0			; + RP0
	DFMP	T0,T4			; *G
	DFDV	T0,T2			;R(G) = XNUM/XDEN
	DFMP	T0,TEMP			;RESULT = F*R
	DFAD	T0,TEMP			; + F
RET:	POP	P,T4			;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	POP	P,T5
	GOODBY	(1)			;RETURN

RP0:	DOUBLE	564154513215,220360747434		;-.161341190239962281D+4
RP1:	DOUBLE	570163061227,221321463447		;-.992259296722360833D+2
RP2:	DOUBLE	577022172714,101054201376		;-.964374927772254698D0
Q0:	DOUBLE	215456407425,323513222652		;.484023570719886887D+4
Q1:	DOUBLE	214427161323,100370362574		;.22337720718962312926D+4
Q2:	DOUBLE	207702765170,172773772374		;.112744743805349493D+3
LN2T32:		205542710300				;22.1807100E0
CON1:		141673317272				;2**-32 * SQRT(3)
LN3D2:	      	200431175236             		;.549306144E0
HALF:	DOUBLE	200400000000,000000000000		;1/2
ONE:	DOUBLE	201400000000,000000000000		;1.0

	RELOC						;DATA
TEMP:	DOUBLE	0,0
	RELOC
	PRGEND
COMMENT |	Not in version 6

TITLE	DINT	DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DINT
EXTERN	DINT.
DINT=DINT.
PRGEND
TITLE	DINT.	DOUBLE PRECISION TRUNCATION TO INTEGER
SUBTTL	CHRIS SMITH/CKS		29-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DINT,.)
	PUSH	P,T2		;SAVE TEMP REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	DMOVE	T0,@(L)		;GET ABS(ARG) IN T0-T1
	JUMPGE	T0,DINT1
	DMOVN	T0,T0
DINT1:	CAMGE	T0,[201400000000] ;IS NUMBER .LT. 1?
	  JRST	[SETZB T0,T1	;YES, DINT IS 0
		 JRST R1]
	CAML	T0,[277400000000] ;IS NUMBER ALL INTEGER?
	  JRST	R		;YES, FAST RETURN

	LDB	T4,[POINT 8,T0,8] ;GET EXPONENT
	TRC	T4,-1		;GET ONE'S COMPLEMENT FOR RIGHT SHIFT
	MOVSI	T2,777000	;GET INITIAL MASK TO GET RID OF FRACTION
	SETZ	T3,
	ASHC	T2,201(T4)	;MOVE MASK OVER INTEGER PART
	AND	T0,T2		;CLEAR FRACTION BITS
	AND	T1,T3

R:	SKIPGE	@(L)		;ATTACH ORIGINAL SIGN
	  DMOVN	T0,T0

R1:	POP	P,T4		;RESTORE TEMP REGISTERS
	POP	P,T3
	POP	P,T2
	POPJ	P,		;RETURN
	PRGEND

End of code not in version 6 |
TITLE	DABS	DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DABS
EXTERN	DABS.
DABS=DABS.
PRGEND
TITLE	DABS.	DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL	/KK/DMN/TWE

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

;DOUBLE PRECISION ABSOLUTE VALUE FUNCTION
;THIS ROUTINE RETURNS THE ABSOLUTE VALUE OF A DOUBLE PRECISION
;ARGUMENT

	SEARCH	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DABS,.)	;[235] ENTRY TO DABS ROUTINE
	DMOVE	T0,@(L)		;GET ARGUMENT
	JUMPGE	T0,DABSX	;IF POSITIVE, RETURN
	DMOVN	T0,T0		;ELSE NEGATE IT
DABSX:	GOODBY	(1)		;RETURN

	PRGEND
TITLE	DMAX1	DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DMAX1
EXTERN	DMAX1.
DMAX1=DMAX1.
PRGEND
TITLE	DMAX1.	DOUBLE PRECISION MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL /TWE/CKS		29-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DMAX1,.)	;[235] ENTRY TO DMAX1 ROUTINE
	PUSH	P,T2		;SAVE AC T2
	PUSH	P,T3		;Save AC t3
	HLLZ	T3,-1(L)	;Get the F10 arg count
	DMOVE	T0,@(L)		;GET FIRST ARGUMENT
	JRST	DMAX.2		;ADDRESS OF NEXT, START CHECKING

DMAX.1:	XMOVEI	T2,@0(L)	;GET ADDRESS OF NEXT ARGUMENT
	CAMLE	T0,(T2)		;IS HIGH ORDER .GT. THIS ONE?
	  JRST	DMAX.2		;YES, GET ADDRESS OF NEXT IN LIST
	CAME	T0,(T2)		;ARE HIGH ORDER WORDS EQUAL?
	  JRST	DMAX.3		;NO, REPLACEMENT NECESSARY
	CAMGE	T1,1(T2) 	;SKIP IF PRESENT ARG IS LARGER
DMAX.3:	DMOVE	T0,(T2)		;PICK UP NEW ARG
DMAX.2:	ADDI	L,1		;Point to next arg
	AOBJN	T3,DMAX.1	;Loop for all args.
	POP	P,T3		;End of arg list
	POP	P,T2		;, restore acs
	GOODBY	(2)		;Return answer in 0 and 1
	PRGEND
TITLE	DMIN1	DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DMIN1
EXTERN	DMIN1.
DMIN1=DMIN1.
PRGEND
TITLE	DMIN1.	DOUBLE PRECISION MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	/TWE/CKS	29-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DMIN1,.)	;[235] ENTRY TO DMIN1 ROUTINE
	PUSH	P,T2		;SAVE AC T2
	PUSH	P,T3		;And T3
	HLLZ	T3,-1(L)	;Get the F10 arg count
	DMOVE	T0,@(L)		;GET FIRST ARGUMENT
	JRST	DMIN.2		;ADDRESS OF NEXT, START CHECKING

DMIN.1:	XMOVEI	T2,@0(L)	;GET ADDRESS OF NEXT ARGUMENT
	CAMGE	T0,(T2)		;IS HIGH ORDER .LT. THIS ONE?
	  JRST	DMIN.2		;YES, GET ADDRESS OF NEXT IN LIST
	CAME	T0,(T2)		;ARE HIGH ORDER WORDS EQUAL?
	  JRST	DMIN.3		;NO, REPLACEMENT NECESSARY
	CAMLE	T1,1(T2) 	;SKIP IF PRESENT ARG IS LARGER
DMIN.3:	DMOVE	T0,(T2)		;PICK UP NEW ARG
DMIN.2:	ADDI	L,1		;Point to next arg
	AOBJN	T3,DMIN.1	;End of arg list?
	POP	P,T3		;Yes, restore T3
	POP	P,T2		; and T2
	GOODBY	(2)		;Return answer in 0 and 1
	PRGEND
TITLE	DSIGN	DOUBLE PRECISION SIGN ROUTINE
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DSIGN
EXTERN	DSIGN.
DSIGN=DSIGN.
PRGEND
TITLE	DSIGN.	DOUBLE PRECISION SIGN ROUTINE
SUBTTL	/KK/DMN/CKS	29-Jan-80

;COPYRIGHT (C) 1980,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

;FROM LIB40 V.006.
;DOUBLE PRECISION TRANSFER OF SIGN
;THIS ROUTINE RETURNS ABSF(ARG1)*SIGN(ARG2)
;THE ROUTINE MAKES USE OF THE FOLLOWING TABLE:
;ARG1	ARG2	RESULT	CHANGE OF SIGN?
;+	+	+	NO
;+	-	-	YES
;-	+	+	YES
;-	-	-	NO

	SEARCH	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DSIGN,.)	;[235] ENTRY TO DSIGN ROUTINE
	DMOVE	T0,@(L)		;GET FIRST ARGUMENT
	SKIPGE	@1(L)		;THEN
	JUMPL	T0,OUT		;CHOOSE
	SKIPL	@1(L)		;THE
	JUMPGE	T0,OUT		;CORRECT
	DMOVN	T0,T0		;SIGN.
OUT:	GOODBY	(2)		;EXIT
	PRGEND
TITLE	DFLOAT	INTEGER TO DOUBLE PRECISION CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DFLOAT
EXTERN	DFLOT.
DFLOAT=DFLOT.
PRGEND
TITLE	DFLOT.	INTEGER TO DOUBLE PRECISION CONVERSION
SUBTTL	ED YOURDON/KK/VAA/TWE/DMN

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	EXTERN	DFL.0
	NOSYM
	TWOSEG	400000
	SALL

	HELLO	(DFLOT.,DFLOAT)	;[235] ENTRY TO DFLOAT ROUTINE
	MOVE	T0,@(L)		;PICK UP THE ARGUMENT
	PJRST	DFL.0		;USE DFL.0 ROUTINE
	PRGEND
TITLE	IDINT	DOUBLE PRECISION TO INTEGER CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	IDINT
EXTERN	IDINT.
IDINT=IDINT.
PRGEND
TITLE	IDINT.	DOUBLE PRECISION TO INTEGER CONVERSION
SUBTTL	TOM EGGERS/DMN

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
	EXTERN	IDF.0
	NOSYM
	TWOSEG	400000
	SALL

	HELLO	(IDINT,.)	;[235] ENTRY TO IDINT
DFIX1:	DMOVE	T0,@(L)		;GET THE ARGUMENT
	PJRST	IDF.0		;USE IDF.0 ROUTINE
	PRGEND
TITLE	DBLE	REAL TO DOUBLE PRECISION CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	DBLE
EXTERN	DBLE.
DBLE=DBLE.
PRGEND
TITLE	DBLE.	REAL TO DOUBLE PRECISION CONVERSION
SUBTTL	ED YOURDON/KK

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

;FROM LIB40 V.005.
;SINGLE PRECISION TO DOUBLE PRECISION CONVERSION
;THIS ROUTINE CONVERTS A SINGLE PRECISION ARGUMENT TO A
;DOUBLE PRECISION ANSWER WITH THE LOW ORDER WORD SET TO 0

	SEARCH	FORPRM
	NOSYM
	TWOSEG	400000
	SALL

	HELLO	(DBLE,.)	;[235] ENTRY TO DOUBLE ROUTINE
	MOVE	T0,@(L)		;GET THE ARGUMENT IN HIGH ORDER
	SETZ	T1,		;SET LOW ORDER PORTION TO ZERO
	GOODBY	(1)		;EXIT

	PRGEND
TITLE	SNGL	DOUBLE PRECISION TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM
NOSYM
ENTRY	SNGL
EXTERN	SNGL.
SNGL=SNGL.
PRGEND
TITLE	SNGL.	DOUBLE PRECISION TO REAL CONVERSION
SUBTTL	/DMN/TWE

;COPYRIGHT (C) 1973,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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.

;DOUBLE PRECISION TO SINGLE PRECISION CONVERSION FUNCTION
;THIS ROUTINE OBTAINS THE MOST SIGNIFICANT PART OF A DOUBLE 
;PRECISION ARGUMENT.

	SEARCH FORPRM
	EXTERN	SNG.0
	NOSYM
	TWOSEG	400000
	SALL

	HELLO	(SNGL,.)	;[235] ENTRY TO SNGL ROUTINE
	DMOVE	T0,@(L)		;GET THE ARGUMENT IN AC 0
	PJRST	SNG.0		;USE AC0 SINGLE ROUTINE
	PRGEND
UNIVERSAL FORDAR	GENERAL DOUBLE PRECISION ARITHMETICS
SUBTTL	TOM EGGERS/DMN/DRT/JNG/SWG	14-Aug-79

;COPYRIGHT (C) 1972,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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 MACROS FOR FORDAR UNIVERSAL FILE

IF1,<				;ONLY ONCE
;	DOUBLE PRECISION FLOAT FUNCTION "DFLOAT"

	DEFINE DFL (X) <
	XALL

	ENTRY	DFL.'X		;ENTRY POINT TO DFL.'X
	SIXBIT	/DFL.'X/
DFL.'X:	MOVEI	X+1,0		;CLEAR LOW ORDER WORD
	ASHC	X,-8		;MAKE ROOM FOR EXPONENT IN HI WORD
	TLC	X,243000	;SET EXP TO 27+8 DECIMAL
	DFAD	X,[EXP 0,0]	;NORMALIZE
	POPJ	P,		;RETURN X=THE DOUBLE PRECISION RESULT
>; END DFL


;	DOUBLE PRECISION FIX FUNCTION "IDINT"
;	DOUBLE TO INTEGER

	DEFINE IDF (X) <
	XALL

	ENTRY	IDF.'X
	SIXBIT	/IDF.'X/
IDF.'X:	PUSH	P,L		;SAVE THE SCRATCH REG
	HLRE	L,X		;GET THE EXPONENT
	ASH	L,-9		;RIGHT 8 BITS
	JUMPGE	X,IDF.XT	;JUMP IF POS.
	DMOVN	X,X		;NEGATE
	TRC	L,-1		;COMPLEMENT THE EXPONENT
IDF.XT:	TLZ	X,777000	;CLEAR THE EXPONENT
	ASHC	X,-201-^D26(L)	;CHANGE FRACTION TO INTEGER
	TLNE	L,400000	;SKIP IF POS.
	MOVN	X,X		;NEGATE
	POP	P,L		;RESTORE THE SCRATCH REG
	POPJ	P,		;RETURN X=FIXED NUMBER
>; END IDF

;	DOUBLE PRECISION TO SINGLE FUNCTION

	DEFINE SNG (X)<
	XALL

	ENTRY	SNG.'X
	SIXBIT	/SNG.'X/
SNG.'X:	JUMPL	X,SNG3		;NEGATIVE ARGUMENT?
	TLNE	X+1,(1B1)	;POSITIVE. ROUND REQUIRED?
	TRON	X,1		;YES, TRY TO ROUND BY SETTING LSB
	  POPJ	P,		;WE WON, FINISHED
	MOVE	X+1,X		;COPY HIGH PART OF ARG
	AND	X,[777000,,1]	;MAKE UNNORMALIZED LSB, SAME EXPONENT
	FAD	X,X+1		;ROUND & RENORMALIZE
	POPJ	P,

;HERE IF ARG IS NEGATIVE
SNG3:	DMOVN	X,X		;MAKE POSITIVE
	TLNE	X+1,(1B1)	;NEED ROUNDING?
	TRON	X,1		;YES, TRY TO DO IT BY SETTING LSB
	JRST	SNG4		;DONE
	MOVN	X+1,X		;MAKE RE-NEGATED COPY OF HIGH PART
	ORCA	X,[777,,-1]	;GET UNNORM NEG LSB WITH SAME EXPONENT
	FADR	X,X+1		;ROUND & NORMALIZE
	POPJ	P,

SNG4:	MOVN	X,X		;RE-NEGATE
	POPJ	P,		;EXIT
>; END SNG

>; END IF1
	PRGEND
TITLE	DFL.0

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	0
	PRGEND


TITLE	DFL.2
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	2
	PRGEND


TITLE	DFL.3
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	3
	PRGEND


TITLE	DFL.4
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	4
	PRGEND


TITLE	DFL.6
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	6
	PRGEND


TITLE	DFL.10
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	10
	PRGEND


TITLE	DFL.12
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	12
	PRGEND


TITLE	DFL.14
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	DFL	14
	PRGEND
TITLE	IDF.0

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	0
	PRGEND


TITLE	IDF.2
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	2
	PRGEND


TITLE	IDF.4
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	4
	PRGEND


TITLE	IDF.6
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	6
	PRGEND


TITLE	IDF.10
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	10
	PRGEND


TITLE	IDF.12
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	12
	PRGEND


TITLE	IDF.14
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	IDF	14
	PRGEND
TITLE	SNG.0

;COPYRIGHT (C) 1979,1981  BY  DIGITAL EQUIPMENT CORPORATION

;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	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	0
	PRGEND


TITLE	SNG.2
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	2
	PRGEND


TITLE	SNG.4
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	4
	PRGEND


TITLE	SNG.6
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	6
	PRGEND


TITLE	SNG.10
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	10
	PRGEND


TITLE	SNG.12
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	12
	PRGEND


TITLE	SNG.14
	SEARCH	FORPRM,FORDAR
	NOSYM
	TWOSEG	400000
	SNG	14

	PRGEND
	TITLE DDIM	Double precision 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) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

	SEARCH	FORPRM
	TWOSEG	400000
	NOSYM
	ENTRY	DDIM
	EXTERN	DDIM.
	DDIM=DDIM.
	PRGEND

	TITLE	DDIM.  Double Precision Positive Difference.
	SUBTTL	C. McCutcheon	6/26/81

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

;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 double precision.

; Passed: A1,A2
; Returns:
; 	A1-A2	if A1 .GT. A2
; 	0	if A1 .LE. A2

	SEARCH	FORPRM
	TWOSEG	400000
	SALL
	HELLO 	(DDIM,.)

	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.

	DFSB	T0,T2			;A1-A2

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

RET0:	SETZB	T0,T1			;Return zero
	JRST	RET
	PRGEND
	TITLE	DINT	Double precision truncation

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

	SEARCH	FORPRM
	TWOSEG	400000
	NOSYM
	ENTRY	DINT
	EXTERN	DINT.
	DINT=DINT.
	PRGEND

	TITLE	DINT.  Double Precision Truncation.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

;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
; 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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DINT,.)	;Enter routine


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

NZERO:	PUSH	P,T2		;Save ac's

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

	HLRZ	T2,T0		;Get exponent
	LSH	T2,-^D9		;Put rightmost
	CAIG	T2,200		;If number .LT. 1.0D0,
	 JRST	ZERO		; return zero.
	CAIL	T2,276		;If exponent .GE. 276 then there is
	 JRST	DONE		; no fraction, return the number passed.

; Now shift out the fractional part of the number.

	ASHC	T0,-276(T2)	;Shift into integer position.
	MOVN	T2,T2		;Negate exponent.
	ASHC	T0,276(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

; No fractional part, return number passed to function

DONE:	DMOVE	T0,@0(L)
	JRST	BYE

ZERO:	SETZB	T0,T1		;Return a 0.
	JRST	BYE
	PRGEND
	TITLE	DNINT	Double precision nearest whole number

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

	SEARCH	FORPRM
	TWOSEG	400000
	NOSYM
	ENTRY	DNINT
	EXTERN	DNINT.
	DNINT=DNINT.
	PRGEND

	TITLE	DNINT.  Double Precision nearest whole number.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

;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 double precision
; nearest whole number to the double precision number passed.

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

	SEARCH	FORPRM
	TWOSEG	400000

	HELLO	(DNINT,.)	;Enter routine

	DMOVE	T0,@0(L)	;Get argument to truncate.
	JUMPN	T0,NZERO	;Is number passed = 0?
	GOODBYE			;Yes.

NZERO:	PUSH	P,T2		;Save an ac

	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.
	HLRZ	T2,T0		;Get exponent
	LSH	T2,-^D9		;Put rightmost
	CAIL	T2,276		;If no fraction possible,
	 JRST	DONE		; return the number passed.

; Now shift out the fractional part of the number.

	DFAD	T0,[200400,,0
			 0,,0]	;Add .5D0 before truncation.
	HLRZ	T2,T0		;Get new exponent.
	LSH	T2,-^D9		;Put rightmost.
	ASHC	T0,-276(T2)	;Shift into integer position.
	MOVN	T2,T2		;Negate exponent.
	ASHC	T0,276(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.
	GOODBYE			;Return.

; No fractional part found, return number passed to function.

DONE:	DMOVE	T0,@0(L)
	JRST	BYE

ZERO:	SETZB	T0,T1		;Return 0.
	JRST	BYE

	PRGEND
	TITLE	DPROD	Double precision product for single precision

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

	SEARCH	FORPRM
	TWOSEG	400000
	NOSYM
	ENTRY	DPROD
	EXTERN	DPROD.
	DPROD=DPROD.
	PRGEND

	TITLE	DPROD.  Double Precision product for single precision.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

;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	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(DPROD,.)	;Enter routine
	MOVE	T1,@1(L)	;Get 2nd argument
	MOVEM	T1,MULT		;Save 2nd argument
	MOVE	T0,@0(L)	;Get 1st argument
	SETZB	T1,MULT+1	;Set both second words to 0

	DFMP	T0,MULT		;Multiply and leave result in AC0.

	GOODBYE			;Return

	RELOC
MULT:	BLOCK	2		;Multiplier
	RELOC
	PRGEND
	TITLE	IDNINT	Integer nearest whole number for double precision

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY ONLY BE USED
;  OR COPIED IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE.
;
;COPYRIGHT (C) 1981 BY DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

	SEARCH	FORPRM
	TWOSEG	400000
	NOSYM
	ENTRY	IDNINT
	EXTERN	IDNIN.
	IDNINT=IDNIN.
	PRGEND

	TITLE	IDNIN.  Integer nearest whole number for double prec.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (C) 1981  BY  DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASS.

;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 to return the nearest
; integer to the double precision number passed.

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

	SEARCH	FORPRM
	TWOSEG	400000
	SALL

	HELLO	(IDNIN.,IDNINT)	;Enter routine

	DMOVE	T0,@0(L)	;Get argument to round.
	JUMPN	T0,NZERO	;If number passed is 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

	DFAD	T0,[200400,,0	;Add 0.5D0 before truncation
			 0,,0]	; to round the number

; Now shift out the fractional part of the number.

	HLRZ	T2,T0		;Get exponent.
	LSH	T2,-^D9		;Put rightmost.
	CAILE	T2,243		;If not representable, (overflow imminent)
	 JRST	DONE		; return largest integer

	TLZ	T0,777000	;Eliminate exponent.
	ASHC	T0,-233(T2)	;Shift into integer position, leaving
				; 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
	GOODBYE			;Return

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

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

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

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

	END