Google
 

Trailing-Edge - PDP-10 Archives - BB-D480F-BB_1985_short - mthdbl.mac
There are 9 other files named mthdbl.mac in the archive. Click here to see a list.

	SEARCH	MTHPRM
	TV	MTHDBL	DOUBLE PRECISION ROUTINES ,2(4011)

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

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

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

	SUBTTL	REVISION HISTORY

COMMENT \

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

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.

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

3055	CKS/JLC	17-Mar-82
	Fix DEXP2 overflow/underflow message output.

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

3216	JLC	3-Sep-82
	Remove extraneous MTHDAR universal file and SEARCHes of it,
	since it was moved to MTHPRM long ago.

3221	BCM	29-Oct-82
	Moved SNGL, SNGL., and SNG.nn routines to the end of file
	MTHSNG and out of MTHDBL.  Resolves forward reference problem.

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

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

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

3244	TGS	07-Mar-83
	Retract edit 3232 to MTHDBL.  DMOD was returning inaccurate results
	in the 2d word.

3245	MRB	11-Apr-83	
	DCOTAN pops to many AC's off the stack if an result to large occurs.

3254	20-19812	MRB	13-Dec-83
	The DEXP2. routine is failing to clear the result field on an
	underflow. (This is edit 3374 to FOROT7.)

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

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

4002	JLC	3-Aug-83
	Move DEXP. to after DTAN., which calls DEXP.
	Make 1.0 a special case for DEXP3.

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

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

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

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	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) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	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) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

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

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

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

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

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

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

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

;REQUIRED (CALLED) ROUTINES:  DSQRT

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  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	MTHPRM
	SEGMENT	CODE
	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:	$LCALL	AOI
;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

	SEGMENT	DATA

TEMP:	DOUBLE	0,0
I:	0
FLAG:	0

	PRGEND
TITLE	DATAN2	ARC TAN FUNCTION
;		(DOUBLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	MAY 9, 1979

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

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

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

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

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

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

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

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

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

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



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

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

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



	HELLO (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
	 $LCALL	RUN,RET
;	  LERR	(LIB,%,<DATAN2: result underflow>,,RET)
				;IF SECOND ARG (X) POSITIVE, RESULT UNDERFLOWS
	DMOVN	T0,MPI		;ELSE RESULT IS PI WITH SIGN OF FIRST ARG
	JRST	YSIGN		;GO ATTACH SIGN AND RETURN

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

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

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

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

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

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

XHI: 	EXP    000000000000,0		; .0000000000000000000    not used
 	EXP    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)


	SEGMENT	DATA

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

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

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

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

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

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

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

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

;DCOSH(X) IS CALCULATED AS FOLLOWS

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

;       Y <= 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	MTHPRM
	SEGMENT	CODE
	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
	$LCALL	ROV
;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)

	SEGMENT	DATA

TEMP:	DOUBLE	0,0

	PRGEND
TITLE	DEXP2.	DOUBLE ** INTEGER EXPONENTIATION
SUBTTL	MARY PAYNE/MHP/CKS

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

;	Mary Payne, July 30, 1982

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

UNDMSG:	$LCALL	RUN		;Result underflow
	SETZB	T0,T1		;[3254]
RETF:	RESFLG	T2		;[3243] reset all flags, t2-t3 are dbl wd PC
RET:	DMOVE	T2,SAVE2	;Restore T2-T4
	MOVE	T4,SAVE4
	POPJ	P,

ONE:	EXP	1.0,0		;D-floating 1.0

	SEGMENT	DATA

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

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

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

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

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


;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	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(DEXP3,.)		;ENTRY TO DEXP3. ROUTINE
	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3		
	DMOVE	T0,@(L)			;GET THE BASE
	CAMN	T0,[1.0]		;IS IT EXACTLY 1.0?
	 JUMPE	T1,POPRET		;YES. JUST RETURN EXACTLY 1.0
	DMOVE	T2,@1(L)		;GET THE EXPONENT
	SETZM	IFLAG			;SET INTEGER EXP FLAG TO 0
	JUMPG	T0,XOK			;IF X IS NOT GT 0
	JUMPE	T0,X0			;AND IF X IS NOT = 0
	PUSH	P,T4			;SAVE ACCUMULATORS
	PUSH	P,T5
	PUSH	P,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:	$LCALL	NNA
