Google
 

Trailing-Edge - PDP-10 Archives - FORTRAN-10_V7wLink_Feb83 - mthsng.mac
There are 9 other files named mthsng.mac in the archive. Click here to see a list.
	SEARCH	MTHPRM
	TV	MTHSNG	SINGLE PRECISION ROUTINES,7(3223)

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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

1345	EDS	16-Mar-1981	Q10-04806
	Fix ASIN and ACOS library error message.

1350	EDS	16-Mar-81	Q10-04761
	Fix TWOSEG and RELOC problems.

1405	DAW	6-Apr-81
	Support extended addressing.

1464	DAW	12-May-81
	Error messages.

1524	DAW	6-Jul-81
	Shorten error message in EXP3. so it doesn't get truncated.

1673	CDM	9-Sept-81
	Add single precision routines.
	(ANINT, NINT)

1677	JLC	10-Sep-81	Q10-6510
	Use seed and multiplier from version 5A for RAN.

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

***** Begin Version 7 *****

3024	CDM	17-Nov-81
	Equate REAL. to FLOAT. to give it integer arguments instead
	of complex.

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

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.

3223	CDM	11-Nov-82
	Replace IFX.n routines due to customer QAR.

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

\
	PRGEND
TITLE	ALOG10	LOG BASE 10 FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	ALOG10
EXTERN	ALG10.
ALOG10=ALG10.
PRGEND
TITLE	ALOG	NATURAL LOG FUNCTION
;		(SINGLE PRECISION FLOATING PRECISION)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	ALOG
EXTERN	ALOG.
ALOG=ALOG.
PRGEND
TITLE	ALOG.	NATURAL LOG FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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


;ALOG10(X) AND ALOG(X) ARE CALCULATED AS FOLLOWS

;       FOR REFERENCE SEE 'COMPUTER APPROXIMATIONS', BY HART ET. AL.
;       WILEY, 1968;  ALGORITHM # 2662.
;       SEE PAGE 227 FOR THE COEFFICIENTS AND PAGE 111 FOR THE 
;       RANGE OF VALIDITY.
;
;       FOR X CLOSE TO ONE -
;       ALOG(X) = L3*Z**7+L4*Z**5+L5*Z**3+L6*Z
;       WHERE Z = (X-1)/(X+1)
;       ALOG10(X) = ALOG(X)*ALOG10(E)

;       FOR X NOT NEAR ONE - 
;       ALOG(X) = (K-1/2)*ALOG(2)+ALOG(F*SQRT(2)) WHERE X = 2**K*F
;       ALOG(F*SQRT(2)) = L3*Z**7+L4*Z**5+L5*Z**3+L6*Z
;       WHERE L3, L4, L5, AND L6 ARE CONSTANTS AND
;       Z = (F-(SQRT(2)/2))/(F+(SQRT(2)/2))

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

;REQUIRED (CALLED) ROUTINES:  NONE

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,ALOG  
;	OR
;  PUSHJ	P,ALOG10

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(ALG10.,ALOG10)	;ENTRY TO LOG TO THE BASE 10 ROUTINE.
	MOVE	T0,@(L)		;GET X IN T0.
	JUMPE	T0,LZERO	;CHECK FOR ZERO ARG.
	MOVEM	T0,TEMP		;SAVE ARGUMENT.
	FUNCT	ALOG.,<TEMP>	;CALC THE BASE E LOG TO THE
	FMPR	T0,LOG10A	;MULTIPLY IT BY LOG10(E).
	GOODBY	(1)		;RETURN

	HELLO	(ALOG,.)	;ENTRY TO LOG TO THE BASE E ROUTINE.
	MOVE	T0,@(L)		;GET X
	JUMPG	T0,ALOGOK	;ARG IS GREATER THAN 0
	JUMPE	T0,LZERO	;CHECK FOR ZERO ARGUMENT
	$LCALL	NAA
;LERR	(LIB,%,<ALOG or ALOG10: negative arg; result=log(ABS(arg))>)
	MOVN	T0,T0		;GET |X|
ALOGOK:	CAMN	T0,ONE		;CHECK FOR 1.0 ARGUMENT
	  JRST	ZERANS		;IT IS 1.0 RETURN ZERO ANS.
	CAML	T0,HI		;IS ARG >/= SQRT(2)?
	  JRST	GEN		;YES GO TO GENERAL CASE
	CAMG	T0,LO		;IS ARG </= 1/SQRT(2)?
	  JRST	GEN		;YES GO TO GENERAL CASE
	SETZM	LS		;IGNORE EXP FOR ARG NEAR 1
	MOVE	T1,T0		;SERIES VARIABLE IS
	FADR	T0,ONE		;(ARG-1)/(ARG+1) IF ARG NEAR 1
	FSBR	T1,ONE		
	JRST	MERGE		;SKIP HANDLING OF EXP IN GENL CASE
GEN:	ASHC	T0,-33		;SEPARATE FRACTION FROM EXPONENT
	ADDI	T0,211000	;FLOAT THE EXPONENT 
	MOVS	T0,T0		;NUMBER NOW IN CORRECT FL. FORMAT
	FADR	T0,BIAS		;REMOVE BIAS + COMPENSATE FOR
				;SQRT(2) FACTOR; T0 HAS EXP - 1/2
	FMPR	T0,LN2		;MULT BY LN2
	MOVEM	T0,LS		;STORE FOR FUTURE - LS=(K-1/2)*LN2
	ASH	T1,-10		;SHIFT FRACTION FOR FLOATING
	TLC	T1,200000	;FLOAT THE FRACTION PART
	MOVE	T0,T1		;COPY FRACTION
	FADR	T1,L1		;SUBTRACT (SQRT(2))/2 FROM T1
	FSBR	T0,L1		; AND ADD IT TO T0
MERGE:	FDVR	T1,T0		;T1 = T1/T0
	MOVEM	T1,LZ		;STORE NEW VARIABLE IN LZ
	FMPR	T1,T1		;CALCULATE Z**2
	MOVE	T0,L3		;PICK UP FIRST CONSTANT
	FMPR	T0,T1		;MULTIPLY BY Z**2
	FADR	T0,L4		;ADD IN NEXT CONSTANT
	FMPR	T0,T1		;MULTIPLY BY Z**2
	FADR	T0,L5		;ADD IN NEXT CONSTANT
	FMPR	T0,T1		;MULTIPLY BY Z**2
	MOVE	T1,LZ		;GET Z INTO T1
	FMPR	T0,T1		;T0 NOW HAS ALL BUT FIRST
				; TERM OF SERIES APPROX.
	FSC	T1,1		;MUL Z BY 2
	FADR	T0,T1		;DO LAST ADD FOR SERIES WITH
				; ROUNDING AND OVERHANG OF
				; AT LEAST 2 BITS
	FADR	T0,LS		;ADD IN EXPONENT PART
	GOODBY	(1)
LZERO:	$LCALL	AZM
;LERR	(LIB,%,ALOG or ALOG10: zero arg; result overflow)
	MOVE	T0,MIFI		;PICK UP MINUS INFINITY
	GOODBY	(1)		;RETURN
ZERANS:	MOVEI	T0,0		;MAKE ANSWER ZERO
	GOODBY	(1)		;RETURN

LOG10A:	177674557305
ONE:	201400000000
L1:	577225754146		;-0.707106781187 = -(SQRT(2)/2)
L3:	177464164321		;.301003281
L4:	177631177674		;.39965794919
L5:	200525253320		;.666669484507
LN2:	200542710300		;0.69314718056
MIFI:	400000000001		;LARGEST NEGATIVE FLOATING NUMBER
HI:	201552023632		;IF ARG IS BETWEEN HI AND LO
LO:	200552023632		;NO SCALING BY SQRT(2) NEEDED.
BIAS:	567377000000		;-(401)/2

	SEGMENT	DATA

TEMP:	0
LS:	0
LZ:	0
	PRGEND
TITLE	AMOD	SINGLE PRECISION REMAINDER FUNCTION
SUBTTL	MARY PAYNE /MHP/CKS	25-Jan-80

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	AMOD
AMOD=AMOD.##
PRGEND
TITLE 	AMOD.	SINGLE PRECISION REMAINDER
SUBTTL  BOB HANEK/CKS

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1983

;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

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

; Bob Hanek, July 19, 1982


	HELLO (AMOD,.)

	DMOVEM	T2, SAV23	;Save registers 2, 3,
	MOVEM	T4, SAV4	;  and 4
	MOVM	T0, @0(L)	; T0 = |A|
	MOVM	T2, @1(L)	; T2 = |B|
;
; Step 1
;
	MOVE	T3, T2		; T3 = |B|
	AND	T3, [777000000000]
				; T3 = Exponent field of B (including bias)
	MOVE	T4, T0		; T4 = |A|
	SUB	T4, T3		; high 9 bits of T4 = c
	JUMPL	T4, TSTSGN	; Done if c < 0

;
; Step 2
;
STEP2:	ASHC	T3, -33		; Get c in the low bits of T4
				;  and m+200 in low bits of T3
	TLZ	T0, 777000	; T0 = I
	TLZ	T2, 777000	; T2 = J
	CAML	T0, T2		; Compare I to J
	SUB	T0, T2		; If I >= J, I <--- I - J
	JRST	STEP5
;
; Steps 3 through 6
;
STEP3:	SETZ	T1,		; T0/T1 = L = 2^35*I -- d = 35
	DIV	T0, T2		; T0 = int(L/J), T1 = L - J*int(L/J)
	MOVE	T0, T1		; T0 = L - J*int(L/J)
STEP5:	SUBI	T4, 43		; c <--- c - d
	JUMPG	T4, STEP3	; If c > 0 go to Step 3
;
; Steps 7 and 8: At this point c = r - d or r = c + d;
;
	SETZ	T1,		; T0/T1 = 2^35*I
	ASHC	T0, (T4)	; T0/T1 = 2^(d+c)*I = 2^r*I
	DIV	T0, T2		; T1 = 2^(p-m)*R
;
; Step 9 - Obtain R in floating point format and check for underflow
;
	MOVE	T0, T1		; Copy fraction into T0
	FSC	T0, (T3)	; Insert biased exponent and normalize
	JFCL	UNDER		; Can underflow

TSTSGN:	SKIPGE	@0(L)		;Remainder in T0. If A < 0
	  MOVN	T0,T0		;  negate it
RESTOR:	DMOVE	T2,SAV23	;Restore registers 2, 3
	MOVE	T4,SAV4		;  and 4
	POPJ	P,		; Return
;
; If processing continues here, the remainder has underflowed.
;
UNDER:	$LCALL	RUN		;Result underflow
	SETZ	T0,		;Store 0 for result
	JRST	RESTOR		;Restore registers and return

	SEGMENT	DATA

SAV23:	BLOCK	2
SAV4:	BLOCK	1

	PRGEND
TITLE	ACOS	ARC SINE AND ARC COSINE FUNCTIONS  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JANUARY 18, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	ACOS
EXTERN	ACOS.
ACOS=ACOS.
PRGEND
TITLE	ASIN	ARC SINE AND ARC COSINE FUNCTIONS  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JANUARY 18, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	ASIN
EXTERN	ASIN.
ASIN=ASIN.
PRGEND
TITLE	ASIN.	ARC SINE AND ARC COSINE FUNCTIONS
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.    	JANUARY 18, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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

;ASIN(X) AND ACOS(X) ARE CALCULATED AS FOLLOWS

;  LET R(Z) = Z*(P0 + Z*(P1 + Z*P2))/(Q0 + Z*(Q1 + Z))
;	(P(I) AND Q(I) ARE GIVEN BELOW)

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

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

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

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

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

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

;REQUIRED (CALLED) ROUTINES:  SQRT

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,ASIN
;	OR
;  PUSHJ	P,ACOS

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(ACOS,.)	;ENTRY TO ACOS ROUTINE
	PUSH	P,T5		;SAVE REGISTER T5
	MOVEI	T5,1		;SET ACOS FLAG
	JRST	ALG		;GO TO MAIN CODE
	
	HELLO	(ASIN,.)	;ENTRY TO ASIN ROUTINE
	PUSH	P,T5		;SAVE REGISTER T5
	MOVEI	T5,0		;SET ASIN FLAG

ALG:	MOVE 	T1,@(L)		;OBTAIN X
	MOVM	T0,T1		;OBTAIN ABS(X) = Y
	CAMG	T0,CHALF	;IF Y < .5
	  JRST	LEHF		;JUMP TO LABEL LEHF
	CAMLE	T0,CONE		;OTHERWISE, IF Y > 1,
	  JRST 	ERR0		;GO TO WRITE ERROR MESSAGE

	PUSH	P,T2		;SAVE REGISTER T2
	MOVN	T0,T0		;OBTAIN NEW Y AND Z; (-Y
	FADRI	T0,201400	;+1)
	FSC	T0,-1		;/2 = Z
	MOVE 	T2,T0		;T2 CONTAINS Z
	PUSH	P,T1		;SAVE REGISTER T1 BEFORE CALL
	MOVEM	T0,TEMP
	FUNCT	SQRT.,<TEMP>	;OBTAIN SQUARE ROOT OF Z
	POP	P,T1		;RESTORE REGISTER T1
	MOVN	T0,T0		;SQRT(Z)=-SQRT(Z)
	FSC	T0,1		;-2*SQRT(Z)=Y
	JRST	SRES		;GO TO OBTAIN RESULT

LEHF:	PUSH	P,T2		;SAVE REGISTER T2
	CAMGE	T0,CEPS		;IF Y < EPS,
	  JRST	MSIGN		;JUMP TO LABEL MSIGN
	MOVE	T2,T0		;MOVE Y TO T2
	FMPR	T2,T0		;Y*Y=Z