;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
	JUMPE	T2,ZERZER		;SEPARATE OUT 0**0 CASE
	$LCALL	ZNI
;LERR	(LIB,%,<DEXP3: base is 0 and exponent is le 0; result = infinity>)
	HRLOI	T0,377777		;RESULT = INFINITY
	HRLOI	T1,377777
POPRET:	POP	P,T3			;RESTORE ACCUMULATORS
	POP	P,T2
	GOODBY	(1)			;RETURN

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

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
OVFL:	$LCALL	ROV			;RESULT OVERFLOW
	HRLOI	T0,377777		; RESULT = INFINITY
	HRLOI	T1,377777
	JRST	RET3			;GO TO RET3
WOK:	CAML	T0,SMALLW		;IF W IS NOT TOO SMALL
	  JRST	WOK2			;THEN PROCEED
UNFL:	$LCALL	RUN			;RESULT UNDERFLOWN
	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:	CAIG	G1,MSKTOP		;[4002] IF OUTSIDE TABLE, -1 IS MASK
	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:	CAIG	G1,MSKTOP		;[4002] IF OUTSIDE TABLE, -1 IS MASK
	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:	CAIG	G1,MSKTOP		;[4002] IF OUTSIDE TABLE, -1 IS MASK
	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
	  JFCL	EXCPT
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

EXCPT:	JUMPG	T5,OVFL			;IF T5 > 0 GO TO OVFL
	JRST	UNFL			;OTHERWISE, GO TO UNFL

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
MSKTOP==MSK26-MASK		;TOP INDEX ALLOWABLE TO MASK TABLE

	SEGMENT	DATA

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

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

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

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

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

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

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

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

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

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

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


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

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

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

;REQUIRED (CALLED) ROUTINES:  NONE

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  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	MTHPRM
	SEGMENT	CODE
	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
	$LCALL	AZM
;LERR	(LIB,%,DLOG or DLOG10: zero arg; result = -infinity)
	HRLZI	T0,400000	;SET RESULT TO
	HRRZI	T1,000001	;LARGE NEGATIVE NUMBER
	GOODBY	(1)
ARGN:	$LCALL	NAA
;	LERR	(LIB,%,DLOG or DLOG10: negative arg;result = LOG(ABS(arg)))
	DMOVN	T0,T0		;ARG = -ARG

STRT:	PUSH	P,T2		;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	PUSH	P,T5
	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

	SEGMENT	DATA

N:		0
FLAG:		0
SAVEW:	DOUBLE	000000000000,000000000000
	PRGEND
TITLE	DMOD	DOUBLE PRECISION
SUBTTL	MARY PAYNE /MHP/CKS	07-Apr-83

;[3244] Replace DMOD routine

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

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

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

NOSYM
ENTRY	DMOD
EXTERN	DMOD.
DMOD=DMOD.
PRGEND
TITLE	DMOD.	DOUBLE PRECISION
SUBTTL	MARY PAYNE /MHP/CKS	07-Apr-83

;[3244] Replace DMOD. routine

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

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

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


	SEARCH	MTHPRM
	SEGMENT	CODE
	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

	SKIPN	@1(L)		;IF ARG2 IS ZERO
	 JRST	DMRETZ		;GO RETURN 0 WITH MSG

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

DMRETZ:	SETZB	T0,T1		;RETURN ZERO IF ARG2=0
	$LCALL	MZZ
	JRST	MRET		;GO RESTORE REGISTERS

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


	SEGMENT	DATA

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

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

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

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

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

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

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

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

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

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

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