SRES:	PUSH	P,T3		;SAVE REGISTER T3
	PUSH	P,T4		;SAVE REGISTER T4
	MOVE	T3,T0		;MOVE Y
	MOVE 	T4,T2		;MOVE Z
	FMPR	T2,PP2		;((P2*Z
	FADR	T2,PP1		;+P1)
	FMPR	T2,T4		;*Z
	FADR	T2,PP0		;+P0)
	FMPR	T2,T4		;*Z = R(Z) NUMERATOR
	FMPR	T0,T2		;*Y
	MOVE 	T2,T4		;(Z
	FADR	T2,Q1		;+Q1)
	FMPR	T2,T4		;*Z
	FADR	T2,Q0		;+Q0 = R(Z) DENOMINATOR
	FDVR	T0,T2		;R(Z)
	FADR	T0,T3		;+Y = RESULT
	POP	P,T4		;RESTORE REGISTER T4
        POP     P,T3		;RESTORE REGISTER T3

MSIGN:  MOVM	T2,T1		;T2 CONTAINS ABS(X)
	JUMPG	T5,COS		;IF FLAG NOT ZERO, GO TO ACOS
	CAMLE	T2,CHALF	;IF ABS(X) GREATER THAN .5
	  FADR	T0,CPI2		;ADD PI/2 TO RESULT
	JUMPGE	T1,RET		;IF X IS POSITIVE, RETURN
	MOVN	T0,T0		;OTHERWISE, RESULT=-RESULT
	JRST	RET		;RETURN

COS:	JUMPL	T1,BCOEF	;IF X IS NEGATIVE, JUMP TO BCOEF
	MOVN	T0,T0		;RESULT=-RESULT
	CAMG	T2,CHALF	;IF ABS(X) <= .5
	  FADR	T0,CPI2		;ADD PI/2 TO RESULT
	JRST	RET		;RETURN

BCOEF:  MOVE	T1,CPI2		;T1 CONTAINS PI/2
	CAMLE	T2,CHALF	;IF ABS(X) > .5
	  FSC	T1,1		;CONSTANT BECOMES PI
	FADR	T0,T1		;ADD CONSTANT TO RESULT

RET:	POP	P,T2		;RESTORE REGISTER T2
	POP	P,T5		;RESTORE REGISTER T5
	GOODBY	(1)		;RETURN TO CALLING PROGRAM

ERR0:	$LCALL	AOI
;LERR	(LIB,%,ASIN or ACOS: ABS(arg) > 1.0; result = +infinity)
	HRLOI	T0,377777	;RESULT=POSITIVE MACHINE INFINITY
	POP	P,T5		;RESTORE REGISTER T5
	GOODBY	(1)		;RETURN TO CALLING PROGRAM

CONE:	201400000000		;1.0
CHALF:	200400000000		;.5
CEPS:	1.0E-8
PP0:	200441171213		;0.564915737
PP1:    600134535161		;-0.409490163
PP2:	173475014642		;1.93496723E-2
CPI:	202622077325		;3.14159265
CPI2:	201622077325		;1.57079633
Q0:	202661665706		;3.38949412
Q1:	575002216372		;-3.98220081

	SEGMENT	DATA

TEMP:	0			;TEMPORARY STORAGE USED FOR SQRT ARG
	PRGEND
TITLE	ATAN2	TWO ARGUMENT ARC TAN FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.      JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	ATAN2
EXTERN	ATAN2.
ATAN2=ATAN2.
PRGEND
TITLE	ATAN	ARC TAN FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	ATAN
EXTERN	ATAN.
ATAN=ATAN.
PRGEND
TITLE	ATAN.	ARC TAN FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	Chris Smith/CKS

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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