;  NOTE THAT COS(X) = COS(-X) AND SIN(X) = -SIN(-X)
;  LET ABS(X) = N*PI + F WHERE ABS(F) <= PI/2
;  WHEN ABS(F) < EPS (SEE CONSTANTS, BELOW), SIN(F) = F
;  OTHERWISE, THE ARGUMENT REDUCTION IS:
;	LET X1 = [ABS(X)], 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	MTHPRM
	SEGMENT	CODE
	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:	$LCALL	ATZ
;LERR	(LIB,%,DSIN or DCOS: ABS(arg) too large; result = zero)
	MOVEI	T0,0		;SET RESULT
	MOVEI	T1,0		;TO ZERO
	POP	P,T2		;RESTORE ACCUMULATOR
	GOODBY	(1)		;
YOK:	PUSH	P,T3		;SAVE ACCUMULATORS
	PUSH	P,T4		;
	PUSH 	P,T5
	DMOVE	T2,T0		;SAVE A COPY OF ABS(X)
	SKIPE	FLAG		;IF THIS IS THE COS ROUTINE
	  DMOVE	T2,XDEN		;RESTORE ABS(X) TO T2
	DFMP	T0,ODPI		;RN = Y/PI
	JFCL
	DMOVE	T4,T0		;SAVE A COPY OF RN
	CAMGE	T0,TWO26	;IF RN IS .LT. 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

	SEGMENT	DATA

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

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

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

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

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

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

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

;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	MTHPRM
	SEGMENT	CODE
	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:	$LCALL	ROV
;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

	SEGMENT	DATA

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

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

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

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

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

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

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

;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	MTHPRM
	SEGMENT	CODE
	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
	$LCALL	NAA
;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) DIGITAL EQUIPMENT CORPORATION 1979, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

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

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

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

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

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

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

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

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