;ATAN(X) is computed as follows:
;
;If X < 0, compute ATAN(|X|) below, then ATAN(X) = -ATAN(-X).
;
;If X > 0, use the identity
;
;	ATAN(X) = ATAN(XHI) + ATAN(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 single precision number,
;and so Z can be calculated without loss of significance.
;
;ATAN(XHI) is found by table lookup.  It is stored as ATANHI + ATANLO to
;provide guard bits for the final addition to ATAN(Z).
;
;ATAN(Z) is evaluated by means of a polynomial approximation from Hart et al.
;(formula 4901).
;
;If X < tan(pi/32), ATAN(X) = ATAN(Z).
;If X > 1/tan(pi/32), ATAN(X) = pi/2 - ATAN(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


	HELLO (ATAN,.)		;ATAN ENTRY

	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3

	MOVE	T0,@(L)		;GET ARGUMENT X
	MOVEM	T0,SGNFLG	;SAVE ARGUMENT SIGN FOR RESULT
	MOVM	T0,T0		;GET |X| 

	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 A COPY 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
	HRRZ	T2,OFFSET-2000+24(T2) ;GET OFFSET INTO XHI TABLES

	MOVE	T1,T0		;GET A COPY OF X
	FSBR	T0,XHI(T2)	;GET X-XHI
	FMPR	T1,XHI(T2)	;    X*XHI
	FADRI	T1,(1.0)	;    1 + X*XHI
	FDVR	T0,T1		;    (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:	MOVM	T1,T0		;GET |Z|
	CAMG	T1,EPS		;IS IT SMALL ENOUGH THAT ATAN(Z) = Z?
	  JRST	SMALLX		;YES, BE FAST, AVOID UNDERFLOW
	FMPR	T1,T1		;GET Z**2
	MOVE	T3,P03		;GET P(Z**2)
	FMPR	T3,T1
	FADR	T3,P02
	FMPR	T3,T1
	FADR	T3,P01
	FMPR	T1,T3
	FMPR	T1,T0		; * Z
	FADR	T0,T1		; + Z = ATAN(Z)

SMALLX:	FAD	T0,ATANLO(T2)	;  + ATAN(XHI) LOW
	FADR	T0,ATANHI(T2)	;  + ATAN(XHI) HI   = ATAN(X)
	SKIPG	SGNFLG		;ATTACH SIGN TO RESULT
	  MOVN	T0,T0
RET:	POP	P,T3		;RETURN
	POP	P,T2
	POPJ	P,

LARGEX:	MOVSI	T1,(-1.0)	;GET -1/X
	FDVRM	T1,T0
	MOVEI	T2,PI2OFFS	;GET OFFSET OF PI/2
	JRST	CALC		;GO COMPUTE PI/2 + ATAN(-1/X)
SUBTTL ATAN2 (Y,X)


;To compute ATAN2(Y,X), let U = |Y| and V = |X|, and compute ATAN(U/V).
;Then find ATAN2(Y,X) based on the signs of Y and X as follows:
;
;	 X	 Y	 ATAN2(Y/X)
;	
;	pos	pos	  ATAN(U/V)
;	pos	neg	 -ATAN(U/V)
;	neg	pos	-(ATAN(U/V) - pi)
;	neg	neg	  ATAN(U/V) - pi
;
;The add of -pi is combined with the add of ATAN(XHI) which is the last step
;of the ATAN algorithm.
;
;The reduced argument for ATAN 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 14 significant bits
;	VLO has at most 13 significant bits
;
;and choose XHI with at most 13 significant bits.  Then VHI*XHI and VLO*XHI can
;be exactly represented as single precision numbers, and the numerator is
;
;	U - V*XHI = (U - VHI*XHI) - VLO*XHI



HELLO (ATAN2,.)			;ATAN2 ENTRY

	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3

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

	MOVM	T3,T0		;GET |Y/X|, ATAN ARG
	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	[MOVE T0,T3	;YES, DO SO
		 JRST SMALL2]
	LSHC	T2,9		;GET INDEX INTO OFFSET TABLES
	ASHC	T2,3
	HRRZ	T2,OFFSET-2000+24(T2) ;GET INDEX INTO XHI TABLES

	PUSH	P,T4		;SAVE ANOTHER REGISTER

	MOVM	T0,@0(L)	;GET |Y| = U
	MOVM	T1,@1(L)	;GET |X| = V
	MOVE	T3,T1		;GET A COPY OF V
	MOVE	T4,T1		;GET ANOTHER
	AND	T3,[777777760000] ;GET HIGH 14 BITS OF V = VHI
	FSBR	T4,T3		;GET LOW 13 BITS OF V = VLO
	FMPR	T3,XHI(T2)	;GET V*XHI = VHI * XHI
	FMPR	T4,XHI(T2)	;	    + VLO * XHI
	FSBR	T0,T3		;GET (U - VHI*XHI)
	FSBR	T0,T4		;		   - VLO*XHI
	MOVM	T3,@0(L)	;GET U
	FMPR	T3,XHI(T2)	;GET U * XHI
	FADR	T1,T3		;GET V + U*XHI
	FDVR	T0,T1		;GET (U - V*XHI) / (V + U*XHI)

	POP	P,T4		;RESTORE T4

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:	MOVN	T0,@1(L)	;GET -X/Y
	FDVR	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,%,<ATAN2: result underflow>,,RET)
				;IF SECOND ARG (X) POSITIVE, RESULT UNDERFLOWS
	MOVN	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,%,<ATAN2: both arguments are zero, result=0.0>)
	SETZ	T0,		;RETURN 0
	JRST	RET

OVER:	MOVE	T0,PI2		;OVERFLOW, RESULT IS PI/2 WITH SIGN OF FIRST ARG
YSIGN:	SKIPGE	@0(L)		;ATTACH SIGN OF FIRST ARGUMENT (Y)
	  MOVN	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.

OFFSET:	DEC	1,1,1,2,2,3,3,4,4,4,5,5,5,6,6,7,7,7,8,8,8,9,9,9
	DEC	10,10,10,10,11,11,11,12,12,12,12,12,13,13,13,13
	DEC	13,13,14,14,14,14,14,14,14,14,14,14,14,14,14

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

XHI:	EXP    000000000000		; .000000000     not used
 	EXP    175657740000		; .105453491     X1
 	EXP    176407740000		; .128875732     X2
 	EXP    176477740000		; .156219482    
 	EXP    176617640000		; .195220947    
 	EXP    176777400000		; .249755859    
 	EXP    177477540000		; .312194824    
 	EXP    177617200000		; .389892578    
 	EXP    177776300000		; .498413086    
 	EXP    200515740000		; .652221680    
 	EXP    200674040000		; .867309570    
 	EXP    201453440000		; 1.17016602    
 	EXP    201645100000		; 1.64501953    
 	EXP    202511040000		; 2.57080078    
 	EXP    203527000000		; 5.35937500     X14

ATANHI:	EXP    000000000000		; .000000000     ATAN(0)
 	EXP    175656261521		; .105065183     ATAN(X1)
 	EXP    176406373155		; .128169263     ATAN(X2)
 	EXP    176475276501		; .154966952    
 	EXP    176612661306		; .192796122    
 	EXP    176765175625		; .244748870    
 	EXP    177465675100		; .302606821    
 	EXP    177574536625		; .371762831    
 	EXP    177731362666		; .462377273    
 	EXP    200447716237		; .577935450    
 	EXP    200555632634		; .714457721    
 	EXP    200672140433		; .863649569    
 	EXP    201406227205		; 1.02459152    
 	EXP    201463117110		; 1.19982255    
 	EXP    201542714675		; 1.38632865     ATAN(X14)
PI2:	EXP    201622077325		; 1.57079633     PI/2
MPI:	EXP    575155700452		;-3.14159268     -PI
 	EXP    575173246104		;-3.03652751     -PI + ATAN(X1)
 	EXP    575176220221		;-3.01342341     -PI + ATAN(X2)
 	EXP    575201554376		;-2.98662573    
 	EXP    575206433526		;-2.94879657    
 	EXP    575215150343		;-2.89684382    
 	EXP    575224470162		;-2.83898586    
 	EXP    575235354335		;-2.76982984    
 	EXP    575251036741		;-2.67921540    
 	EXP    575267664122		;-2.56365722    
 	EXP    575311247221		;-2.42713496    
 	EXP    575334330561		;-2.27794310    
 	EXP    575361014154		;-2.11700118    
 	EXP    576016720235		;-1.94177012    
 	EXP    576076516022		;-1.75526401     -PI + ATAN(X14)
 	EXP    576155700453		;-1.57079633     -PI/2

ATANLO:	EXP    000000000000		; .000000000    
 	EXP    636247477474		;-.313207918E-09
 	EXP    636051035473		;-.428319397E-09
 	EXP    635130365140		;-.770380426E-09
 	EXP    637267032057		;-.149588704E-09
 	EXP    142516146051		; .607905122E-09
 	EXP    634037666255		;-.174675323E-08
 	EXP    636153733615		;-.367500207E-09
 	EXP    635370616245		;-.478798078E-09
 	EXP    635035566432		;-.877241210E-09
 	EXP    143723624120		; .170180781E-08
 	EXP    144654652224		; .312016779E-08
 	EXP    632252164020		;-.497345720E-08
 	EXP    145452272305		; .434176811E-08
 	EXP    145613704216		; .576086101E-08
 	EXP    143420550604		; .992093574E-09
 	EXP    147735722717		; .278181351E-07
 	EXP    150564207746		; .433374110E-07
 	EXP    147566433274		; .218018803E-07
 	EXP    147760532442		; .289103999E-07
 	EXP    150515527434		; .388444175E-07
 	EXP    150504133030		; .377392659E-07
 	EXP    147677716232		; .260713819E-07
 	EXP    147427462256		; .162747642E-07
 	EXP    147525537265		; .198887566E-07
 	EXP    147516656470		; .194903134E-07
 	EXP    147773114124		; .295199429E-07
 	EXP    147623410142		; .234877224E-07
 	EXP    150504167752		; .377458393E-07
 	EXP    147450401401		; .172587422E-07
 	EXP    147500703763		; .186778351E-07
 	EXP    634357227174		;-.992093574E-09

;COEFFICIENTS OF APPROXIMATION POLYNOMIAL ATAN(X) = X*P(X**2)

P01:	EXP    600252525261	;-.333333308    
P02:	EXP    176631445546	; .199987124    
P03:	EXP    601337626575	;-.140725380    

EPS:	EXP 163721135503	;LARGEST X WITH ATAN(X)=X
MINX:	EXP 175623327343	;TAN(PI/32)
MAXX:	EXP 204504715427	;1/TAN(PI/32)

	SEGMENT	DATA

SGNFLG:	BLOCK	1		;SIGN TO BE ATTACHED TO RESULT

	PRGEND
TITLE	COSH	HYPERBOLIC COSINE FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	COSH
EXTERN	COSH.
COSH=COSH.
PRGEND
TITLE	COSH.	HYPERBOLIC COSINE FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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

;COSH(X) IS CALCULATED AS FOLLOWS

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

;	IF X <= 88.029678, RESULT 1 IS CALCULATED, AS
;		COSH(X) = (EXP(X) + 1/EXP(X))/2.0
;	IF 88.028678 < W < 128 * LN(2)
;	RESULT 2 IS OBTAINED, AS
;		COSH(X) = (V/2)*EXP(X - LN(V))
;	IF X >= 128 * LN(2)
;		COSH(X) = +MACHINE INFINITY
;		AND AN ERROR MESSAGE IS RETURNED

;THE RANGE OF DEFINITION FOR COSH IS ABS(X)<=88.722,
;  AND ARGUMENTS OUT OF THIS RANGE WILL CAUSE AN ERROR
;  MESSAGE TO BE TYPED.  A RESULT OF + MACHINE INFINITY
;  WILL BE RETURNED

;REQUIRED (CALLED) ROUTINES:  EXP

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,COSH

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(COSH,.)	;ENTRY TO HYPERBOLIC COSINE ROUTINE.
	MOVE	T0,@(L)		;OBTAIN THE ARGUMENT
	MOVM	T1,T0		;X=ABS(ARG)		
	CAMG	T1,EPS		;IF ABS(X) <= EPS
	  JRST	SMALL		;ANSWER IS ONE
	CAMLE	T1,XMAX  	;IF X > 88.029,
	  JRST	OV88		;GO TO OV88
;				;OTHERWISE, OBTAIN RESULT 1
	MOVEM	T1,TEMP
	FUNCT	EXP.,<TEMP>     ;EXP(X)
	CAML	T1,TWO14	;IF EXP(X) .GT. 2**14
	  JRST	HALVE		;NEGLECT EXP(-X)
	MOVSI	T1,201400	;1.0
	FDVR	T1,T0		;/EXP(X).
	FADR	T0,T1		;+ EXP(-X).
HALVE:  FSC     T0,-1           ;/2.0
	GOODBY	(1)		;RETURN

OV88:	CAMGE	T1,XXMAX	;TOO LARGE?
	  JRST	EXPP		;NO.  GO TO EXPP
OVFL:	$LCALL	ROV
;LERR	(LIB,%,COSH: result overflow)
	HRLOI	T0,377777	;ANSWER = +INFINITY.
	GOODBY	(1)		;RETURN

EXPP:	FSBR	T1,LN2VE	;X-LN(V)
	MOVEM	T1,TEMP
	FUNCT	EXP.,<TEMP>     ;OBTAIN RESULT 2
	MOVE	T1,T0		;SAVE A COPY OF EXP (2 - LN2VE)
	FMPR	T0,CON1		;MULTIPLY BY (LN2VE - LN2)
	FADR	T0,T1		;SUM = (1/2)*EXP(W)
	  JFCL	OVFL		;(ROUNDING ERRORS MIGHT CAUSE OVERFLOW)
	GOODBY	(1)		;COSH RETURN

SMALL:	MOVSI	T0,201400	;RESULT IS ONE
	GOODBY	(1)		;RETURN

LN2VE:	200542714000		;LN(V)=.693161011
XMAX:	207540074620		;88.029678
XXMAX:	207542710300		;128 * LN(2)
CON1:	160720040562		;LN2VE - LN(2)

	SEGMENT	DATA

TEMP:	0			;TEMPORARY STORAGE FOR EXP CALL
TWO14:	217400000000		;2**14
EPS:	163400000000		;2**(-14)
	PRGEND
TITLE	EXP	EXPONENTIAL FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979


;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	EXP
EXTERN	EXP.
EXP=EXP.
PRGEND
TITLE	EXP.	EXPONENTIAL FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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

;EXP(X) IS CALCULATED AS FOLLOWS

;  IF X < -89.4159863, EXP = 0
;  IF X >  88.0296919, EXP = +MACHINE INFINITY
;  OTHERWISE,
;	THE ARGUMENT REDUCTION IS:
;		LET N  = THE NEAREST INTEGER TO X/LN(2)
;		THE REDUCED ARGUMENT IS  G = X-N*LN(2)
;		THE CALCULATION IS:
;		EXP = R(G)*2**(N+1)
;	            WHERE R(G) = 0.5 + G*P/(Q - G*P)
;	            P = P1*G**2 + 0.25
;	            Q = Q1*G**2 + 0.5
;	        P1 AND Q1 ARE GIVEN BELOW AS XP1 AND XQ1

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

;REQUIRED (CALLED) ROUTINES:  NONE

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

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,EXP

;THE ANSWER IS RETURNED IN ACCUMULATOR T0


	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(EXP,.)		;ENTRY TO EXPONENTIAL ROUTINE
	MOVE	T0,@(L)		;OBTAIN THE ARGUMENT X
	CAMGE	T0,E77		;IF X<  -89.415...
	  JRST	OUT2		;GO TO EXIT.
	CAMG	T0,E7		;IF X LE 88.029...
	  JRST	EXP1		;GO TO ALGORITHM.  OTHERWISE.
	$LCALL	ROV
;LERR	(LIB,%,EXP: result overflow)
	HRLOI	T0, 377777	;EXP = +MACHINE INFINITY
	GOODBY	(1)		;RETURN

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

EXP1:	PUSH	P,T2		;SAVE ACCUMULATORS		
	PUSH	P,T3
	SETZ	T1,0		;T0,T1=X=N*LN(2)+G
	MOVE	T2,T0	
	SETZ	T3,0		;T2,T3=X
	DFMP	T2,RLN2		;T2,T3=N+G/LN(2)
	FIXR	T2,T2		;T4=N
	MOVEM	T2,TEMP		;SAVE N
	FLTR	T2,T2		;FLOAT N
	SETZ	T3,0		;T2,T3=FLOAT(N)
	DFMP	T2,LN2		;T2,T3=N*LN(2)
	DFSB	T0,T2		;T0,T1=G
	MOVE	T2,T0		;GET ABS(G)
	JUMPGE	T2,GPOS		;IF G IS NEGATIVE NEGATE G
	DMOVN	T0,T0	
GPOS:	TLNE	T1,(1B1)	;IF 2ND BIT OF 2ND WORD IS ON
	ADDI	T0,1		;ROUND G
	TLO	T0,(1B9)	;GUARD AGAINST SPILL INTO EXPONENT
	JUMPGE	T2,GOTG		;REGAIN SIGN OF G
	MOVN	T0,T0
GOTG:	MOVE	T1,T0
	FMPR	T1,T1		;G**2
	JFCL
	MOVE	T2,T1
	FMPR	T2,XP1		;P1*G**2
	JFCL
	FADRI	T2,177400	;+0.25 = P
	FMPR	T0,T2		;*G
	JFCL
	FMPR	T1,XQ1		;Q1*G**2
	JFCL
	FADRI	T1,200400	;+0.5 = Q
	FSBR	T1,T0		;Q - G*P
	FDVR	T0,T1		;G*P/(Q - G*P)
	FADRI	T0,200400	;+ 0.5 = R(G)
	MOVE	T3,TEMP		;RETRIEVE N
	FSC	T0,1(T3)	;R(G)*2**(N+1)
	POP	P,T3		;RESTORE ACCUMULATORS
	POP	P,T2

RET:	GOODBY	(1)		;RETURN

E7:	207540074636		;88.0296919
E77:	570232254037		;-89.4159863
RLN2:	DOUBLE 201561250731,112701375747	;1./LN(2)=1.44269504
LN2:	DOUBLE 200542710277,276434757260	;LN(2)
XP1:	171420514076		;.00416028863
XQ1:	174631375331		;.0499871789

	SEGMENT	DATA

TEMP:	0
 	PRGEND
TITLE	EXP1.	INTEGER ** INTEGER EXPONENENTIATION
SUBTTL	CHRIS SMITH/CKS		28-Jan-80

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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

	HELLO	(EXP1,.)
	MOVEI	T0,1		;SET RESULT TO 1
	MOVE	T1,@1(L)	;GET EXPONENT

	PUSH	P,T2		;SAVE TEMP REGISTER

	SKIPN	T2,@0(L)	;GET BASE
	  JRST	E1ZERO		;BASE IS ZERO, GO HANDLE
	JUMPG	T1,E1POS	;EXPONENT POSITIVE, GO HANDLE
	JUMPE	T1,E1RET	;EXPONENT ZERO, RETURN 1

	TRNE	T1,1		;EXPONENT NEGATIVE, IS IT ODD?
	  MOVE	T0,T2		;YES, RESULT IS -1 IF BASE IS -1
	CAME	T2,[-1]		;IS BASE -1?
	CAIN	T2,1		;OR +1?
	  JRST	E1RET		;YES, RESULT IS 1 OR -1
E1RTZ:	SETZ	T0,		;ELSE RESULT IS 0
	JRST	E1RET		;RETURN

E1LP:	IMUL	T2,T2		;SQUARE BASE
	  JFCL	E1OVFL		;CATCH OVERFLOW
E1POS:	TRNE	T1,1		;CHECK LOW BIT OF EXPONENT
	  IMUL	T0,T2		;MULTIPLY ANSWER BY BASE
	   JFCL	E1OVFL		;CATCH OVERFLOW
	LSH	T1,-1		;DIVIDE EXPONENT BY 2
	JUMPN	T1,E1LP		;HANDLE ALL BITS OF EXPONENT

E1RET:	POP	P,T2		;RESTORE TEMP REGISTERS
	POPJ	P,		;DONE

E1ZERO:	JUMPG	T1,E1RTZ	;BASE 0, RESULT IS 0 IF EXPONENT POSITIVE
	JUMPE	T1,E1ZZZ	;RESULT IS INDETERMINATE IF EXPONENT ZERO
				;ELSE FALL INTO OVERFLOW

E1OVFL:	$LCALL	ROV		;RESULT OVERFLOW
	HRLOI	T0,377777	;GUESS POSITIVE RESULT
	SKIPL	@0(L)		;CHECK SIGN OF BASE
	  JRST	E1RET		;BASE NONNEGATIVE, GO RETURN +INFINITY
	MOVE	T1,@1(L)	;BASE NEGATIVE, CHECK FOR ODD EXPONENT
	TRNE	T1,1		;ODD EXPONENT?
	  MOVSI	T0,400000	;NEGATIVE**ODD, RETURN -INFINITY
	JRST	E1RET

E1ZZZ:	$LCALL	ZZZ		;0**0 IS INDETERMINATE
	SETZ	T0,		;RETURN 0
	JRST	E1RET

	PRGEND
TITLE	EXP2.	REAL**INTEGER EXPONENTIATION
SUBTTL	MARY PAYNE/MHP/CKS

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1983

;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	(EXP2,.)
	MOVSI	T0,(1.0)	;Floating 1 to T0
	MOVM	T1,@1(L)	;|exponent| to T1
	JUMPE	T1,EXP0		;Exponent = 0 is special
	PUSH	P,T2		;Save a register
	SKIPE	T2,@0(L)	;Base to T2. If base not 0
	  JRST	STEP1		;  go to main flow
	JRST	BASE0		;Else to special code

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

	POP	P,T2		;Restore T2
	SKIPL	@1(L)		;If exponent > 0
	  POPJ	P,		;  return
	MOVSI	T1,(1.0)	;Else get reciprocal of
	FDVRM	T1,T0		;result; Underflow impossible
	  JOV	OVMSG		;  On overflow get message
	POPJ	P,		;Else return

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

BASE0:	POP	P,T2		;Base is 0, exponent is not.
	SKIPL	@1(L)		;If exponent > 0
	  JRST	ZERO		;  result is 0
	$LCALL	ROV		;Else overflow
	HRLOI	T0,377777	;Store +biggest
	POPJ	P,

ZERO:	SETZ	T0,		;Result is 0
	POPJ	P,		;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 T1 is not 0. Moreover, if T1 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, T1 = 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 T1
;makes it zero, we join the handling at OVER: for overflow/underflow
;on the MUL of T0 by T2. 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:	FMPR	T0,T2		;No over/underflow. Hence flags
				;  from square of T2 still valid
	LSH	T1,-1		;Discard low bit of exponent
	JUMPE	T1,OVER		;If T1 = 0, T0 has wrapped final
				;  result or its reciprocal
				;  which may be in range

	POP	P,T2		;Restore T2. Final product surely
	JSP	T1,.+1		;over/underflows. Get exception flags
	TLNE	T1,(1B11)	;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	OVMSG		;Else reciprocal gives overflow

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

OVER:	POP	P,T2		;Restore register
	JSP	T1,.+1		;Get exception flags
	TLNE	T1,(1B11)	;If underflow flag set
	  JRST	UNDER		;  underflow on product
	SKIPL	@1(L)		;Else, overflow on result if
	  JRST	OVMSG		;  exponent > 0. Get message
	MOVSI	T1,(1.0)	;For exponent < 0, get
	FDVRM	T1,T0		;reciprocal of wrapped overflow
	  JOV	RET		;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

OVMSG:	$LCALL	ROV		;Result overflow
	JUMPL	T0,NEGOV	;If result > 0
	HRLOI	T0,377777	;Store +BIGGEST
	POPJ	P,		;  and return

NEGOV:	MOVE	T0,[400000000001] ;If result < 0, store -BIGGEST
	POPJ	P,		;  and return

UNDMSG:	$LCALL	RUN		;result underflow
	SETZ	T0,
RET:	POPJ	P,

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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


;EXP3 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 
; ABS(X)**Y IS CALCULATED.)
;  -129.0 <= 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 EXP3 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,EXP3

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(EXP3,.)		;ENTRY TO EXP3. ROUTINE
      	MOVE	T0,@(L)			;GET THE BASE
	MOVE	T1,@1(L)		;GET THE EXPONENT