;	THE IDENTITIES
;	 TAN(PI/2.0-G)=1.0/TAN(G)
;	 TAN(N*PI+H)=TAN(H), -PI/2.0 < H <= PI/2.0
;	 TAN(-X)=-TAN(X)
;	 COTAN(X)=1.0/TAN(X)
;	 COTAN(-X)=-COTAN(X)
;	ARE USED TO REDUCE THE ORIGINAL PROBLEM TO A PROBLEM WITH
;	AN ARGUMENT GREATER THAN -PI/2.0 AND LESS THAN OR EQUAL TO
;	PI/2.0, PI=3.14159265358979324.
;	THEN DEFINE N AND F SO THAT
;	 X=N*PI/4.0+F, 0.0 <= F <= PI/4.0, WITH PI=3.14159265358979324.
;	THEN
;	 TAN(F) = F, IF F<2**(-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	MTHPRM
	SEGMENT	CODE
	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:	$LCALL	ROV
;LERR	(LIB,%,DCOTAN: Result overflow)
	HRLOI	T0,377777	;SET RESULT TO MACHINE
	HRLOI	T1,377777	;INFINITY
	JUMPGE	T4,ERR3		;[3245]IF ARG IS NEGATIVE
	  DMOVN	T0,T0		;RESULT = - RESULT
ERR3:	POP	P,T5		;[3245]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:	$LCALL	RUN
;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 T2
	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

	SEGMENT	DATA

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


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

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

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

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

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

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

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

;TANH(X) IS CALCULATED AS FOLLOWS

;  TANH(X) = -TANH(-X).  LET F = ABS(X)
;	IF F > 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	MTHPRM
	SEGMENT	CODE
	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

	SEGMENT	DATA

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

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

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

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

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

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

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

;DEXP(X) IS CALCULATED AS FOLLOWS

;  IF X <= -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	MTHPRM
	SEGMENT	CODE
	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:	$LCALL	ROV
;LERR	(LIB,%,DEXP: result overflow)
	HRLOI	T0,377777		;DEXP = +MACHINE INFINITY
	HRLOI	T1,377777
	GOODBY	(1)			;RETURN

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

EXP1:	JUMPE	T0,[MOVSI T0,(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 NOT, 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

	SEGMENT	DATA

SAVEN:	0
	PRGEND
TITLE	DABS	DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL	H. P. WEISS/HPW		11-DEC-73

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

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	DABS
EXTERN	DABS.
DABS=DABS.
PRGEND
TITLE	DABS.	DOUBLE PRECISION ABSOLUTE VALUE
SUBTTL	/KK/DMN/TWE

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

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1973, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

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

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1973, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

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

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1973, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	DSIGN
EXTERN	DSIGN.
DSIGN=DSIGN.
PRGEND
TITLE	DSIGN.	DOUBLE PRECISION SIGN ROUTINE
SUBTTL	/KK/DMN/CKS	29-Jan-80

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

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

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

;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	MTHPRM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1973, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	DFLOAT
EXTERN	DFLOT.
DFLOAT=DFLOT.
PRGEND
TITLE	DFLOT.	INTEGER TO DOUBLE PRECISION CONVERSION
SUBTTL	ED YOURDON/KK/VAA/TWE/DMN

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

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

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

	SEARCH	MTHPRM
	EXTERN	DFL.0
	NOSYM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1973, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	IDINT
EXTERN	IDINT.
IDINT=IDINT.
PRGEND
TITLE	IDINT.	DOUBLE PRECISION TO INTEGER CONVERSION
SUBTTL	TOM EGGERS/DMN

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

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

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

	SEARCH	MTHPRM
	EXTERN	IDF.0
	NOSYM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1973, 1985
;ALL RIGHTS RESERVED.
;
;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY  IN  ACCORDANCE  WITH  THE  TERMS  OF  SUCH LICENSE AND WITH THE
;INCLUSION OF THE ABOVE COPYRIGHT NOTICE.  THIS SOFTWARE OR ANY  OTHER
;COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY
;OTHER PERSON.  NO TITLE TO AND OWNERSHIP OF THE  SOFTWARE  IS  HEREBY
;TRANSFERRED.

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

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

SEARCH	MTHPRM
NOSYM
ENTRY	DBLE
EXTERN	DBLE.
DBLE=DBLE.
PRGEND
TITLE	DBLE.	REAL TO DOUBLE PRECISION CONVERSION
SUBTTL	ED YOURDON/KK

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

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

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

;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	MTHPRM
	NOSYM
	SEGMENT	CODE
	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	DFL.0

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

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

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

	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	0
	PRGEND


TITLE	DFL.2
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	2
	PRGEND


TITLE	DFL.3
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	3
	PRGEND


TITLE	DFL.4
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	4
	PRGEND


TITLE	DFL.6
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	6
	PRGEND


TITLE	DFL.10
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	10
	PRGEND


TITLE	DFL.12
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	12
	PRGEND


TITLE	DFL.14
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	DFL	14
	PRGEND
TITLE	IDF.0

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

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

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

	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	0
	PRGEND


TITLE	IDF.2
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	2
	PRGEND


TITLE	IDF.4
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	4
	PRGEND


TITLE	IDF.6
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	6
	PRGEND


TITLE	IDF.10
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	10
	PRGEND


TITLE	IDF.12
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	12
	PRGEND


TITLE	IDF.14
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	IDF	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) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	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
	JFCL	EXCEP			;Can overflow and underflow

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

RET0:	SETZB	T0,T1			;Return zero
	JRST	RET

EXCEP:	JUMPE	T0,UNDER
	$LCALL	ROV,RET			;Result overflow
UNDER:	$LCALL	ROV,RET			;Result underflow

	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) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

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

; Fortran Instrinsic Function to return a Double Precision
; truncation of the Double Precision number passed.

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(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) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

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

; Fortran Instrinsic Function to return 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	MTHPRM
	SEGMENT	CODE

	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) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(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.
	JFCL	EXCEP		;Over/underflow can occur

RET:	POPJ	P,

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

	SEGMENT	DATA

MULT:	BLOCK	2		;Multiplier
	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) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

	SEARCH	MTHPRM
	SEGMENT	CODE
	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) DIGITAL EQUIPMENT CORPORATION 1981, 1985
;ALL RIGHTS RESERVED.

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

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

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

; Fortran Intrinsic Function 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	MTHPRM
	SEGMENT	CODE
	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:	$LCALL	ROV
;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