STRT:	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,T2			;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
	MOVM	T2,T1			;GET A COPY OF ABS(Y)
	MOVE	T4,T2			;PUT ABS(Y) IN T4
	SETZ	T3,			;SET T3 TO ZERO
	LSHC	T3,11			;GET EXPONENT OF ABS(Y)
	SUBI	T3,170			;GET SHIFTING FACTOR
	LSH	T2,(T3)			;SHIFT OFF EXP AND PART OF INTEGER
	TLNE	T2,400000		;IF Y IS ODD
	  SETOM	IFLAG			;SET IFLAG TO ONES
	LSH	T2,1			;SHIFT OFF REST OF INTEGER PART
	MOVM	T0,T0			;GET ABS(X) IN T0
	JUMPN	T2,ERR0			;ERROR IF Y NOT AN INTEGER.
	CAME	T0,A1			;IF X NOT -1
	  JRST	YINT			;  REJOIN MAIN FLOW.
	JRST	RET2			;ELSE RESULT = + OR - 1
	
ERR0:	$LCALL	NNA
;LERR	(LIB,%,<EXP3: negative ** non-integer; ABS(base) used instead of base>)
	JRST	CONT			;GO TO CONT
X0:	JUMPG	T1,RET1			;0 ** POSITIVE IS +1
	JUMPE	T1,ZERZER		;0 ** 0 IS AN ERROR
	$LCALL	ZNI
;LERR	(LIB,%,<EXP3: 0.0 ** negative; result = infinity>)
	HRLOI	T0,377777		;RESULT = INFINITY
	GOODBY				;RETURN
ZERZER:	$LCALL	ZZZ
;LERR	(LIB,%,<EXP3: 0.0 ** 0.0 is undefined; result = 0.0>)
	GOODBY

XOK:	CAMN	T0,A1			;IF X = 1, RESULT IS 1
	  JRST	RET1			;RETURN
	PUSH	P,T2			;SAVE ACCUMULATORS
	PUSH	P,T3
	PUSH	P,T4
YINT:	CAMN	T1,A1			;IF Y = 1.0
	  JRST	RET2			;GO TO RET2
	JUMPN	T1,CONT			;IF Y IS NOT 0, GO TO CONT
	MOVE	T0,A1			;SET RESULT TO 1.0
	JRST	RET2			;GO TO RET2
CONT:	PUSH	P,T5
	MOVE	T2,T0			;OBTAIN THE EXPONENT
	ASH	T2,-33			;SHIFT MANTISSA OFF
	SUBI	T2,200			;SUBTRACT 128 FROM EXPONENT
	MOVEM	T2,M
	MOVE	T2,T0			;OBTAIN FRACTIONAL PART
	AND	T2,MASK1		;EXTRACT MANTISSA
	IOR	T2,MASK2		;SET EXPONENT TO 0
	MOVEI	T5,1			;NP = 1
	CAMLE	T2,A1+10		;IF F GT A1(9)
	  JRST	NXT1			;THEN GO TO NXT1
	ADDI	T5,10			;OTHERWISE NP=9
NXT1:	CAMLE	T2,A1+3(T5)		;IF F GT A1(P+4)
	  JRST	NXT2			;THEN GO TO NXT2
	ADDI	T5,4			;OTHERWISE NP = NP+4
NXT2:	CAMLE	T2,A1+1(T5)		;IF F GT A1(P+2)
	  JRST	NXT3			;THEN GO TO NXT3
	ADDI	T5,2			;OTHERWISE NP=NP+2
NXT3:	MOVE	T3,T5			;NA = NP
	ADDI	T3,1			; + 1
	ASH	T3,-1			; / 2
	MOVE	T0,T2			;Z1 =
	FADR	T0,A1(T5)		;F + A1(P+1)
	MOVE	T4,T2			;Z =
	FSBR	T4,A1(T5)		;[F-A1(P+1)
	FSBR	T4,A2-1(T3)		; -A2(NA)]
	FDVR	T4,T0			;/Z1
	FSC	T4,1			;Z = Z+Z
	MOVE	T0,T4			;RZ =
	FMPR	T0,T0			;Z**2
	MOVE	T2,T0			;SAVE A COPY OF Z**2
	FMPR	T0,RP2			;*RP2
	FADR	T0,RP1			;+RP1
	FMPR	T0,T2			; * Z**2
	FMPR	T0,T4			; * Z
	FADR	T0,T4			; + Z
	MOVE	T4,T0			; U2 =
	FMPR	T0,XK			;RZ*XK
	FADR	T4,T0			;+ RZ
	MOVE	T2,M			;GET U1
	ASH	T2,4			;M*16
	SUB	T2,T5			;-NP
	FLTR	T2,T2			;FLOAT
	FSC	T2,-4			; /16 = U1
	SETZB	T3,T5			;
	DFAD	T4,T2			;U = U1+U2
	SETZ	T2,			;
	DFMP	T4,T1			;W = Y*U
	  JFCL	EXCEP			;CAN OVERFLOW OR UNDERFLOW
	CAMGE	T4,BIGW			;IF W IS NOT TOO BIG
	  JRST	WOK			;THEN PROCEED
OVFL:	$LCALL	ROV
;LERR	(LIB,%,EXP3: result overflow)
	HRLOI	T0,377777		; RESULT = INFINITY
	JRST	RET			; RETURN

EXCEP:	JUMPN	T4,OVER			;JUMP IF W = Y*U OVERFLOWS
	MOVE	T0,A1			;IF W UNDERFLOWS, RESULT IS 1
	JRST	RET3			;RETURN
OVER:	JUMPL	T4,UNDFL		;W NEGATIVE, RESULT UNDERFLOWS
	JRST	OVFL			;W POSITIVE, RESULT OVERFLOWS

WOK:	CAML	T4,SMALLW		;IF W IS NOT TOO SMALL
	  JRST	WOK2			;THEN PROCEED
UNDFL:	$LCALL	RUN
;LERR	(LIB,%,EXP3: result underflow)
	SETZ	T0,			; RESULT = 0
	JRST	RET			; RETURN

WOK2:	DMOVE	T0,T4			;SAVE A COPY OF DABS(W)
	JUMPGE	T4,WPOS
	  DMOVN	T0,T4
WPOS:	LDB	T1,[POINT 8,T0,8]	;GET BIASED EXPONENT OF ABS(W)
	SUBI	T1,175			;GET SHIFTING FACTOR
	JUMPGE	T1,GETW1		;IF T1 IS .GE. 0
	       				;THEN GO TO GETW1
	SETZB	T1,T3			;M1 = 0;  NP1+1 = 0
	PUSHJ	P,SNG.4##		;W2 = W ROUNDED TO SINGLE PRECISION
	JRST	REST
GETW1:	AND	T0,MASK(T1)		;GET W1
	JUMPG	T4,GETW2		;IF W IS NEGATIVE
	DMOVN	T4,T4			;GET |W|.
	SETZ	T1,			;ZERO EXTEND |W1|
	DFSB	T4,T0			; AND SUBTRACT FROM |W|.
	PUSHJ	P,SNG.4##		;ROUND |W2| TO SINGLE
	MOVN	T4,T4			;AND NEGATE TO GET W2
	FSC	T0,4			;GET 16*|W1|
	FIX	T3,T0			;  AND CONVERT TO INTEGER
	MOVN	T3,T3			;AND MAKE NEGATIVE
	JRST	GETSCL			;GO GET M1 AND NP1+1
GETW2:	SETZ	T1,
	DFSB	T4,T0			;W2 = W-W1
	PUSHJ	P,SNG.4##		;ROUND W2 TO SINGLE PRECISION
	FSC	T0,4			;W1*16
	FIX	T3,T0			; INT(W1*16)
GETSCL:	MOVE	T5,T3
	JUMPGE	T3,NPOS
	  ADDI	T3,17
NPOS:	ASH	T3,-4			; /16
	CAIL	T5,0			;IF N1 IS GE 0
	  ADDI	T3,1			;M1 = M1+1
	MOVE	T1,T3
	ASH	T3,4			; M1 * 16
	SUB	T3,T5			; - N1 = NP1
REST:	MOVE	T0,T4			;GET(2**W2)-1
	FMPR	T0,Q5			;Z = Q5*W2
	FADR	T0,Q4			; +Q4
	FMPR	T0,T4			;*W2
	FADR	T0,Q3			;+Q3
	FMPR	T0,T4			;*W2
	FADR	T0,Q2			;+Q2
	FMPR	T0,T4			;*W2
	FADR	T0,Q1			;+Q1
	FMPR	T0,T4			;*W2
	FMPR	T0,A1(T3)		; Z = Z*A1(NP1+1)
	FADR	T0,A1(T3)		; + A1(NP1+1)
	FSC	T0,(T1)			; ADD M1 TO THE EXP OF Z

RET:	POP	P,T5
RET2:	SKIPE	IFLAG			;IF IFLAG IS ON
	  MOVN	T0,T0			;  NEGATE RESULT
RET3:	POP	P,T4
	POP	P,T3
	POP	P,T2
RET1:	GOODBY	(1)			;RETURN

BIGW:		127.0E0			;UPPER BOUND FOR WW
SMALLW:	-129.0E0			;LOWER BOUND FOR W
XK:	0.4426950409E0			;LOG2(E)-1.0
RP1:	0.8333332862E-1			;
RP2:	0.1250648500E-1			;
Q1:	0.6931471806E+0			;USED TO DETERMINE
Q2:	0.2402265061E+0			;2**W2
Q3:	0.5550404881E-1			;
Q4:	0.9616206596E-2			;
Q5:	0.1305255159E-2			;
A1:	1.0E0				;A1(I), I=1,17 =
A12:	200752225751			;2**((1-I))/16. THIS
A13:	200725403067			;TABLE IS SEARCHED 
A14:	200701463367			;TO DETERMINE P.
A15:	200656423746			;
A16:	200634222141			;
A17:	200612634521			;
A18:	200572042435			;
A19:	200552023632			;
A110:	200532540767			;
A111:	200513773265			;
A112:	200475724623			;
A113:	200460337603			;
A114:	200443417234			;
A115:	200427127017			;
A116:	200413253033			;
A117:	0.5E0				;
A2:	633244441546			;
A22:	144605245161			;
A23:	633243741100			;
A24:	637260060666			;
A25:	144510250245			;
A26:	140443204215			;
A27:	633165304342			;
A28:	143763044173			;
MASK1:	000777777777			;MASK FOR MANTISSA
MASK2:	200000000000			;MASK FOR EXPONENT
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			;

	SEGMENT	DATA

IFLAG:	0				;ODD INT EXP  FLAG
M:	0				;
	PRGEND

TITLE   RAN	RANDOM NUMBER GENERATOR
;                    (INTEGER)
SUBTTL  IMSL, INC.    JANUARY 23, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED
;ONLY IN ACCORRANCE  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.

;SETRAN MAY BE USED AS FOLLOWS
;  THE SEED FOR RAN MAY BE INITIALIZED TO AN INTEGER VALUE IN THE
;  EXCLUSIVE RANGE (1,2147483647) BY A CALL TO SETRAN(I) WITH I = SEED.
;  SETRAN PROTECTS AGAINST AN ILLEGAL SEED BY SETTING SEED EQUAL TO
;  MOD(ABS(SEED),2147483648).
;  A CALL TO SETRAN WITH I = 0 WILL SET THE SEED TO THE DEFAULT VALUE
;  USED BY RAN. THE DEFAULT VALUE OF SEED IS 524287.
;  UPON EACH CALL TO RAN A NEW SEED IS GENERATED FOR SUBSEQUENT USE.
;RAN(I) IS CALCULATED AS FOLLOWS
;  FIRST THE NEW SEED IS CALCULATED AS
;  SEED = MOD(SEED*630360016,2147483647)
;  THE RANDOM NUMBER IS THEN CALCULATED AS
;  RAN = SEED/2147483648
;A CALL TO SAVRAN, CALL SAVRAN(I), RETURNS THE LAST SEED USED BY 
;  RAN IN I.
;

;REQUIRED (CALLED) ROUTINES:  NONE

;THE ROUTINES SETRAN AND SAVRAN HAVE THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,SETRAN
;       OR
;  PUSHJ	P,SAVRAN

;THE RESULT FOR SAVRAN IS RETURNED IN ARG

;THE ROUTINE RAN HAS THE FOLLOWING CALLING SEQUENCE:

;  PUSHJ	P,RAN

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(SETRAN)	;ENTRY TO SETRAN
	MOVM	T0,@(L)		;OBTAIN ABSOLUTE VALUE OF ARGUMENT
	AND	T0,CON4		;T0=MOD(T0,2**31)
	CAIN	T0,0		;IF T0 IS ZERO
	  MOVE	T0,JSEED	;SET T0 TO DEFAULT SEED
	MOVEM	T0,ISEED	;MOVE NEW SEED TO MEMORY
	GOODBY	(1)		;RETURN

	HELLO	(RAN)		;ENTRY TO RAN
    	MOVE	T0,ISEED	;OBTAIN CURRENT SEED
	MUL	T0,CON3		;SEED=SEED*16807
	DIV	T0,CON4		;MOD(SEED,(2**31)-1)
	MOVEM	T1,ISEED	;MOVE NEW SEED TO MEMORY
	MOVSI	T0,237000	;FLOAT SEED AND DIVIDE BY 2**31
	DFAD	T0,DZERO	;NORMALIZE RESULT
	GOODBY	(1)		;RETURN

   	HELLO	(SAVRAN)	;ENTRY TO SAVRAN
	MOVE	T0,ISEED	;OBTAIN CURRENT SEED
	MOVEM	T0,@(L)		;MOVE CURRENT SEED TO ARGUMENT
	GOODBY	(1)		;RETURN

DZERO:	EXP	0,0		;ZEROS
CON3:	^D630360016		;MULTIPLIER FROM VERSION 5A
CON4:	017777777777		;(2**31) - 1

	SEGMENT	DATA

ISEED:	^D524287		;SEED FROM VERSION 5A
JSEED:	^D524287		;SEED FROM VERSION 5A

 	PRGEND
TITLE	RANS	SHUFFLED RANDOM NUMBER GENERATOR  
;		(SINGLE PRECISION)
SUBTTL	IMSL, INC.	FEBRUARY 26, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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

;RANS IS A PRIME-MODULUS RANDOM NUMBER GENERATOR WITH SHUFFLING WHICH
;COMPUTES PSEUDO RANDOM NUMBERS AS FOLLOWS:
;  ON THE INITIAL REFERENCE TO RANS, RAN IS CALLED 128 TIMES TO 
;  GENERATE S(1),S(2),...,S(128) (UNIFORM RANDOM DEVIATES IN (0,1)) AND A
;  NEW SEED X(0). X(0) IS OBTAINED BY CALLING SAVRAN AFTER S(128) HAS BEEN
;  GENERATED.

;  THEN X(I+1) = (16807*X(I) MOD ((2**31)-1)
;       J = (X(I+1) MOD 128) +1
;       T = S(J)
;       S(J) = X(I+1)/(2**31)
;  RANS RETURNS T AS ITS FUNCTION VALUE.

;REQUIRED (CALLED) ROUTINES: RAN

;REGISTER T2 WAS SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE;

;	PUSHJ	P,RAN

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(RANS)		;ENTRY TO RANS
     	PUSH	P,T2			;
	SKIPN	T1,IOPT			;IF THIS IS NOT INITIAL ENTRY
 	  JRST	GGUW			;THEN GO GENERATE AND SHUFFLE
	SETZM	IOPT			;
	MOVE	T2,SETCNT		;SET COUNTERS

GGUB:	FUNCT	RAN,<IOPT>		;GENERATE A RANDOM NUMBER
	MOVEM	T0,WK-1(T2)
	AOBJN	T2,GGUB			;RETURN TO GGUB UNTIL 128
					;NUMBERS HAVE BEEN GENERATED

GGUW:	FUNCT	SAVRAN,<ICEED>	;GET LAST SEED GENERATED
	MUL	T0,CON3			;SEED = SEED*16807
	DIV	T0,CON4			;MOD(SEED,(2**31)-1)
	MOVEM	T1,ICEED		;MOVE NEW SEED TO MEMORY
	FUNCT	SETRAN,<ICEED>	;SET NEW SEED
	ANDI	T1,177			;SEED = MOD(SEED,128.)
	ADDI	T1,1			;J = SEED+1
	MOVE	T0,ICEED		;RETRIEVE SEED
	FLTR	T0,T0			;FLOAT SEED
	FSC	T0,-37			;RESULT = SEED/2**31
	MOVE	T2,T0			;
	MOVE	T0,WK-1(T1)		;RESULT = WK(J)
	MOVEM	T2,WK-1(T1)		;WK(J) = T2
	POP	P,T2
	GOODBY	(1)

SETCNT:	777600000001			;SET COUNTERS
CON4:	017777777777			;(2**31)-1
CON3:	40647				;16807

	SEGMENT	DATA

ICEED:	0				;
IOPT:	1				;OPTION SWITCH
WK:	BLOCK	200			;ALLOCATE WORK AREA
	PRGEND
TITLE	COS	SINE AND COSINE FUNCTIONS
;	        (SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	COS
EXTERN	COS.
COS=COS.
PRGEND
TITLE	SIN	SINE AND COSINE FUNCTIONS
;	        (SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	SIN
EXTERN	SIN.
SIN=SIN.
PRGEND
TITLE	COSD	SINE AND COSINE FUNCTIONS
;	        (SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	COSD
EXTERN	COSD.
COSD=COSD.
PRGEND
TITLE	SIND	SINE AND COSINE FUNCTIONS
;	        (SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	SIND
EXTERN	SIND.
SIND=SIND.
PRGEND
TITLE	SIN.	SINE AND COSINE FUNCTIONS
;	        (SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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:
;	N = THE NEAREST INTEGER TO ABS(X)/PI, FOR SIN, OR
;	    THIS INTEGER + 1/2 FOR COS	

;	THEN THE REDUCED ARGUMENT F = ABS(X)-N*PI

;	LET G = F**2
;	THEN R(G) = ((((R5*G + R4)*G + R3)*G + R2)*G + R1)*G
;	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 IS GIVEN BY [|X|/PI] < 2**26 FOR RADIANS, AND BY
;  [|X|/180] < 2**20 FOR DEGREES.
;  SIN(X) = COS(X) = 0.0 AND AN
;  ERROR MESSAGE WILL RESULT FOR ARGUMENTS OUT OF THAT RANGE
 
;REQUIRED (CALLED) ROUTINES:  NONE

;REGISTERS T2 AND T3 WERE SAVED, USED, AND RESTORED
 
;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:
 
;  XMOVEI	L,ARG
;  PUSHJ	P,SIN
;             OR
;  PUSHJ	P,COS
 
;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL
	HELLO	(COSD,.)
	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	MOVM	T0,@(L)		;GET DABS(X) IN T0 AND T2
	SETZ	T1,
	CAMGE	T0,K180		;IS ARG REDUCTION NECESSARY?
	  JRST	GTCOSD		;NO, SKIP IT
	DMOVE	T2,T0
	FDVRI	T2,(180.0)	;GET X / 180
	CAML	T2,[225400000000] ;MORE THAN 20 BITS?
	  JRST	OV		;YES, CAN'T RETURN ACCURATE RESULT
	FIX	T4,T2		;GET N
	FLTR	T2,T4
	FADRI	T2,(0.5)
	FMPR	T2,K180
	FSBR	T0,T2		;GET |X| - ([|X|/180] + 1/2) * 180
	DFMP	T0,PI180	;CONVERT TO RADIANS
	AOJA	T4,GETABS	;CONTINUE BELOW

GTCOSD:	FSBRI	T0,(90.0)	;OFFSET ARG TO [-90,90]
	MOVEI	T4,1		;SET SIGN TO MINUS
	JRST	GETRAD		;GO CONVERT TO RADIANS

	HELLO	(SIND,.)
	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	MOVM	T0,@(L)		;GET DABS(X) IN T0 AND T2
	SETZB	T1,T4		;CLEAR SIGN FLAG
	CAMGE	T0,[90.0]	;IS ARG REDUCTION NECESSARY?
	  JRST	GTSIND		;NO, SKIP IT
	DMOVE	T2,T0
	FDVRI	T2,(180.0)	;GET X / 180
	CAML	T2,[225400000000] ;MORE THAN 20 BITS?
	  JRST	OV		;YES, CAN'T RETURN ACCURATE RESULT
	FIXR	T4,T2		;GET N
	FLTR	T2,T4
	FMPR	T2,K180
	FSBR	T0,T2		;GET X - [X/180]*180

GTSIND:	SKIPGE	@(L)		;CHECK ARG SIGN
	  TRC	T4,1		;NEGATIVE, COMPLEMENT RESULT SIGN

GETRAD:	DFMP	T0,PI180	;CONVERT TO RADIANS
	JRST	GETABS		;CONTINUE BELOW
	HELLO	(COS,.)
	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	MOVM	T0,@(L)		;GET DABS(X) IN T0-T1
	SETZ	T1,
	CAMG	T0,PI		;IS ARG REDUCTION NECESSARY?
	  JRST	GETCOS		;NO, SKIP IT
	DMOVE	T2,T0		;COPY X
	DFMP	T2,ONEDPI	;GET X / PI
	CAML	T2,[233400000000] ;MORE THAN 26 BITS?
	  JRST	OV		;YES, CAN'T RETURN ACCURATE RESULT
	FIX	T4,T2		;GET [|X|/PI]
	FLTR	T2,T4
	FADRI	T2,(0.5)	;GET [|X|/PI] + 1/2
	SETZ	T3,		;CLEAR LOW WORD
	DFMP	T2,PI		;GET |X| - ([|X|/PI] + 1/2) * PI
	DFSB	T0,T2
	AOJA	T4,GETABS	;CONTINUE ARG REDUCTION BELOW

GETCOS:	DFSB	T0,PI2		;OFFSET ARG TO [-PI/2,PI/2]
	MOVEI	T4,1		;SET SIGN TO MINUS
	JRST	GETABS		;CONTINUE BELOW

	HELLO	(SIN,.)
	PUSH	P,T2		;SAVE REGISTERS
	PUSH	P,T3
	PUSH	P,T4
	MOVM	T0,@(L)		;GET DABS(X) IN T0-T1
	SETZB	T1,T4		;CLEAR SIGN FLAG
	CAMG	T0,PI2		;IS ARG REDUCTION NECESSARY?
	  JRST	GETSIN		;NO, SKIP IT
	DMOVE	T2,T0		;GET DBLE(X) IN T2
	DFMP	T2,ONEDPI	;GET X / PI
	CAML	T2,[233400000000] ;MORE THAN 26 BITS?
	  JRST	OV		;YES, CAN'T RETURN ACCURATE RESULT
	FIXR	T4,T2		;GET [X/PI] = N
	FLTR	T2,T4
	SETZ	T3,		;CLEAR LOW WORD
	DFMP	T2,PI		;GET X - [X/PI] * PI
	DFSB	T0,T2

GETSIN:	SKIPGE	@(L)		;CHECK SIGN OF ARG
	  TRC	T4,1		;NEGATIVE, COMPLEMENT RESULT SIGN

GETABS:	JUMPGE	T0,RND		;IF X POSITIVE, GO ROUND IT
	DMOVN	T0,T0		;GET ABS(X)
	TRC	T4,1		;COMPLEMENT SIGN OF RESULT

RND:	TLNN	T1,(1B1)	;NEED ROUNDING?
	  JRST	CALC		;NO, WON
	ADDI	T0,1		;YES, INCREMENT FRACTION OF HIGH WORD
	TLO	T0,400		;FIX POSSIBLE SPILL INTO EXPONENT
CALC:	CAMGE	T0,EPS		;IF F .LT. EPS,
	  JRST	RET1		;SIN OR COS = F (WITHIN SIGN)
	MOVE	T3,T0
	FMPR	T3,T3		;OBTAIN G
	MOVE	T2,T3		
	FMPR	T2,R5		;((((G*R5
	FADR	T2,R4		;+R4)
	FMPR	T2,T3		;*G
	FADR	T2,R3		;+R3)
	FMPR	T2,T3		;*G
	FADR	T2,R2		;+R2)
	FMPR	T2,T3		;*G
	FADR	T2,R1		;+R1)
	FMPR	T2,T3		;*G = R(G)
	FMPR	T2,T0		;*F
	FADR	T0,T2		;+F

RET1:	TRNE	T4,1		;IF SIGN FLAG SET,
	  MOVN	T0,T0		;NEGATE RESULT
RET:	POP	P,T4		;RESTORE ACCUMULATORS
	POP	P,T3
	POP	P,T2
	GOODBY(1)		;RETURN

OV:	MOVEI	T0,0		;RETURN ZERO
	$LCALL	ATZ
;LERR	(LIB,%,SIN or COS: ABS(arg) too large; result = zero)
	JRST	RET

PI:	EXP 202622077325,021026430215 	;PI
PI2:	EXP 201622077325,021026430215   ;PI/2
ONEDPI:	EXP 177505746033,162344202513	;1/PI
PI180:	EXP 173435750650,224516471053   ;PI/180
K180:	EXP 180.0,0			;180
R1:	601252525253		;-.166666666
R2:	172421042056		;.833333072E-2
R3:	613137720533		;-.198408328E-3
R4:	156561327224		;.275239711E-5
R5:	630145743634		;-.238683464E-7
EPS:	163552023642		;0.863167530E-4
 	PRGEND
TITLE	SINH	HYPERBOLIC SINE FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	SINH
EXTERN	SINH.
SINH=SINH.
PRGEND
TITLE	SINH.	HYPERBOLIC SINE FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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


;SINH IS CALCULATED AS FOLLOWS:

;  LET V BE APPROXIMATELY 2 SO THAT LN V AND ABS(X)+LN V
;  CAN BE EXACTLY REPRESENTED WHEN X IS EXACTLY REPRESENTABLE.
;  THEN, LETTING W = ABS(X), AND NOTING THAT -SINH(-X)=SINH(X),
;  FOR
;	0 <= W < EPS, SINH = W*SIGN(X)
;	EPS <= W <= 1, SINH = (W*P4(W**2))*SIGN(X), ALGORITHM 1
;	1 < W <= 88.029678, SINH = SIGN(X)*(EXP(W)-EXP(-W))/2, ALGORITHM 2
;	88.029678 < W < 128 * LN(2)
;		SINH = SIGN(X)*((V/2)*EXP(W - LN V)), ALGORITHM 3
;	W >= 128 * LN(2), SINH = SIGN(X) * MACHINE INFINITY

;	LET Z = W**2.  THEN
;		P4(Z) = 1 + Z*(C1 + Z*(C2 + Z*(C3 + C4*Z)))
;		    WHERE THE C(I) ARE GIVEN BELOW

;THE RANGE OF DEFINITION FOR SINH IS ABS(X) <= 88.722,
;  AND ARGUMENTS OUT OF THIS RANGE WILL CAUSE AN ERROR MESSAGE
;  TO BE TYPED.  A RESULT OF SIGN(X)*MACHINE INFINITY WILL BE RETURNED

;REQUIRED (CALLED) ROUTINES:  EXP

;REGISTER T2 IS SAVED, USED, AND RESTORED

;THE ROUTINE IS CALLED IN THE FOLLOWING MANNER:

;  XMOVEI	L,ARG
;  PUSHJ	P,SINH

;THE ANSWER IS RETURNED IN AC T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(SINH,.)
	MOVM	T0,@(L)		;ABS(X) = W
	CAMGE	T0,EPS		;IF W < EPS,
	  JRST	EXPP1		;RETURN--SINH = X
	CAMLE	T0,EIGHT8	;IF W > 88.029678,
	  JRST	ALG3		;GO TO ALGORITHM 3
	CAMG 	T0,ONE		;IF W <= 1,
	  JRST	ALG1		;GO TO ALGORITHM 1
;				 OTHERWISE, ALGORITHM 2:
	MOVEM	T0,TEMP
	FUNCT	EXP.,<TEMP>	;EXP(ABS(X))
	CAML	T0,TWO14	;IF EXP(ABS(X)).GT.2**14
	JRST	HALVE		;NEGLECT EXP(-ABS(X))
	HRLZI	T1,576400	;-1.0 
	FDVR	T1,T0		;-EXP(-ABS(X))
	FADR	T0,T1		;EXP(ABS(X))-EXP(-ABS(X))
HALVE:	FSC	T0,-1	        ;/2.0
	JRST	EXPP1
;				 ALGORITHM 1
ALG1:	PUSH	P,T2
	MOVE	T1,T0		;W
	FMPR	T0,T0		;W**2
	MOVE	T2,T0
	FMPR	T0,C4		;C4*W**2
	FADR	T0,C3		;+C3
	FMPR	T0,T2		;*W**2
	FADR	T0,C2		;+C2
	FMPR	T0,T2		;*W**2
	FADR	T0,C1		;+C1
	FMPR	T0,T2		;*W**2
	FMPR	T0,T1		;*W
	FADR	T0,T1		;+W = W*P4(Z)
	POP	P,T2		;RESTORE T2
	JRST	EXPP1
;				 ALGORITHM 3
ALG3:	CAMGE	T0,XXMAX	;IF ARGUMENT IS NOT TOO LARGE
	  JRST	EXPP		;CALCULATE AT EXPP
OVFL:	$LCALL	ROV
;LERR	(LIB,%,SINH: result overflow)
	HRLOI	T0,377777	;SINH = +MACHINE INFINITY.
	JRST	EXPP1		
EXPP:	FSBR	T0,LN2VE	;W-LN(V)
	MOVEM	T0,TEMP
	FUNCT	EXP.,<TEMP>	;EXP(W - LN V)
	MOVE	T1,T0		;SAVE A COPY OF EXP(2 - LN2VE)
	FMPR	T0,CON1		;MULTIPLY BY (LN2VE - LN2)
	FADR	T0,T1		;SUM = (1/2)*EXP(2)
	  JFCL	OVFL		;(ROUNDING ERRORS MIGHT CAUSE OVERFLOW)
EXPP1:	SKIPGE	@(L)		;CONSIDER SIGN OF SINH
	  MOVN	T0,T0	
	GOODBY	(1)		;SINH RETURN

LN2VE:	200542714000		;LN(V)=.693161011
EIGHT8:	207540074620		;88.029678
XXMAX:	207542710300
C1:	176525252524 		;1.666666643E-1
C2:	172421042352		;8.333352593E-3
C3:	164637771324		;1.983581245E-4
C4:	156572227373		;2.818523951E-6
CON1:	160720040562		;LN2VE - LN(2)
ONE:	1.0E0
EPS:	164400000000
TWO14:	217400000000

	SEGMENT	DATA

TEMP:	0			;TEMPORARY STORAGE FOR EXP ARG
	PRGEND
TITLE	SQRT	SQUARE ROOT FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	SQRT
EXTERN	SQRT.
SQRT=SQRT.
PRGEND
TITLE	SQRT.	SQUARE ROOT FUNCTION
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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

;SQRT(X) IS CALCULATED AS FOLLOWS

;  THE SQUARE ROOT OF ABS(X) IS CALCULATED; THERE IS AN ERROR MESSAGE NOTING
;  WHEN X IS NEGATIVE.  THE INITIAL GUESS IS OPTIMUM FOR NUMBERS 
;  BETWEEN 0.5 AND 1.0.  THIS GUESS IS FOLLOWED BY TWO ITERATIONS
;  OF NEWTON'S METHOD.

;THE RANGE OF DEFINITION FOR SQRT IS THE NON-NEGATIVE REPRESENTABLE
;  REAL NUMBERS

;REQUIRED (CALLED) ROUTINES:  NONE

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,SQRT

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(SQRT,.)	;ENTRY TO SQRT ROUTINE
	SKIPG	T1,@(L)		;OBTAIN X AND CHECK FOR NEGATIVITY
	  JRST	SQRTLE		;NO, HANDLE NON-POSITIVE ARGUMENT

SQRTP:	MOVE	T0,T1		;SAVE NUMBER
	LSH	T1,-1		;DIVIDE EXP BY 2
	TLZE	T1,400		;WAS EXPONENT ODD?
	  JRST	SQRT2		;YES

;HERE WHEN EXPONENT IS EVEN.  T1 CONTAINS AN UNNORMALIZED FLOATING
;  POINT NUMBER, THE FRACTION PART OF WHICH IS 1/2 THE FRACTION OF
;  THE ARGUMENT.  OUR INITIAL GUESS IS MADE BY A LINEAR APPROXIMATION
;  USING Y0 = SE (X + C), WHERE SE AND C ARE CONSTANTS USED FOR
;  EVEN EXPONENTS IN X.

	ADD	T1,[XWD 267,607000]	;COMPUTE LINEAR APPROXIMATION
	FMPRI	T1,301454		;RESCALE EXPONENT
	JRST	SQRT3

;HERE, WITH ODD EXPONENT, USE Y0 = SO * (X+C).

SQRT2:	ADD	T1,[XWD 267,607000]	;LINEAR APPROXIMATION
	FMPRI	T1,301650		;RESCALE EXPONENT

SQRT3:	FDV	T0,T1		;ORIGINAL / INITIAL GUESS
	FAD	T1,T0		;AVERAGE THEM
	FSC	T1,-1

	MOVM	T0,@(L)		;GET ORIGINAL NUMBER
	FDV	T0,T1		;SECOND ITERATION
	FADR	T0,T1
	FSC	T0,-1		;AVERAGE THIRD GUESS WITH SECOND

	GOODBY	(1)		;SQRT RETURN

SQRTLE:	JUMPE	T1,ZERO
	$LCALL	NAA
;LERR	(LIB,%,<SQRT: negative arg; result = SQRT(ABS(arg))>)
	MOVM	T1,T1		;GET MAGNITUDE
	JRST	SQRTP		;POSITIVE NOW

ZERO:	MOVEI	T0,0		;HERE ON ZERO ARG. RETURN ZERO
	GOODBY	(1)		;SQRT RETURN

	PRGEND
TITLE   COTAN	TANGENT AND COTANGENT FUNCTIONS
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC.      JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	COTAN
EXTERN	COTAN.
COTAN=COTAN.
PRGEND
TITLE	TAN	TANGENT AND COTANGENT FUNCTIONS  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC       JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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	TAN
EXTERN	TAN.
TAN=TAN.
PRGEND
TITLE	TAN.	TANGENT AND COTANGENT FUNCTIONS
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL  IMSL, INC       JUNE 1, 1979

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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.14159265.
;	THEN DEFINE N AND F SO THAT
;	 X=N*PI/4.0+F, 0.0 <= F <= PI/4.0, WITH PI=3.14159265.
;	THEN
;	 TAN(F) = F, IF F < 2**(-14)
;		= F*R(F**2), OTHERWISE
;	 WHERE
;	   R(F**2) = (P0+F**2*(P1+F**2*P2))/(Q0+F**2*(Q1+F**2))
;	 AND
;	   P0=62.604
;	   P1=-6.9716
;	   P2=6.7309
;	   Q0=P0
;	   Q1=-27.839
;	THE RESULT IS THEN RECIPROCATED, IF NECESSARY, AND GIVEN
;	THE APPROPRIATE SIGN.

; COEFFICIENTS ARE DERIVED FROM THOSE GIVEN IN CODY AND WAITE, "SOFTWARE
; MANUAL FOR ELEMENTARY FUNCTIONS", FOR MACHINES WITH 25-32 BIT PRECISION

;THE RANGE OF DEFINITION FOR TAN IS ABS(X) <= 2**26*(PI/2)
;  AND FOR COTAN(X), 2**(-126) * (1/2+2**(-27)) < ABS(X) <= 2**26 * (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,TAN
;	OR
;  PUSHJ	P,COTAN

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(COTAN,.)	;ENTRY TO COTAN ROUTINE
	PUSH	P,T2
	MOVM	T2,@(L)		;OBTAIN Y = ABS(X)
	MOVEI	T1,1		;1 = COTAN FLAG
	CAML	T2,EPS1		;IF Y IS NOT TOO SMALL FOR COTAN,
	  JRST 	ALG		;CONTINUE
	$LCALL	ROV
;LERR	(LIB,%,COTAN: result overflow)
	HRLOI	T0,377777	;ANSWER = + MACHINE INFINITY
	SKIPGE	T2,@(L)		;IF X IS NEGATIVE
	MOVN	T0,T0		;RESULT = -RESULT
	POP	P,T2		;RESTORE ACCUMULATOR
	GOODBY	(1)		;RETURN

	HELLO	(TAN,.)		;ENTRY TO TAN ROUTINE
	PUSH	P,T2
	MOVM    T2,@(L)		;OBTAIN Y = ABS(X)
	SETZ	T1,
ALG:	CAMG	T2,XMAX		;IF Y IS NOT TOO LARGE,
	  JRST	ALGC		;CONTINUE ALGORITHM
	$LCALL	ATZ
;LERR	(LIB,%,TAN or COTAN: ABS(arg) too large; result = zero)
	SETZ	T0,
	POP	P,T2		;RESTORE ACCUMULATOR
	GOODBY	(1)		;RETURN

ALGC:	PUSH	P,T3		;SAVE ACCUMULATORS
	PUSH	P,T4
	PUSH	P,T5
	SETZ	T3,
	CAMG	T2,PID4         ;SKIP ARGUMENT REDUCTION IF
	  JRST  COMP1           ;ARG .LT. PI/4
	DFMP	T2,FORDPI       ;T2+T3=(N+G)
	FIX	T4,T2           
	MOVE	T0,T4           ;SAVE N IN T0
	TRNE	T4,1            ;TEST N
	  JRST	ODD	        ;JUMP TO ODD IF N IS ODD
	FLTR	T4,T4
	SETZ	T5,
	DFSB	T2,T4	        ;T2,T3 = G
	DFMP	T2,PIDD4        ;T4 = PI/4*G = F
	TLNE	T3,(1B1)	;ROUND T2 AND
	  ADDI	T2,1		;GUARD AGAINST SPILL
	TLO	T2,(1B9)	;INTO EXPONENT
	JRST	COMP
ODD:	ADDI	T4,1	        ;T4 = N+1
	FLTR	T4,T4
	SETZ	T5,
	DFSB	T4,T2	        ;T2,T3 = 1-G
	DFMP	T4,PIDD4        ;T4 = PI/4(1-G) = PI/4-F
	TLNE	T5,(1B1)	;ROUND T2 AND
	  ADDI	T4,1		;GUARD AGAINST SPILL
	TLO	T4,(1B9)	;INTO EXPONENT
	MOVE	T2,T4	        ;T2 = PI/4-F
COMP:	MOVE	T3,T0	        ;MOVE N TO T3
COMP1:	CAMGE	T2,EPS		;IF F IS LESS THAN EPS
	  JRST	LT              ;GO TO LT
	MOVE	T5,T2           ;OTHERWISE
	FMPR 	T5,T5           ;G=F*F
	MOVE	T4,T5		;RESULT = F*(R-Q)/Q + F
	FMPR	T4,R2		;FORM R
	FADR	T4,R1
	FMPR	T4,T5
	MOVE	T0,T5
	FADR	T0,Q1		;FORM Q
	FMPR	T5,T0
	FSBR	T4,T5		;FORM R-Q
	FADR	T5,R0		;COMPLETE Q BY ADDING R0
	FDVR	T4,T5		;(R-Q)/Q
	FMPR	T4,T2		;*F
	FADR	T2,T4		;+F
LT:	MOVE	T0,T2
   	ANDI	T3,3
	CAIGE 	T3,2            ;IS N .GE. 2 ?
	  JRST	NLE2            ;NO, GO TO NLE2
	MOVN	T0,T0           ;RESULT = -RESULT
	MOVN	T3,T3
	ADDI	T3,3            ;N =3-N
NLE2:	CAMN	T3,T1           ;IF N IS EQUAL TO IFLAG
	  JRST	SSGN            ;GO TO SSGN
	HRLZI	T3,201400       ;PUT 1 IN T3
	FDVR	T3,T0           ;RESULT = 1/RESULT
	MOVE	T0,T3
SSGN:	SKIPGE	T2,@(L)		;IF X IS NEGATIVE
	  MOVN	T0,T0		;RESULT = -RESULT
	POP	P,T5
	POP	P,T4
	POP	P,T3
	POP	P,T2
      	GOODBY	(1)		;RETURN

EPS1:   002400000001            ;2**(-126)*(1/2 + 2**(-27))
XMAX:   233622077325		;PI/2 * 2**26
PID4:	200622077325            ;PI/4
PIDD4:  DOUBLE 200622077325,021026430216  ;PI/4  (DOUBLE PRECISION)
FORDPI:	DOUBLE 201505746033,162344202512  ;4/PI  (DOUBLE PRECISION)
EPS:    163400000000            ;2**(-14)
R0:	206764652343            ;62.60411195
R1:	574101637670            ;-6.971684006
R2:     175423545327            ;6.730910259*(10**(-2))
Q1:	572102441002            ;-27.83972122
	PRGEND
TITLE	TANH	HYPERBOLIC TANGENT FUNCTION  
;		(SINGLE PRECISION FLOATING POINT)
SUBTTL	IMSL, INC.	JUNE 1, 1979


;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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 > 9.8479016, TANH = 1.0*SIGN(X)
;	IF F > LN(3)/2 AND F <= 9.8479016, TANH = RESULT 1 =
;		SIGN(X)*(1 - 2/(EXP(2*F) + 1)))
;	IF F < 2**(-15), 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*(A + B*G)/(C + G).  A, B, AND C APPEAR BELOW

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

;REQUIRED (CALLED) ROUTINES:  EXP

;REGISTER T2 WAS SAVED, USED, AND RESTORED

;THE ROUTINE HAS THE FOLLOWING CALLING SEQUENCE:

;  XMOVEI	L,ARG
;  PUSHJ	P,TANH

;THE ANSWER IS RETURNED IN ACCUMULATOR T0

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(TANH,.)	;ENTRY TO TANH ROUTINE
	MOVM	T0,@(L)		;ABS(X) = F
	CAMLE	T0,AMAX		;IF F > 9.84...,
	  JRST	CAREA		;GO TO COMPLETION AREA
	CAMGE	T0,A215		;IF F < 2**(-15),
	  JRST	COMP1		;GO TO COMPLETION AREA
	PUSH	P,T2		;SAVE A REGISTER.
	MOVE	T1,T0		;COPY F INTO T1
	MOVE	T2,T0		;  AND T2.
	CAMLE	T2,ALN3		;IF F > (LN 3)/2,
	  JRST	RES1		;GO TO OBTAIN RESULT 1
;				 OBTAIN RESULT 2
	FMPR	T1,T1		;R(G) CALCULATION.  F*F = G
	MOVE	T0,T1
	FMPR	T0,B		;*B
	FADR	T0,A		;+A
	FMPR	T0,T1		;*G = G*(A + B*G)
	FADR	T1,C		
	FDVR	T0,T1		;R(G) = G*(A + B*G) / (C + G)
	FMPR	T0,T2		;*F
	FADR	T0,T2		;F + F*R(G) IS IN T0
	JRST	COMPL		;GO TO COMPLETION AREA
;				 COMPUTE RESULT 1
RES1:	FSC	T1,1		;2*F
	MOVEM	T1,TEMP
	FUNCT	EXP.,<TEMP>	;EXP(2*F)
	FADRI	T0,201400	;+1.0
	MOVSI	T1,575400	;OBTAIN -2.0
	FDVR	T1,T0		;/EXP(2*F)+1)
	FADRI	T1,201400	;+1. RESULT 1 COMPLETE
	MOVE	T0,T1		;PUT IT IN T0

COMPL:	POP	P,T2		;RESTORE ACCUMULATORS
COMP1:	SKIPGE	@(L)		;CHANGE TANH SIGN IF NECESSARY
	  MOVN	T0,T0
	GOODBY	(1)		;TANH RETURN

CAREA:	MOVSI	T0,201400	;RESULT = 1.0 IS IN T0
	JRST	COMP1

AMAX:	204473104012		;9.8479016
A215:	162400000000		;2**(-15)
ALN3:	200431175237		;0.549306145
A:	577132164714		;-0.823772813
B:	607011671163		;-0.383101067E-2
C:	202474250317		;2.47131965

	SEGMENT	DATA

TEMP:	0			;TEMPORARY STORAGE FOR EXP ARG
	PRGEND
TITLE	REAL.	Integer to real conversion
SUBTTL	H. P. WEISS/HPW/CDM	17-Nov-81

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	REAL.
REAL.=FLOAT.##
PRGEND
TITLE	FLOAT	INTEGER TO REAL CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	FLOAT
FLOAT=FLOAT.##
PRGEND
TITLE	FLOAT.	INTEGER TO REAL CONVERSION
SUBTTL	D. TODD /DRT 15-FEB-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(FLOAT,.)	;ENTRY TO FLOAT ROUTINE
	MOVE	T0,@(L)		;GET THE ARGUMENT
	FLTR	T0,T0		;FLOAT IT
	POPJ	17,
	PRGEND
TITLE	IFIX	REAL TO INTEGER CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	IFIX
IFIX=IFIX.##
PRGEND
TITLE	INT	REAL TO INTEGER CONVERSION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	INT
INT=INT.##
PRGEND
TITLE	IFIX.	REAL TO INTEGER CONVERSION
SUBTTL	D. TODD /DRT/ED YOURDON/KK/TWE/DMN Feb-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(INT,.)		;ENTRY TO INT ROUTINE.
	JRST	IFIX1		;GO TO MAIN ROUTINE.

	HELLO	(IFIX,.)	;ENTRY TO IFIX ROUTINE
IFIX1:	MOVE	T0,@(L)		;GET THE ARGUMENT
	FIX	T0,T0		;FIX IT
	POPJ	17,
	PRGEND

TITLE	ABS	SINGLE PRECISION ABSOLUTE VALUE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	ABS
ABS=ABS.##
PRGEND
TITLE	IABS	SINGLE PRECISION ABSOLUTE VALUE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	IABS
IABS=IABS.##
PRGEND
TITLE	ABS.	SINGLE PRECISION ABSOLUTE VALUE FUNCTION
SUBTTL	D. TODD /DRT/ED YOURDON/KK/TWE	15-Feb-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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

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

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


	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(ABS,.)		;ENTRY TO ABS ROUTINE
	MOVM	T0,@(L)		;GET /ARG/.
	GOODBY	(1)		;RETURN

	HELLO	(IABS,.)	;ENTRY TO IABS ROUTINE.
	MOVM	T0,@(L)		;GET /ARG/.  IF OVERFLOW, SET ANSWER TO +INF
	JOV	[HRLOI T0,377777
		 GOODBY (1)]
	GOODBY	(1)		;RETURN

	PRGEND
TITLE	DIM	SINGLE PRECISION POSITIVE DIFFERENCE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	DIM
DIM=DIM.##
PRGEND
TITLE	DIM.	SINGLE PRECISION DIFFERENCE FUNCTION
SUBTTL	D. TODD /DRT/ED YOURDON/TWE	15-Feb-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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

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

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(DIM,.)		;ENTRY TO DIM ROUTINE
	MOVE	T0,@0(L)	;PICK UP FIRST ARGUMENT
	CAMG	T0,@1(L)	;IF A .GT. B, GO TO SUBTRACT.
	TDZA	T0,T0		;O'E, ZERO A AND GO TO EXIT.
	FSBR	T0,@1(L)	;CALC A - B.
	 JFCL	EXCEP		;CAN UNDERFLOW AND OVERFLOW
RET:	POPJ	P,		;RETURN

EXCEP:	JUMPE	T0,UNDER	;UNDERFLOW?
	$LCALL	ROV,RET		;NO, OVERFLOW
UNDER:	$LCALL	RUN,RET		;YES, UNDERFLOW

	PRGEND
TITLE	IDIM	INTEGER POSITIVE DIFFERENCE FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	IDIM
IDIM=IDIM.##
PRGEND
TITLE	IDIM.	INTEGER POSITIVE DIFFERENCE FUNCTION
SUBTTL	D. TODD /DRT/ED YOURDON/KK/TWE	15-Feb-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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


;IDIM RETURNS THE POSITIVE DIFFERENCE OF A AND B:
;IF (A-B) .LE. 0, THEN DIM(A,B)=0
;IF(A-B) .G. 0 , THEN DIM(A,B)=(A-B)

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(IDIM,.)	;ENTRY TO IDIM ROUTINE
	MOVE	T0, @(L)	;PICK UP FIRST ARGUMENT
	CAMG	T0,@1(L)	;IF A .LE. B,
	TDZA	T0,T0		;SET ANSWER TO 0
	SUB	T0,@1(L)	;IF A .GT. B,
	JFCL	OVER		;OVERFLOW CAN OCCUR
RET:	GOODBY	(2)		;RETURN

OVER:	$LCALL	ROV		;RESULT OVERFLOW
	HRLOI	T0,377777	;STORE +BIGGEST
	JRST	RET		;RETURN

	PRGEND
TITLE	SIGN	SINGLE PREICISION XFER OF SIGN FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	SIGN
SIGN=SIGN.##
PRGEND
TITLE	ISIGN	SINGLE PREICISION XFER OF SIGN FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	ISIGN
ISIGN=ISIGN.##
PRGEND
TITLE	SIGN.	SINGLE PRECISION AND INTEGER XFER OF SIGN FUNCTION
SUBTTL	D. TODD /DRT/ED YOURDON/KK/TWE	12-Feb-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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.32(323)

;SIGN(A,B) AND ISIGN(A,B)
;IF B .GE. 0, THEN ABSF(A) IS RETURNED IN ACCUMULATOR 0
;IF B .L. 0, THEN -ABSF(A) IS RETURNED IN ACCUMULATOR 0

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(SIGN,.)	;ENTRY TO SIGN
	MOVM	T0,@(L)		;GET MAGNITUDE OF FIRST ARGUMENT
	SKIPGE	@1(L)		;IS SECOND ARGUMENT POSITIVE?
	MOVN	T0,T0		;NO, NEGATE RESULT
	GOODBY	(2)		;RETURN

	HELLO	(ISIGN,.)	;ENTRY TO ISIGN
	MOVM	T0,@(L)		;GET MAGNITUDE OF FIRST ARGUMENT
	SKIPGE	@1(L)		;IS SECOND ARGUMENT POSITIVE?
	MOVN	T0,T0		;NO, NEGATE RESULT
	JFCL	OVER		;CAN OVERFLOW
RET:	GOODBY	(2)		;RETURN

OVER:	$LCALL	ROV		;RESULT OVERFLOW
	HRLOI	T0,377777	;RETURN +BIGGEST
	JRST	RET		;RETURN

	PRGEND
TITLE	AINT	FLOATING POINT TRUNCATION FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	AINT
AINT=AINT.##
PRGEND
TITLE	AINT.	FLOATING POINT TRUNCATION FUNCTION
SUBTTL	ED YOURDON /KK/TWE	15-FEB-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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.32(323)

;FLOATING POINT TRUNCATION FUNCTION.
;TRUNCATES FRACTIONAL PART OF FLOATING POINT NUMBER
;AND RETURNS ANSWER AS A FLOATING POINT NUMBER.
;THE ANSWER IS RETURNED IN AC 0.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(AINT,.)	;ENTRY TO AINT ROUTINE.
	MOVM	T0,@(L)		;GET ABS(ARG)
	CAML	T0,MOD1		;IS ABS(ARG) .LT. 2**26?
	  JRST	AINT1		;NO, NO FRACTION BITS, EXIT.
	FAD	T0,MOD1		;YES, REMOVE
	FSB	T0,MOD1		;THE FRACTION BITS.
AINT1:	SKIPGE	@(L)		;SET THE
	  MOVN	T0,T0		;CORRECT SIGN AND
	GOODBY	(1)		;RETURN

MOD1:	233400000000		;2**26

	PRGEND
TITLE	MOD	INTEGER MOD FUNCTION
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	MOD
MOD=MOD.##
PRGEND
TITLE	MOD.	INTEGER MOD FUNCTION
SUBTTL	D. TODD /DRT/ED YOURDON/KK 15-Feb-1973

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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.32(323)

;INTEGER MOD FUNCTION
;MOD(A,B) = A-[A/B]*B, WHERE [A/B] IS THE GREATEST (IN
;MAGNITUDE) INTEGER IN A/B.  THAT IS, THE MOD FUNCTION 
;RETURNS THE REMAINDER OF THE QUOTIENT OF A AND B.  HENCE,
;9 MOD 2 IS 1, AND SO FORTH.

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(MOD,.)		;ENTRY TO MOD ROUTINE
	MOVE	T0,@(L)		;FIRST ARG TO AC 0.
	IDIV	T0,@1(L)	;DIVIDE, REMAINDER IN AC 1.
	MOVE	T0,T1		;PUT THE ANSWER IN AC 0.
	GOODBY	(2)		;RETURN

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	AMAX1
AMAX1=AMAX1.##
PRGEND
TITLE	MAX1	MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	MAX1
MAX1=MAX1.##
PRGEND
TITLE	AMAX0	MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	AMAX0
AMAX0=AMAX0.##
PRGEND
TITLE	MAX0	MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	MAX0
MAX0=MAX0.##
PRGEND
TITLE	MAX.	MAXIMUM OF A SERIES OF ARGUMENTS
SUBTTL	D. TODD /DRT/HPW	11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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.
;AMAX1, MAX1, AMAX0, AND MAX0 CALCULATE THE MAXIMUM OF A
;STRING OF ARGUMENTS

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(AMAX1,.)	;ENTRY TO AMAX1 ROUTINE.
	PUSHJ	P,MAX.		;FIND THE MAXIMUM AND
	GOODBY	(2)		;AMAX1 RETURN

	HELLO	(MAX0,.)	;ENTRY TO MAX0 ROUTINE.
	PUSHJ	P,MAX.		;FIND THE MAXIMUM AND
	GOODBY	(2)		;MAX0 RETURN

	HELLO	(MAX1,.)	;ENTRY TO MAX1 ROUTINE.
	PUSHJ	P,MAX.		;FIND THE MAXIMUM AND
	FIX	T0,T0		;FIX IT
	GOODBY	(2)		;MAX1 RETURN

	HELLO	(AMAX0,.)	;ENTRY TO AMAX0 ROUTINE.
	PUSHJ	P,MAX.		;FIND THE MAXIMUM AND
	FLTR	T0,T0		;FLOAT IT
	GOODBY	(2)		;AMAX0 RETURN

MAX.:	PUSH	P,T3		;Save an ac
	HLLZ	T3,-1(L)	;Get the F10 arg count
	JRST	MAX.2		;DON'T COMPARE FIRST TIME
MAX.1:	CAMGE	T0,@(L)		;IS CURRENT ARG .GT. NEXT ARG?
MAX.2:	MOVE	T0,@(L)		;NO, PUT ARG IN A.
	ADDI	L,1		;Bump arg pointer
	AOBJN	T3,MAX.1	;Continue thru the arg list
	POP	P,T3		;Restore ac
	POPJ	P,		;RETURN
	PRGEND
TITLE	AMIN1	MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	AMIN1
AMIN1=AMIN1.##
PRGEND
TITLE	MIN1	MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	MIN1
MIN1=MIN1.##
PRGEND
TITLE	AMIN0	MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	AMIN0
AMIN0=AMIN0.##
PRGEND
TITLE	MIN0	MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	H. P. WEISS/HPW		11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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	MIN0
MIN0=MIN0.##
PRGEND
TITLE	MIN.	MINIMUM OF A SERIES OF ARGUMENTS
SUBTTL	D. TODD /DRT/HPW	11-DEC-73

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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.
;AMIN1, MIN1, AMIN0, AND MIN0 CALCULATE THE MINIMUM OF A
;STRING OF ARGUMENTS

	SEARCH	MTHPRM
	SEGMENT	CODE
	NOSYM
	SALL

	HELLO	(AMIN1,.)	;ENTRY TO AMIN1 ROUTINE.
	PUSHJ	P,MIN.		;FIND THE MINIMUM AND
	GOODBY	(2)		;AMIN1 RETURN

	HELLO	(MIN0,.)	;ENTRY TO MIN0 ROUTINE.
	PUSHJ	P,MIN.		;FIND THE MINIMUM AND
	GOODBY	(2)		;MIN0 RETURN

	HELLO	(MIN1,.)	;ENTRY TO MIN1 ROUTINE.
	PUSHJ	P,MIN.		;FIND THE MINIMUM AND
	FIX	T0,T0		;FIX IT
	GOODBY	(2)		;MIN1 RETURN

	HELLO	(AMIN0,.)	;ENTRY TO AMIN0 ROUTINE.
	PUSHJ	P,MIN.		;FIND THE MINIMUM AND
	FLTR	T0,T0		;FLOAT IT
	GOODBY	(2)		;AMIN0 RETURN

MIN.:	PUSH	P,T3		;Save an ac
	HLLZ	T3,-1(L)	;Get the F10 arg count
	JRST	MIN.2		;DON'T COMPARE FIRST TIME
MIN.1:	CAMLE	T0,@(L)		;IS CURRENT ARG .LT. NEXT ARG?
MIN.2:	MOVE	T0,@(L)		;NO, PUT ARG IN A.
	ADDI	L,1		;Bump arg pointer
	AOBJN	T3,MIN.1	;Continue thru the arg list
	POP	P,T3		;Restore ac
	POPJ	P,		;EXIT.

	PRGEND
	TITLE ANINT	Floating point 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, 1983

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

	TITLE	ANINT.	Floating point nearest whole number.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

;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 single precision
; nearest whole number of the single precision number passed.
; (Both the number passed and returned are floating point.)

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

	SEARCH	MTHPRM
	SEGMENT	CODE
	SALL

	HELLO	(ANINT,.)	;Enter routine

	MOVM	T0,@0(L)	;Get argument to truncate.
	JUMPN	T0,NZERO	;Number passed = 0?
	GOODBYE			;Yes - Return

; Now shift out the fractional part of the number.

NZERO:	TLNN	T0,200000	;Is num .LT. 1/2?
	 JRST	ZERO		;Yes, go return zero.

	FAD	T0,[200400,,0]	;Add .5 before truncation.

	HLRZ	T1,T0		;Get new exponent.
	LSH	T1,-^D9		;Put rightmost.
	CAILE	T1,233		;If no fractional part,
	 JRST	DONE		; return the number passed.
				;ANINT(232777,,777777) must fall through.

	LSH	T0,-233(T1)	;Shift into integer position.
	MOVN	T1,T1		;Negate exponent.
	LSH	T0,233(T1)	;Shift back to where found.

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

	GOODBYE			;Return

ZERO:	TDZA	T0,T0		;Set result to zero and skip
DONE:	 MOVE	T0,@0(L)	;No fractional part, return
	GOODBYE			;Return

	PRGEND
	TITLE NINT	Integer nearest whole numer of floating

;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, 1983

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

	TITLE	NINT.	Integer nearest whole number of floating.
	SUBTTL	C. McCutcheon - 6/25/81

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1981, 1983

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

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

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

; Fortran Instrinsic Function to return the nearest integer
; to the single precision floating point 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	(NINT,.)	;Enter routine

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

NZERO:	FAD	T0,[200400,,0]	;Add .5 before truncation.

; Now shift out the fractional part of the number.

	HLRZ	T1,T0		;Get exponent.
	LSH	T1,-^D9		;Put rightmost.
	CAILE	T1,243		;If not representable as integer,
	 JRST	DONE		; return largest integer.

	FIX	T0,T0		;Make it into an integer

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

	GOODBYE			;Return

; Number not representable as integer, return largest integer.

DONE:	$LCALL	ROV
;LERR	(LIB,%,<NINT: result overflow>)
	HRLOI	T0,377777	;Largest positive number.
	SKIPGE	@0(L)		;If original number .LT. zero,
	 HRLZI	T0,400000	;Largest negative number.
	GOODBYE			;Return

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1973, 1983

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

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1973, 1983

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

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

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

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

	SEARCH MTHPRM
	EXTERN	SNG.0
	NOSYM
	SEGMENT	CODE
	SALL

	HELLO	(SNGL,.)	;[235] ENTRY TO SNGL ROUTINE
	DMOVE	T0,@(L)		;GET THE ARGUMENT IN AC 0
	PJRST	SNG.0		;USE AC0 SINGLE ROUTINE
	PRGEND
TITLE	SNG.0

;COPYRIGHT (C) DIGITAL EQUIPMENT CORPORATION 1979, 1983

;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
	SNG	0
	PRGEND


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


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


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


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


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


TITLE	SNG.14
	SEARCH	MTHPRM
	NOSYM
	SEGMENT	CODE
	SNG	14
	PRGEND			;[3223]
	TITLE	IFX.0		;[3223] Replace routines
	SEARCH	MTHPRM
	ENTRY	IFX.0
	SEGMENT	CODE
	NOSYM
IFX.0:	FIX	0,0
	POPJ	17,
	PRGEND
	TITLE	IFX.1
	SEARCH	MTHPRM
	ENTRY	IFX.1
	SEGMENT	CODE
	NOSYM
IFX.1:	FIX	1,1
	POPJ	17,
	PRGEND
	TITLE	IFX.2
	SEARCH	MTHPRM
	ENTRY	IFX.2
	SEGMENT	CODE
	NOSYM
IFX.2:	FIX	2,2
	POPJ	17,
	PRGEND
	TITLE	IFX.3
	SEARCH	MTHPRM
	ENTRY	IFX.3
	SEGMENT	CODE
	NOSYM
IFX.3:	FIX	3,3
	POPJ	17,
	PRGEND
	TITLE	IFX.4
	SEARCH	MTHPRM
	ENTRY	IFX.4
	SEGMENT	CODE
	NOSYM
IFX.4:	FIX	4,4
	POPJ	17,
	PRGEND
	TITLE	IFX.5
	SEARCH	MTHPRM
	ENTRY	IFX.5
	SEGMENT	CODE
	NOSYM
IFX.5:	FIX	5,5
	POPJ	17,
	PRGEND
	TITLE	IFX.6
	SEARCH	MTHPRM
	ENTRY	IFX.6
	SEGMENT	CODE
	NOSYM
IFX.6:	FIX	6,6
	POPJ	17,
	PRGEND
	TITLE	IFX.7
	SEARCH	MTHPRM
	ENTRY	IFX.7
	SEGMENT	CODE
	NOSYM
IFX.7:	FIX	7,7
	POPJ	17,
	PRGEND
	TITLE	IFX.10
	SEARCH	MTHPRM
	ENTRY	IFX.10
	SEGMENT	CODE
	NOSYM
IFX.10:	FIX	10,10
	POPJ	17,
	PRGEND
	TITLE	IFX.11
	SEARCH	MTHPRM
	ENTRY	IFX.11
	SEGMENT	CODE
	NOSYM
IFX.11:	FIX	11,11
	POPJ	17,
	PRGEND
	TITLE	IFX.12
	SEARCH	MTHPRM
	ENTRY	IFX.12
	SEGMENT	CODE
	NOSYM
IFX.12:	FIX	12,12
	POPJ	17,
	PRGEND
	TITLE	IFX.13
	SEARCH	MTHPRM
	ENTRY	IFX.13
	SEGMENT	CODE
	NOSYM
IFX.13:	FIX	13,13
	POPJ	17,
	PRGEND
	TITLE	IFX.14
	SEARCH	MTHPRM
	ENTRY	IFX.14
	SEGMENT	CODE
	NOSYM
IFX.14:	FIX	14,14
	POPJ	17,
	PRGEND
	TITLE	IFX.15
	SEARCH	MTHPRM
	ENTRY	IFX.15
	SEGMENT	CODE
	NOSYM
IFX.15:	FIX	15,15
	POPJ	17,
	PRGEND

	TITLE	MTHNRM	UNNORMALIZE UNDERFLOWED RESULT
	SEARCH	MTHPRM

	ENTRY	SPRUNM,DPRUNM

	SEGMENT	CODE

SPRUNM:	DMOVEM	T2,SAVE2	;SAVE T2 AND T3
	DMOVE	T2,@(16)	;GET RESULT STORED BY THE HARDWARE, HAS
				;CORRECT FRACTION AND EXPONENT TOO LARGE
				;BY 400
	HLRE	T1,T2		;GET EXPONENT AND SIGN AND SOME FRACTION BITS
	ASH	T1,-9		;GET RID OF FRACTION BITS
	TSCE	T1,T1		;GET ABS(EXPONENT), SKIP IF POSITIVE FRACTION
	  TLOA	T2,777000	;NEGATIVE FRACTION, SET EXPONENT TO ALL ONES
	TLZ	T2,777000	;POSITIVE FRACTION, SET EXPONENT TO ALL ZEROS
	CAME	T1,[377,,377]	;SUPPRESS ZERO-BIT SHIFT (-0 IS -256)
	  ASHC	T2,400001(T1)	;UNNORMALIZE FRACTION, KEEP 1 BIT FOR ROUNDING
	ADDI	T2,1		;ROUND HIGH WORD OF FRACTION
	ASH	T2,-1		;DISCARD ROUNDING BIT
	MOVE	T0,T2		;RETURN IN T0
	DMOVE	T2,SAVE2	;RESTORE T2 AND T3
	POPJ	P,

DPRUNM:	DMOVEM	T2,SAVE2	;SAVE T2 AND T3
	DMOVE	T2,@(16)	;GET RESULT STORED BY THE HARDWARE, HAS
				;CORRECT FRACTION AND EXPONENT TOO LARGE
				;BY 400
	HLRE	T1,T2		;GET EXPONENT AND SIGN AND SOME FRACTION BITS
	ASH	T1,-9		;GET RID OF FRACTION BITS
	TSCE	T1,T1		;GET ABS(EXPONENT), SKIP IF POSITIVE FRACTION
	  TLOA	T2,777000	;NEGATIVE FRACTION, SET EXPONENT TO ALL ONES
	TLZ	T2,777000	;POSITIVE FRACTION, SET EXPONENT TO ALL ZEROS
	CAME	T1,[377,,377]	;SUPPRESS ZERO-BIT SHIFT (-0 IS -256)
	  ASHC	T2,400001(T1)	;UNNORMALIZE FRACTION, KEEP 1 BIT FOR ROUNDING
	TLO	T3,(1B0)	;PREVENT INTEGER OVERFLOW WHEN WE ROUND
	ADDI	T3,1		;ROUND LOW WORD
	TLZN	T3,(1B0)	;DID FRACTION OVERFLOW INTO SIGN BIT?
	  ADDI	T2,1		;YES, PROPAGATE CARRY TO HIGH WORD
	ASHC	T2,-1		;DISCARD ROUNDING BIT
	DMOVE	T0,T2		;RETURN IN T0 AND T1
	DMOVE	T2,SAVE2	;RESTORE T2 AND T3
	POPJ	P,

	SEGMENT	DATA

SAVE2:	BLOCK	2		;FOR SAVING T2,T3

	END