;
;	File:		Power.a
;
;	Contains:	Routines to emulate exponential functions
;
;	Originally Written by:	Motorola Inc.
;	Adapted to Apple/MPW:	Jon Okada
;
;	Copyright:	© 1990,1991, 1993 by Apple Computer, Inc., all rights reserved.
;
;   This file is used in these builds:   Mac32
;
;	Change History (most recent first):
;
;		<SM2>	  2/3/93	CSS		Update from Horror:
;							<H2> 12/21/91 jmp (BG,Z4) (for JOkada) Corrected three constants under labels
;								      	  "EM1TINY:" and "EM12TINY" to obtain correct rounding behavior
;									  for FETOXM1 with small arguments.
; ———————————————————————————————————————————————————————————————————————————————————————
;	Pre-Horror ROM comments begin here.
; ———————————————————————————————————————————————————————————————————————————————————————
;		 <3>	 4/13/91	BG		Modified FTWOTOX emulation to not signal inexact on exact cases.
;		 <2>	 3/30/91	BG		Rolling in Jon Okada's latest changes.
;		 <1>	12/14/90	BG		First checked into TERROR/BBS.
;

;  power.a

;  Based upon Motorola files 'setox.sa' and 'stwotox.sa'.

;  setox

;  CHANGE LOG:
;  07 Jan 91	JPO	Deleted constants HUGE and TINY (not referenced).
;			  Moved constants and table EXPTBL to file
;			  'constants.a'.
;  28 Mar 91	JPO	Modified 'stwotox' to deliver exact results for
;			  integral input.  Streamlined some instruction
;			  streams throughout the file.
;  11 Dec 91	JPO	Corrected three constants under labels "EM1TINY:"
;			  and "EM12TINY" to obtain correct rounding behavior
;			  for FETOXM1 with small arguments.
;

*
*	setox.sa 3.1 12/10/90
*
*	The entry point setox computes the exponential of a value.
*	setoxd does the same except the input value is a denormalized
*	number.	setoxm1 computes exp(X)-1, and setoxm1d computes
*	exp(X)-1 for denormalized X.
*
*	INPUT
*	-----
*	Double-extended value in memory location pointed to by address
*	register a0.
*
*	OUTPUT
*	------
*	exp(X) or exp(X)-1 returned in floating-point register fp0.
*
*	ACCURACY and MONOTONICITY
*	-------------------------
*	The returned result is within 0.85 ulps in 64 significant bit, i.e.
*	within 0.5001 ulp to 53 bits if the result is subsequently rounded
*	to double precision. The result is provably monotonic in double
*	precision.
*
*	SPEED
*	-----
*	Two timings are measured, both in the copy-back mode. The
*	first one is measured when the function is invoked the first time
*	(so the instructions and data are not in cache), and the
*	second one is measured when the function is reinvoked at the same
*	input argument.
*
*	The program setox takes approximately 210/190 cycles for input
*	argument X whose magnitude is less than 16380 log2, which
*	is the usual situation.	For the less common arguments,
*	depending on their values, the program may run faster or slower --
*	but no worse than 10% slower even in the extreme cases.
*
*	The program setoxm1 takes approximately ???/??? cycles for input
*	argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes
*	approximately ???/??? cycles. For the less common arguments,
*	depending on their values, the program may run faster or slower --
*	but no worse than 10% slower even in the extreme cases.
*
*	ALGORITHM and IMPLEMENTATION NOTES
*	----------------------------------
*
*	setoxd
*	------
*	Step 1.	Set ans := 1.0
*
*	Step 2.	Return	ans := ans + sign(X)*2^(-126). Exit.
*	Notes:	This will always generate one exception -- inexact.
*
*
*	setox
*	-----
*
*	Step 1.	Filter out extreme cases of input argument.
*		1.1	If |X| >= 2^(-65), go to Step 1.3.
*		1.2	Go to Step 7.
*		1.3	If |X| < 16380 log(2), go to Step 2.
*		1.4	Go to Step 8.
*	Notes:	The usual case should take the branches 1.1 -> 1.3 -> 2.
*		 To avoid the use of floating-point comparisons, a
*		 compact representation of |X| is used. This format is a
*		 32-bit integer, the upper (more significant) 16 bits are
*		 the sign and biased exponent field of |X|; the lower 16
*		 bits are the 16 most significant fraction (including the
*		 explicit bit) bits of |X|. Consequently, the comparisons
*		 in Steps 1.1 and 1.3 can be performed by integer comparison.
*		 Note also that the constant 16380 log(2) used in Step 1.3
*		 is also in the compact form. Thus taking the branch
*		 to Step 2 guarantees |X| < 16380 log(2). There is no harm
*		 to have a small number of cases where |X| is less than,
*		 but close to, 16380 log(2) and the branch to Step 9 is
*		 taken.
*
*	Step 2.	Calculate N = round-to-nearest-int( X * 64/log2 ).
*		2.1	Set AdjFlag := 0 (indicates the branch 1.3 -> 2 was taken)
*		2.2	N := round-to-nearest-integer( X * 64/log2 ).
*		2.3	Calculate	J = N mod 64; so J = 0,1,2,..., or 63.
*		2.4	Calculate	M = (N - J)/64; so N = 64M + J.
*		2.5	Calculate the address of the stored value of 2^(J/64).
*		2.6	Create the value Scale = 2^M.
*	Notes:	The calculation in 2.2 is really performed by
*
*			Z := X * constant
*			N := round-to-nearest-integer(Z)
*
*		 where
*
*			constant := single-precision( 64/log 2 ).
*
*		 Using a single-precision constant avoids memory access.
*		 Another effect of using a single-precision "constant" is
*		 that the calculated value Z is
*
*			Z = X*(64/log2)*(1+eps), |eps| <= 2^(-24).
*
*		 This error has to be considered later in Steps 3 and 4.
*
*	Step 3.	Calculate X - N*log2/64.
*		3.1	R := X + N*L1, where L1 := single-precision(-log2/64).
*		3.2	R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
*	Notes:	a) The way L1 and L2 are chosen ensures L1+L2 approximate
*		 the value	-log2/64	to 88 bits of accuracy.
*		 b) N*L1 is exact because N is no longer than 22 bits and
*		 L1 is no longer than 24 bits.
*		 c) The calculation X+N*L1 is also exact due to cancellation.
*		 Thus, R is practically X+N(L1+L2) to full 64 bits.
*		 d) It is important to estimate how large can |R| be after
*		 Step 3.2.
*
*			N = rnd-to-int( X*64/log2 (1+eps) ), |eps|<=2^(-24)
*			X*64/log2 (1+eps)	=	N + f,	|f| <= 0.5
*			X*64/log2 - N	=	f - eps*X 64/log2
*			X - N*log2/64	=	f*log2/64 - eps*X
*
*
*		 Now |X| <= 16446 log2, thus
*
*			|X - N*log2/64| <= (0.5 + 16446/2^(18))*log2/64
*					<= 0.57 log2/64.
*		 This bound will be used in Step 4.
*
*	Step 4.	Approximate exp(R)-1 by a polynomial
*			p = R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
*	Notes:	a) In order to reduce memory access, the coefficients are
*		 made as "short" as possible: A1 (which is 1/2), A4 and A5
*		 are single precision; A2 and A3 are double precision.
*		 b) Even with the restrictions above,
*			|p - (exp(R)-1)| < 2^(-68.8) for all |R| <= 0.0062.
*		 Note that 0.0062 is slightly bigger than 0.57 log2/64.
*		 c) To fully utilize the pipeline, p is separated into
*		 two independent pieces of roughly equal complexities
*			p = [ R + R*S*(A2 + S*A4) ]	+
*				[ S*(A1 + S*(A3 + S*A5)) ]
*		 where S = R*R.
*
*	Step 5.	Compute 2^(J/64)*exp(R) = 2^(J/64)*(1+p) by
*				ans := T + ( T*p + t)
*		 where T and t are the stored values for 2^(J/64).
*	Notes:	2^(J/64) is stored as T and t where T+t approximates
*		 2^(J/64) to roughly 85 bits; T is in extended precision
*		 and t is in single precision. Note also that T is rounded
*		 to 62 bits so that the last two bits of T are zero. The
*		 reason for such a special form is that T-1, T-2, and T-8
*		 will all be exact --- a property that will give much
*		 more accurate computation of the function EXPM1.
*
*	Step 6.	Reconstruction of exp(X)
*			exp(X) = 2^M * 2^(J/64) * exp(R).
*		6.1	If AdjFlag = 0, go to 6.3
*		6.2	ans := ans * AdjScale
*		6.3	Restore the user FPCR
*		6.4	Return ans := ans * Scale. Exit.
*	Notes:	If AdjFlag = 0, we have X = Mlog2 + Jlog2/64 + R,
*		 |M| <= 16380, and Scale = 2^M. Moreover, exp(X) will
*		 neither overflow nor underflow. If AdjFlag = 1, that
*		 means that
*			X = (M1+M)log2 + Jlog2/64 + R, |M1+M| >= 16380.
*		 Hence, exp(X) may overflow or underflow or neither.
*		 When that is the case, AdjScale = 2^(M1) where M1 is
*		 approximately M. Thus 6.2 will never cause over/underflow.
*		 Possible exception in 6.4 is overflow or underflow.
*		 The inexact exception is not generated in 6.4. Although
*		 one can argue that the inexact flag should always be
*		 raised, to simulate that exception cost to much than the
*		 flag is worth in practical uses.
*
*	Step 7.	Return 1 + X.
*		7.1	ans := X
*		7.2	Restore user FPCR.
*		7.3	Return ans := 1 + ans. Exit
*	Notes:	For non-zero X, the inexact exception will always be
*		 raised by 7.3. That is the only exception raised by 7.3.
*		 Note also that we use the FMOVEM instruction to move X
*		 in Step 7.1 to avoid unnecessary trapping. (Although
*		 the FMOVEM may not seem relevant since X is normalized,
*		 the precaution will be useful in the library version of
*		 this code where the separate entry for denormalized inputs
*		 will be done away with.)
*
*	Step 8.	Handle exp(X) where |X| >= 16380log2.
*		8.1	If |X| > 16480 log2, go to Step 9.
*		(mimic 2.2 - 2.6)
*		8.2	N := round-to-integer( X * 64/log2 )
*		8.3	Calculate J = N mod 64, J = 0,1,...,63
*		8.4	K := (N-J)/64, M1 := truncate(K/2), M = K-M1, AdjFlag := 1.
*		8.5	Calculate the address of the stored value 2^(J/64).
*		8.6	Create the values Scale = 2^M, AdjScale = 2^M1.
*		8.7	Go to Step 3.
*	Notes:	Refer to notes for 2.2 - 2.6.
*
*	Step 9.	Handle exp(X), |X| > 16480 log2.
*		9.1	If X < 0, go to 9.3
*		9.2	ans := Huge, go to 9.4
*		9.3	ans := Tiny.
*		9.4	Restore user FPCR.
*		9.5	Return ans := ans * ans. Exit.
*	Notes:	Exp(X) will surely overflow or underflow, depending on
*		 X's sign. "Huge" and "Tiny" are respectively large/tiny
*		 extended-precision numbers whose square over/underflow
*		 with an inexact result. Thus, 9.5 always raises the
*		 inexact together with either overflow or underflow.
*
*
*	setoxm1d
*	--------
*
*	Step 1.	Set ans := 0
*
*	Step 2.	Return	ans := X + ans. Exit.
*	Notes:	This will return X with the appropriate rounding
*		 precision prescribed by the user FPCR.
*
*	setoxm1
*	-------
*
*	Step 1.	Check |X|
*		1.1	If |X| >= 1/4, go to Step 1.3.
*		1.2	Go to Step 7.
*		1.3	If |X| < 70 log(2), go to Step 2.
*		1.4	Go to Step 10.
*	Notes:	The usual case should take the branches 1.1 -> 1.3 -> 2.
*		 However, it is conceivable |X| can be small very often
*		 because EXPM1 is intended to evaluate exp(X)-1 accurately
*		 when |X| is small. For further details on the comparisons,
*		 see the notes on Step 1 of setox.
*
*	Step 2.	Calculate N = round-to-nearest-int( X * 64/log2 ).
*		2.1	N := round-to-nearest-integer( X * 64/log2 ).
*		2.2	Calculate	J = N mod 64; so J = 0,1,2,..., or 63.
*		2.3	Calculate	M = (N - J)/64; so N = 64M + J.
*		2.4	Calculate the address of the stored value of 2^(J/64).
*		2.5	Create the values Sc = 2^M and OnebySc := -2^(-M).
*	Notes:	See the notes on Step 2 of setox.
*
*	Step 3.	Calculate X - N*log2/64.
*		3.1	R := X + N*L1, where L1 := single-precision(-log2/64).
*		3.2	R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
*	Notes:	Applying the analysis of Step 3 of setox in this case
*		 shows that |R| <= 0.0055 (note that |X| <= 70 log2 in
*		 this case).
*
*	Step 4.	Approximate exp(R)-1 by a polynomial
*			p = R+R*R*(A1+R*(A2+R*(A3+R*(A4+R*(A5+R*A6)))))
*	Notes:	a) In order to reduce memory access, the coefficients are
*		 made as "short" as possible: A1 (which is 1/2), A5 and A6
*		 are single precision; A2, A3 and A4 are double precision.
*		 b) Even with the restriction above,
*			|p - (exp(R)-1)| <	|R| * 2^(-72.7)
*		 for all |R| <= 0.0055.
*		 c) To fully utilize the pipeline, p is separated into
*		 two independent pieces of roughly equal complexity
*			p = [ R*S*(A2 + S*(A4 + S*A6)) ]	+
*				[ R + S*(A1 + S*(A3 + S*A5)) ]
*		 where S = R*R.
*
*	Step 5.	Compute 2^(J/64)*p by
*				p := T*p
*		 where T and t are the stored values for 2^(J/64).
*	Notes:	2^(J/64) is stored as T and t where T+t approximates
*		 2^(J/64) to roughly 85 bits; T is in extended precision
*		 and t is in single precision. Note also that T is rounded
*		 to 62 bits so that the last two bits of T are zero. The
*		 reason for such a special form is that T-1, T-2, and T-8
*		 will all be exact --- a property that will be exploited
*		 in Step 6 below. The total relative error in p is no
*		 bigger than 2^(-67.7) compared to the final result.
*
*	Step 6.	Reconstruction of exp(X)-1
*			exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ).
*		6.1	If M <= 63, go to Step 6.3.
*		6.2	ans := T + (p + (t + OnebySc)). Go to 6.6
*		6.3	If M >= -3, go to 6.5.
*		6.4	ans := (T + (p + t)) + OnebySc. Go to 6.6
*		6.5	ans := (T + OnebySc) + (p + t).
*		6.6	Restore user FPCR.
*		6.7	Return ans := Sc * ans. Exit.
*	Notes:	The various arrangements of the expressions give accurate
*		 evaluations.
*
*	Step 7.	exp(X)-1 for |X| < 1/4.
*		7.1	If |X| >= 2^(-65), go to Step 9.
*		7.2	Go to Step 8.
*
*	Step 8.	Calculate exp(X)-1, |X| < 2^(-65).
*		8.1	If |X| < 2^(-16312), goto 8.3
*		8.2	Restore FPCR; return ans := X - 2^(-16382). Exit.
*		8.3	X := X * 2^(140).
*		8.4	Restore FPCR; ans := ans - 2^(-16382).
*		 Return ans := ans*2^(140). Exit
*	Notes:	The idea is to return "X - tiny" under the user
*		 precision and rounding modes. To avoid unnecessary
*		 inefficiency, we stay away from denormalized numbers the
*		 best we can. For |X| >= 2^(-16312), the straightforward
*		 8.2 generates the inexact exception as the case warrants.
*
*	Step 9.	Calculate exp(X)-1, |X| < 1/4, by a polynomial
*			p = X + X*X*(B1 + X*(B2 + ... + X*B12))
*	Notes:	a) In order to reduce memory access, the coefficients are
*		 made as "short" as possible: B1 (which is 1/2), B9 to B12
*		 are single precision; B3 to B8 are double precision; and
*		 B2 is double extended.
*		 b) Even with the restriction above,
*			|p - (exp(X)-1)| < |X| 2^(-70.6)
*		 for all |X| <= 0.251.
*		 Note that 0.251 is slightly bigger than 1/4.
*		 c) To fully preserve accuracy, the polynomial is computed
*		 as	X + ( S*B1 +	Q ) where S = X*X and
*			Q	=	X*S*(B2 + X*(B3 + ... + X*B12))
*		 d) To fully utilize the pipeline, Q is separated into
*		 two independent pieces of roughly equal complexity
*			Q = [ X*S*(B2 + S*(B4 + ... + S*B12)) ] +
*				[ S*S*(B3 + S*(B5 + ... + S*B11)) ]
*
*	Step 10.	Calculate exp(X)-1 for |X| >= 70 log 2.
*		10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all practical
*		 purposes. Therefore, go to Step 1 of setox.
*		10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical purposes.
*		 ans := -1
*		 Restore user FPCR
*		 Return ans := ans + 2^(-126). Exit.
*	Notes:	10.2 will always create an inexact and return -1 + tiny
*		 in the user rounding precision and mode.
*
*

*		Copyright (C) Motorola, Inc. 1990
*			All Rights Reserved
*
*	THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA 
*	The copyright notice above does not evidence any  
*	actual or intended publication of such source code.

* setox	IDNT	2,1 Motorola 040 Floating Point Software Package


ADJFLAG	equ L_SCR2
SCALE	equ FP_SCR1
ADJSCALE equ FP_SCR2
SC	equ FP_SCR3
ONEBYSC	equ FP_SCR4


setoxd:
*--entry point for EXP(X), X is denormalized
	MOVE.L	(a0),d0
	ANDI.L	#$80000000,d0
	ORI.L	#$00800000,d0		...sign(X)*2^(-126)
	MOVE.L	d0,-(sp)
	FMOVE.S	#"$3F800000",fp0
	fmove.l	d1,fpcr
	FADD.S	(sp)+,fp0
	bra		t_frcinx


setox:
*--entry point for EXP(X), here X is finite, non-zero, and not NaN's

*--Step 1.
	MOVE.L	(a0),d0		...load part of input X
	ANDI.L	#$7FFF0000,d0	...biased expo. of X
	CMPI.L	#$3FBE0000,d0	...2^(-65)
	BGE.B	EXPC1		...normal case
	BRA.W	EXPSM

EXPC1:
*--The case |X| >= 2^(-65)
	MOVE.W	4(a0),d0	...expo. and partial sig. of |X|
	CMPI.L	#$400CB167,d0	...16380 log2 trunc. 16 bits
	BLT.B	EXPMAIN		...normal case
	BRA.W	EXPBIG

EXPMAIN:
*--Step 2.
*--This is the normal branch:	2^(-65) <= |X| < 16380 log2.
	FMOVE.X	(a0),fp0		...load input from (a0)

	FMOVE.X	fp0,fp1
	FMUL.S	#"$42B8AA3B",fp0	...64/log2 * X
	fmovem.x fp2/fp3,-(a7)		...save fp2
;	MOVE.L	#0,ADJFLAG(a6)		; DELETED <3/28/91, JPO>			<T3>
	clr.l	ADJFLAG(a6)		; <3/28/91, JPO>				<T3>
	FMOVE.L	fp0,d0			...N = int( X * 64/log2 )
	LEA	EXPTBL,a1
	FMOVE.L	d0,fp0			...convert to floating-format

	MOVE.L	d0,L_SCR1(a6)		...save N temporarily
	ANDI.L	#$3F,d0			...D0 is J = N mod 64
	LSL.L	#4,d0
	ADDA.L	d0,a1			...address of 2^(J/64)
	MOVE.L	L_SCR1(a6),d0
	ASR.L	#6,d0			...D0 is M
	ADDI.W	#$3FFF,d0		...biased expo. of 2^(M)
	MOVE.W	L2,L_SCR1(a6)		...prefetch L2, no need in CB

EXPCONT1:
*--Step 3.
*--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
*--a0 points to 2^(J/64), D0 is biased expo. of 2^(M)
	FMOVE.X	fp0,fp2
	FMUL.S	#"$BC317218",fp0	...N * L1, L1 = lead(-log2/64)
	FMUL.X	L2,fp2			...N * L2, L1+L2 = -log2/64
	FADD.X	fp1,fp0	 		...X + N*L1
	FADD.X	fp2,fp0			...fp0 is R, reduced arg.
*	MOVE.W	#$3FA5,EXPA3		...load EXPA3 in cache

*--Step 4.
*--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
*-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
*--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
*--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))]

	FMOVE.X	fp0,fp1
	FMUL.X	fp1,fp1	 		...fp1 IS S = R*R

	FMOVE.S	#"$3AB60B70",fp2	...fp2 IS A5
*	MOVE.W	#0,2(a1)		...load 2^(J/64) in cache

	FMUL.X	fp1,fp2	 		...fp2 IS S*A5
	FMOVE.X	fp1,fp3
	FMUL.S	#"$3C088895",fp3	...fp3 IS S*A4

	FADD.D	EXPA3,fp2		...fp2 IS A3+S*A5
	FADD.D	EXPA2,fp3		...fp3 IS A2+S*A4

	FMUL.X	fp1,fp2	 		...fp2 IS S*(A3+S*A5)
	MOVE.W	d0,SCALE(a6)		...SCALE is 2^(M) in extended
	clr.w	SCALE+2(a6)
	move.l	#$80000000,SCALE+4(a6)
	clr.l	SCALE+8(a6)

	FMUL.X	fp1,fp3	 		...fp3 IS S*(A2+S*A4)

	FADD.S	#"$3F000000",fp2	...fp2 IS A1+S*(A3+S*A5)
	FMUL.X	fp0,fp3	 		...fp3 IS R*S*(A2+S*A4)

	FMUL.X	fp1,fp2	 		...fp2 IS S*(A1+S*(A3+S*A5))
	FADD.X	fp3,fp0	 		...fp0 IS R+R*S*(A2+S*A4),
*					...fp3 released

	FMOVE.X	(a1)+,fp1		...fp1 is lead. pt. of 2^(J/64)
	FADD.X	fp2,fp0	 		...fp0 is EXP(R) - 1
*					...fp2 released

*--Step 5
*--final reconstruction process
*--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) )

	FMUL.X	fp1,fp0	 		...2^(J/64)*(Exp(R)-1)
	fmovem.x (a7)+,fp2/fp3		...fp2 restored
	FADD.S	(a1),fp0		...accurate 2^(J/64)

	FADD.X	fp1,fp0	 		...2^(J/64) + 2^(J/64)*...
	MOVE.L	ADJFLAG(a6),d0

*--Step 6
	TST.L	D0
	BEQ.B	NORMAL
ADJUST:
	FMUL.X	ADJSCALE(a6),fp0
NORMAL:
	FMOVE.L	d1,FPCR	 		...restore user FPCR
	FMUL.X	SCALE(a6),fp0		...multiply 2^(M)
	bra	t_frcinx

EXPSM:
*--Step 7
	FMOVEM.X (a0),fp0		...in case X is denormalized
	FMOVE.L	d1,FPCR
	FADD.S	#"$3F800000",fp0	...1+X in user mode
	bra	t_frcinx

EXPBIG:
*--Step 8
	CMPI.L	#$400CB27C,d0		...16480 log2
	BGT.B	EXP2BIG
*--Steps 8.2 -- 8.6
	FMOVE.X	(a0),fp0		...load input from (a0)

	FMOVE.X	fp0,fp1
	FMUL.S	#"$42B8AA3B",fp0	...64/log2 * X
	fmovem.x fp2/fp3,-(a7)		...save fp2
	MOVE.L	#1,ADJFLAG(a6)
	FMOVE.L	fp0,d0			...N = int( X * 64/log2 )
	LEA	EXPTBL,a1
	FMOVE.L	d0,fp0			...convert to floating-format
	MOVE.L	d0,L_SCR1(a6)		...save N temporarily
	ANDI.L	#$3F,d0			 ...D0 is J = N mod 64
	LSL.L	#4,d0
	ADDA.L	d0,a1			...address of 2^(J/64)
	MOVE.L	L_SCR1(a6),d0
	ASR.L	#6,d0			...D0 is K
	MOVE.L	d0,L_SCR1(a6)		...save K temporarily
	ASR.L	#1,d0			...D0 is M1
	SUB.L	d0,L_SCR1(a6)			...a1 is M
	ADDI.W	#$3FFF,d0		...biased expo. of 2^(M1)
	MOVE.W	d0,ADJSCALE(a6)		...ADJSCALE := 2^(M1)
	clr.w	ADJSCALE+2(a6)
	move.l	#$80000000,ADJSCALE+4(a6)
	clr.l	ADJSCALE+8(a6)
	MOVE.L	L_SCR1(a6),d0		...D0 is M
	ADDI.W	#$3FFF,d0		...biased expo. of 2^(M)
	BRA.W	EXPCONT1		...go back to Step 3

EXP2BIG:
*--Step 9
	FMOVE.L	d1,FPCR
	MOVE.L	(a0),d0
	bclr.b	#sign_bit,(a0)		...setox always returns positive
;	CMPI.L	#0,d0			; DELETED <3/28/91, JPO>			<T3>
;	BLT	t_unfl			; DELETED <3/28/91, JPO>			<T3>
	tst.l	d0			; <3/28/91, JPO>				<T3>
	bmi	t_unfl			; <3/28/91, JPO>				<T3>
	BRA	t_ovfl


setoxm1d:
*--entry point for EXPM1(X), here X is denormalized
*--Step 0.
	bra	t_extdnrm



setoxm1:
*--entry point for EXPM1(X), here X is finite, non-zero, non-NaN

*--Step 1.
*--Step 1.1
	MOVE.L	(a0),d0	 	...load part of input X
	ANDI.L	#$7FFF0000,d0	...biased expo. of X
	CMPI.L	#$3FFD0000,d0	...1/4
	BGE.B	EM1CON1		...|X| >= 1/4
	BRA.W	EM1SM

EM1CON1:
*--Step 1.3
*--The case |X| >= 1/4
	MOVE.W	4(a0),d0	...expo. and partial sig. of |X|
	CMPI.L	#$4004C215,d0	...70log2 rounded up to 16 bits
	BLE.B	EM1MAIN		...1/4 <= |X| <= 70log2
	BRA.W	EM1BIG

EM1MAIN:
*--Step 2.
*--This is the case:	1/4 <= |X| <= 70 log2.
	FMOVE.X	(a0),fp0	...load input from (a0)

	FMOVE.X	fp0,fp1
	FMUL.S	#"$42B8AA3B",fp0	...64/log2 * X
	fmovem.x fp2/fp3,-(a7)		...save fp2
*	MOVE.W	#$3F81,EM1A4		...prefetch in CB mode
	FMOVE.L	fp0,d0			...N = int( X * 64/log2 )
	LEA	EXPTBL,a1
	FMOVE.L	d0,fp0			...convert to floating-format

	MOVE.L	d0,L_SCR1(a6)		...save N temporarily
	ANDI.L	#$3F,d0			...D0 is J = N mod 64
	LSL.L	#4,d0
	ADDA.L	d0,a1			...address of 2^(J/64)
	MOVE.L	L_SCR1(a6),d0
	ASR.L	#6,d0			...D0 is M
	MOVE.L	d0,L_SCR1(a6)		...save a copy of M
*	MOVE.W	#$3FDC,L2		...prefetch L2 in CB mode

*--Step 3.
*--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
*--a0 points to 2^(J/64), D0 and a1 both contain M
	FMOVE.X	fp0,fp2
	FMUL.S	#"$BC317218",fp0	...N * L1, L1 = lead(-log2/64)
	FMUL.X	L2,fp2			...N * L2, L1+L2 = -log2/64
	FADD.X	fp1,fp0			...X + N*L1
	FADD.X	fp2,fp0			...fp0 is R, reduced arg.
*	MOVE.W	#$3FC5,EM1A2		...load EM1A2 in cache
	ADDI.W	#$3FFF,d0		...D0 is biased expo. of 2^M

*--Step 4.
*--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
*-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6)))))
*--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
*--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))]

	FMOVE.X	fp0,fp1
	FMUL.X	fp1,fp1			...fp1 IS S = R*R

	FMOVE.S	#"$3950097B",fp2	...fp2 IS a6
*	MOVE.W	#0,2(a1)		...load 2^(J/64) in cache

	FMUL.X	fp1,fp2			...fp2 IS S*A6
	FMOVE.X	fp1,fp3
	FMUL.S	#"$3AB60B6A",fp3	...fp3 IS S*A5

	FADD.D	EM1A4,fp2		...fp2 IS A4+S*A6
	FADD.D	EM1A3,fp3		...fp3 IS A3+S*A5
	MOVE.W	d0,SC(a6)		...SC is 2^(M) in extended
	clr.w	SC+2(a6)
	move.l	#$80000000,SC+4(a6)
	clr.l	SC+8(a6)

	FMUL.X	fp1,fp2			...fp2 IS S*(A4+S*A6)
	MOVE.L	L_SCR1(a6),d0		...D0 is M
	NEG.W	D0			...D0 is -M
	FMUL.X	fp1,fp3			...fp3 IS S*(A3+S*A5)
	ADDI.W	#$3FFF,d0		...biased expo. of 2^(-M)
	FADD.D	EM1A2,fp2		...fp2 IS A2+S*(A4+S*A6)
	FADD.S	#"$3F000000",fp3	...fp3 IS A1+S*(A3+S*A5)

	FMUL.X	fp1,fp2			...fp2 IS S*(A2+S*(A4+S*A6))
	ORI.W	#$8000,d0		...signed/expo. of -2^(-M)
	MOVE.W	d0,ONEBYSC(a6)		...OnebySc is -2^(-M)
	clr.w	ONEBYSC+2(a6)
	move.l	#$80000000,ONEBYSC+4(a6)
	clr.l	ONEBYSC+8(a6)
	FMUL.X	fp3,fp1			...fp1 IS S*(A1+S*(A3+S*A5))
*					...fp3 released

	FMUL.X	fp0,fp2			...fp2 IS R*S*(A2+S*(A4+S*A6))
	FADD.X	fp1,fp0			...fp0 IS R+S*(A1+S*(A3+S*A5))
*					...fp1 released

	FADD.X	fp2,fp0			...fp0 IS EXP(R)-1
*					...fp2 released
	fmovem.x (a7)+,fp2/fp3		...fp2 restored

*--Step 5
*--Compute 2^(J/64)*p

	FMUL.X	(a1),fp0		...2^(J/64)*(Exp(R)-1)

*--Step 6
*--Step 6.1
	MOVE.L	L_SCR1(a6),d0		...retrieve M
	CMPI.L	#63,d0
	BLE.B	MLE63
*--Step 6.2	M >= 64
	FMOVE.S	12(a1),fp1		...fp1 is t
	FADD.X	ONEBYSC(a6),fp1		...fp1 is t+OnebySc
	FADD.X	fp1,fp0			...p+(t+OnebySc), fp1 released
	FADD.X	(a1),fp0		...T+(p+(t+OnebySc))
	BRA.B	EM1SCALE
MLE63:
*--Step 6.3	M <= 63
	CMPI.L	#-3,d0
	BGE.B	MGEN3
MLTN3:
*--Step 6.4	M <= -4
	FADD.S	12(a1),fp0		...p+t
	FADD.X	(a1),fp0		...T+(p+t)
	FADD.X	ONEBYSC(a6),fp0		...OnebySc + (T+(p+t))
	BRA.B	EM1SCALE
MGEN3:
*--Step 6.5	-3 <= M <= 63
	FMOVE.X	(a1)+,fp1		...fp1 is T
	FADD.S	(a1),fp0		...fp0 is p+t
	FADD.X	ONEBYSC(a6),fp1		...fp1 is T+OnebySc
	FADD.X	fp1,fp0			...(T+OnebySc)+(p+t)

EM1SCALE:
*--Step 6.6
	FMOVE.L	d1,FPCR
	FMUL.X	SC(a6),fp0

	bra	t_frcinx

EM1SM:
*--Step 7	|X| < 1/4.
	CMPI.L	#$3FBE0000,d0		...2^(-65)
	BGE.B	EM1POLY

EM1TINY:
*--Step 8	|X| < 2^(-65)
;	CMPI.L	#$00330000,d0		...2^(-16312) - DELETED <12/11/91, JPO>		<Z4><H2>

	cmpi.l	#$00470000,d0		; compare |X| with 2^(-16312)  <12/11/91, JPO>	<Z4><H2>

	BLT.B	EM12TINY
*--Step 8.2
;	MOVE.L	#$80010000,SC(a6)	...SC is -2^(-16382) - DELETED <12/11/91, JPO>	<Z4><H2>

	move.l	#$00010000,SC(a6)	; SC is +2^(-16382)  <12/11/91, JPO>		<Z4><H2>

	move.l	#$80000000,SC+4(a6)
	clr.l	SC+8(a6)
	FMOVE.X	(a0),fp0
	FMOVE.L	d1,FPCR
	FADD.X	SC(a6),fp0

	bra	t_frcinx

EM12TINY:
*--Step 8.3
	FMOVE.X	(a0),fp0
	FMUL.D	TWO140,fp0
;	MOVE.L	#$80010000,SC(a6)	; DELETED <12/11/91, JPO>			<Z4><H2>

	move.l	#$00010000,SC(a6)	; SC is +2^(-16382)  <12/11/91, JPO>		<Z4><H2>

	move.l	#$80000000,SC+4(a6)
	clr.l	SC+8(a6)
	FADD.X	SC(a6),fp0
	FMOVE.L	d1,FPCR
	FMUL.D	TWON140,fp0

	bra	t_frcinx

EM1POLY:
*--Step 9	exp(X)-1 by a simple polynomial
	FMOVE.X	(a0),fp0		...fp0 is X
	FMUL.X	fp0,fp0			...fp0 is S := X*X
	fmovem.x fp2/fp3,-(a7)		...save fp2
	FMOVE.S	#"$2F30CAA8",fp1	...fp1 is B12
	FMUL.X	fp0,fp1			...fp1 is S*B12
	FMOVE.S	#"$310F8290",fp2	...fp2 is B11
	FADD.S	#"$32D73220",fp1	...fp1 is B10+S*B12

	FMUL.X	fp0,fp2			...fp2 is S*B11
	FMUL.X	fp0,fp1			...fp1 is S*(B10 + ...

	FADD.S	#"$3493F281",fp2	...fp2 is B9+S*...
	FADD.D	EM1B8,fp1		...fp1 is B8+S*...

	FMUL.X	fp0,fp2			...fp2 is S*(B9+...
	FMUL.X	fp0,fp1			...fp1 is S*(B8+...

	FADD.D	EM1B7,fp2		...fp2 is B7+S*...
	FADD.D	EM1B6,fp1		...fp1 is B6+S*...

	FMUL.X	fp0,fp2			...fp2 is S*(B7+...
	FMUL.X	fp0,fp1			...fp1 is S*(B6+...

	FADD.D	EM1B5,fp2		...fp2 is B5+S*...
	FADD.D	EM1B4,fp1		...fp1 is B4+S*...

	FMUL.X	fp0,fp2			...fp2 is S*(B5+...
	FMUL.X	fp0,fp1			...fp1 is S*(B4+...

	FADD.D	EM1B3,fp2		...fp2 is B3+S*...
	FADD.X	EM1B2,fp1		...fp1 is B2+S*...

	FMUL.X	fp0,fp2			...fp2 is S*(B3+...
	FMUL.X	fp0,fp1			...fp1 is S*(B2+...

	FMUL.X	fp0,fp2			...fp2 is S*S*(B3+...)
	FMUL.X	(a0),fp1		...fp1 is X*S*(B2...
	
	FMUL.S	#"$3F000000",fp0	...fp0 is S*B1
	FADD.X	fp2,fp1			...fp1 is Q
*					...fp2 released

	fmovem.x (a7)+,fp2/fp3		...fp2 restored

	FADD.X	fp1,fp0			...fp0 is S*B1+Q
*					...fp1 released

	FMOVE.L	d1,FPCR
	FADD.X	(a0),fp0

	bra	t_frcinx

EM1BIG:
*--Step 10	|X| > 70 log2
	MOVE.L	(a0),d0
;	CMPI.L	#0,d0			; DELETED <3/28/91, JPO>			<T3>
	tst.l	d0			; <3/28/91, JPO>				<T3>
	BGT.W	EXPC1
*--Step 10.2
	FMOVE.S	#"$BF800000",fp0	...fp0 is -1
	FMOVE.L	d1,FPCR
	FADD.S	#"$00800000",fp0	...-1 + 2^(-126)

	bra	t_frcinx




;  stwotox

;  CHANGE LOG:
;  07 Jan 91	JPO	Deleted constants BOUNDS1, BOUNDS2, HUGE, and TINY
;			  (not referenced).  Changed constant labels EXPA1-EXPA5
;			  to EXPT1-EXPT5 and table label EXPTBL to EXP2TBL.
;			  Moved constants and table EXP2TBL to file 'constants.a'.
;			  Changed local variable names X and N to XPWR and NPWR,
;			  respectively.  Deleted local variable names XDCARE and
;			  XFRAC (not referenced).  Renamed label "EXPBIG" to
;			  "EXPBIG2".  Deleted unreferenced label "EXPSM".
;

*
*	stwotox.sa 3.1 12/10/90
*
*	stwotox  --- 2**X
*	stwotoxd --- 2**X for denormalized X
*	stentox  --- 10**X
*	stentoxd --- 10**X for denormalized X
*
*	Input: Double-extended number X in location pointed to
*		by address register a0.
*
*	Output: The function values are returned in Fp0.
*
*	Accuracy and Monotonicity: The returned result is within 2 ulps in
*		64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
*		result is subsequently rounded to double precision. The
*		result is provably monotonic in double precision.
*
*	Speed: The program stwotox takes approximately 190 cycles and the
*		program stentox takes approximately 200 cycles.
*
*	Algorithm:
*
*	twotox
*	1. If |X| > 16480, go to ExpBig2.
*
*	2. If |X| < 2**(-70), go to Exp2Sm.
*
*	3. Decompose X as X = N/64 + r where |r| <= 1/128. Furthermore
*		decompose N as
*		 N = 64(M + M') + j,  j = 0,1,2,...,63.
*
*	4. Overwrite r := r * log2. Then
*		2**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
*		Go to expr to compute that expression.
*
*	tentox
*	1. If |X| > 16480*log_10(2) (base 10 log of 2), go to ExpBig2.
*
*	2. If |X| < 2**(-70), go to Exp2Sm.
*
*	3. Set y := X*log_2(10)*64 (base 2 log of 10). Set
*		N := round-to-int(y). Decompose N as
*		 N = 64(M + M') + j,  j = 0,1,2,...,63.
*
*	4. Define r as
*		r := ((X - N*L1)-N*L2) * L10
*		where L1, L2 are the leading and trailing parts of log_10(2)/64
*		and L10 is the natural log of 10. Then
*		10**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
*		Go to expr to compute that expression.
*
*	expr
*	1. Fetch 2**(j/64) from table as Fact1 and Fact2.
*
*	2. Overwrite Fact1 and Fact2 by
*		Fact1 := 2**(M) * Fact1
*		Fact2 := 2**(M) * Fact2
*		Thus Fact1 + Fact2 = 2**(M) * 2**(j/64).
*
*	3. Calculate P where 1 + P approximates exp(r):
*		P = r + r*r*(A1+r*(A2+...+r*A5)).
*
*	4. Let AdjFact := 2**(M'). Return
*		AdjFact * ( Fact1 + ((Fact1*P) + Fact2) ).
*		Exit.
*
*	ExpBig2
*	1. Generate overflow by Huge * Huge if X > 0; otherwise, generate
*		underflow by Tiny * Tiny.
*
*	Exp2Sm
*	1. Return 1 + X.
*

*		Copyright (C) Motorola, Inc. 1990
*			All Rights Reserved
*
*	THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA 
*	The copyright notice above does not evidence any  
*	actual or intended publication of such source code.

* STWOTOX	IDNT	2,1 Motorola 040 Floating Point Software Package



;N	equ	L_SCR1		; renamed <1/7/91, JPO>

NPWR	equ	L_SCR1		; <1/7/91, JPO>

;X	equ	FP_SCR1		; renamed <1/7/91, JPO>

XPWR	equ	FP_SCR1		; <1/7/91, JPO>

;XDCARE	equ	X+2		; removed <1/7/91, JPO>
;XFRAC	equ	X+4		; removed <1/7/91, JPO>

ADJFACT	equ	FP_SCR2

FACT1	equ	FP_SCR3
FACT1HI	equ	FACT1+4
FACT1LOW equ	FACT1+8

FACT2	equ	FP_SCR4
FACT2HI	equ	FACT2+4
FACT2LOW equ	FACT2+8


stwotoxd:
*--ENTRY POINT FOR 2**(X) FOR DENORMALIZED ARGUMENT

	fmove.l	d1,fpcr			...set user's rounding mode/precision
	Fmove.S	#"$3F800000",FP0	...RETURN 1 + X
	move.l	(a0),d0
	or.l	#$00800001,d0
	fadd.s	d0,fp0
	bra	t_frcinx


stwotox:
*--ENTRY POINT FOR 2**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
	FMOVEM.X (a0),FP0	...LOAD INPUT, do not set cc's

	MOVE.L	(A0),D0
	MOVE.W	4(A0),D0
	FMOVE.X	FP0,XPWR(a6)		; <1/7/91, JPO>
	ANDI.L	#$7FFFFFFF,D0

	CMPI.L	#$3FB98000,D0		...|X| >= 2**(-70)?
	BGE.B	TWOOK1
	BRA.W	EXPBORS

TWOOK1:
	CMPI.L	#$400D80C0,D0		...|X| > 16480?
	BLE.B	TWOMAIN
	BRA.W	EXPBORS
	

TWOMAIN:
*--USUAL CASE, 2^(-70) <= |X| <= 16480

	FMOVE.X	FP0,FP1
	FMUL.S	#"$42800000",FP1  ...64 * X
	
	FMOVE.L	FP1,NPWR(a6)	...N = ROUND-TO-INT(64 X)  <1/7/91, JPO>

; filter out case where X is integral <3/28/91, JPO>						<T3> thru next <T3>
	fmove.l	fpsr,d0		; check for inexact conversion <3/28/91, JPO>
	andi.w	#$0200,d0	; <3/28/91, JPO>
	bne.b	TWOMAIN1	; inexact <3/28/91, JPO>

	move.l	NPWR(a6),d0	; check if conversion result is a multiple of 64 <3/28/91, JPO>
	bftst	d0{26:6}	; <3/28/91, JPO>
	bne.b	TWOMAIN1	; no <3/28/91, JPO>

	asr.l	#6,d0		; yes. get integer equivalent of X <3/28/91, JPO>
	add.l	#$3fff,d0	; add bias <3/28/91, JPO>
	bmi.b	@1		; extended subnormal <3/28/91, JPO>

	cmp.l	#$7fff,d0	; overflow? <3/28/91, JPO>
	bge.b	TWOMAIN1	; yes <3/28/91, JPO>

	move.w	d0,XPWR(a6)		; normal exact result <3/28/91, JPO>
	move.w	#$3fff,ADJFACT(a6)	; adjust factor is unity <3/28/91, JPO>
	bra.b	@2			; continue below <3/28/91, JPO>
@1:				; subnormal result (may underflow) <3/28/91, JPO>
	add.l	#$3fff,d0	; second bias <3/28/91, JPO>
	clr.l	XPWR(a6)	; result is smallest positive normal <3/28/91, JPO>
	move.w	d0,ADJFACT(a6)	; adjust factor is denorm factor <3/28/91, JPO>

@2:					; result may depend on rounding modes <3/28/91, JPO>
	move.l	#$80000000,XPWR+4(a6)	; prepare rest of result <3/28/91, JPO>
	clr.l	XPWR+8(a6)		; <3/28/91, JPO>
	fmove.x	XPWR(a6),fp0		; fp0 <- result <3/28/91, JPO>
	move.l	#$80000000,ADJFACT+4(a6) ; prepare rest of adjust factor <3/28/91, JPO>
	clr.l	ADJFACT+8(a6)		; <3/28/91, JPO>
	fmove.l	d1,FPCR			; restore user's rounding <3/28/91, JPO>
	fmul.x	ADJFACT(a6),fp0		; final result to fp0 <3/28/91, JPO>
	rts				; return <3/28/91, JPO>

TWOMAIN1:			; label ADDED <3/28/91, JPO>					<T3>
	MOVE.L	d2,-(sp)
	LEA	EXP2TBL,a1 	...LOAD ADDRESS OF TABLE OF 2^(J/64)  <1/7/91, JPO>
	FMOVE.L	NPWR(a6),FP1	...N --> FLOATING FMT  <1/7/91, JPO>
	MOVE.L	NPWR(a6),D0	; <1/7/91, JPO>
	MOVE.L	D0,d2
	ANDI.L	#$3F,D0		...D0 IS J
	ASL.L	#4,D0		...DISPLACEMENT FOR 2^(J/64)
	ADDA.L	D0,a1		...ADDRESS FOR 2^(J/64)
	ASR.L	#6,d2		...d2 IS L, N = 64L + J
	MOVE.L	d2,D0
	ASR.L	#1,D0		...D0 IS M
	SUB.L	D0,d2		...d2 IS M', N = 64(M+M') + J
	ADDI.L	#$3FFF,d2
	MOVE.W	d2,ADJFACT(a6) 	...ADJFACT IS 2^(M')
	MOVE.L	(sp)+,d2
*--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
*--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
*--ADJFACT = 2^(M').
*--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.

	FMUL.S	#"$3C800000",FP1	...(1/64)*N
	MOVE.L	(a1)+,FACT1(a6)
	MOVE.L	(a1)+,FACT1HI(a6)
	MOVE.L	(a1)+,FACT1LOW(a6)
	MOVE.W	(a1)+,FACT2(a6)
	clr.w	FACT2+2(a6)

	FSUB.X	FP1,FP0	 		...X - (1/64)*INT(64 X)

	MOVE.W	(a1)+,FACT2HI(a6)
	clr.w	FACT2HI+2(a6)
	clr.l	FACT2LOW(a6)
	ADD.W	D0,FACT1(a6)
	
	FMUL.X	LOG2,FP0		...FP0 IS R
	ADD.W	D0,FACT2(a6)

	BRA.W	expr

EXPBORS:
*--FPCR, D0 SAVED
	CMPI.L	#$3FFF8000,D0
	BGT.B	EXPBIG2		; label RENAMED <1/7/91, JPO>

;EXPSM:				; label DELETED <1/7/91, JPO>
*--|X| IS SMALL, RETURN 1 + X

	FMOVE.L	d1,FPCR			;restore users exceptions
	FADD.S	#"$3F800000",FP0	...RETURN 1 + X

	bra	t_frcinx

EXPBIG2:				; label RENAMED <1/7/91, JPO>
*--|X| IS LARGE, GENERATE OVERFLOW IF X > 0; ELSE GENERATE UNDERFLOW
*--REGISTERS SAVE SO FAR ARE FPCR AND  D0
	MOVE.L	XPWR(a6),D0	; <1/7/91, JPO>
;	CMPI.L	#0,D0		; DELETED <3/28/91, JPO>			<T3>
;	BLT.B	EXPNEG		; DELETED <3/28/91, JPO>			<T3>
	tst.l	d0		; <3/28/91, JPO>				<T3>
	bmi.b	EXPNEG		; <3/28/91, JPO>				<T3>

	bclr.b	#7,(a0)		;t_ovfl expects positive value
	bra	t_ovfl

EXPNEG:
	bclr.b	#7,(a0)		;t_unfl expects positive value
	bra	t_unfl


stentoxd:
*--ENTRY POINT FOR 10**(X) FOR DENORMALIZED ARGUMENT

	fmove.l	d1,fpcr			...set user's rounding mode/precision
	Fmove.S	#"$3F800000",FP0	...RETURN 1 + X
	move.l	(a0),d0
	or.l	#$00800001,d0
	fadd.s	d0,fp0
	bra	t_frcinx


stentox:
*--ENTRY POINT FOR 10**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
	FMOVEM.X (a0),FP0	...LOAD INPUT, do not set cc's

	MOVE.L	(A0),D0
	MOVE.W	4(A0),D0
	FMOVE.X	FP0,XPWR(a6)		; <1/7/91, JPO>
	ANDI.L	#$7FFFFFFF,D0

	CMPI.L	#$3FB98000,D0		...|X| >= 2**(-70)?
	BGE.B	TENOK1
	BRA.B	EXPBORS

TENOK1:
	CMPI.L	#$400B9B07,D0		...|X| <= 16480*log2/log10 ?
	BLE.B	TENMAIN
	BRA.B	EXPBORS

TENMAIN:
*--USUAL CASE, 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10

	FMOVE.X	FP0,FP1
	FMUL.D	L2TEN64,FP1	...X*64*LOG10/LOG2
	
	FMOVE.L	FP1,NPWR(a6)	...N=INT(X*64*LOG10/LOG2)  <1/7/91, JPO>
	MOVE.L	d2,-(sp)
	LEA	EXP2TBL,a1 	...LOAD ADDRESS OF TABLE OF 2^(J/64)  <1/7/91, JPO>
	FMOVE.L	NPWR(a6),FP1	...N --> FLOATING FMT  <1/7/91, JPO>
	MOVE.L	NPWR(a6),D0	; <1/7/91, JPO>
	MOVE.L	D0,d2
	ANDI.L	#$3F,D0		...D0 IS J
	ASL.L	#4,D0		...DISPLACEMENT FOR 2^(J/64)
	ADDA.L	D0,a1		...ADDRESS FOR 2^(J/64)
	ASR.L	#6,d2		...d2 IS L, N = 64L + J
	MOVE.L	d2,D0
	ASR.L	#1,D0		...D0 IS M
	SUB.L	D0,d2		...d2 IS M', N = 64(M+M') + J
	ADDI.L	#$3FFF,d2
	MOVE.W	d2,ADJFACT(a6) 	...ADJFACT IS 2^(M')
	MOVE.L	(sp)+,d2

*--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
*--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
*--ADJFACT = 2^(M').
*--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.

	FMOVE.X	FP1,FP2

	FMUL.D	L10TWO1,FP1		...N*(LOG2/64LOG10)_LEAD
	MOVE.L	(a1)+,FACT1(a6)

	FMUL.X	L10TWO2,FP2		...N*(LOG2/64LOG10)_TRAIL

	MOVE.L	(a1)+,FACT1HI(a6)
	MOVE.L	(a1)+,FACT1LOW(a6)
	FSUB.X	FP1,FP0			...X - N L_LEAD
	MOVE.W	(a1)+,FACT2(a6)

	FSUB.X	FP2,FP0			...X - N L_TRAIL

	clr.w	FACT2+2(a6)
	MOVE.W	(a1)+,FACT2HI(a6)
	clr.w	FACT2HI+2(a6)
	clr.l	FACT2LOW(a6)

	FMUL.X	LOG10,FP0		...FP0 IS R
	
	ADD.W	D0,FACT1(a6)
	ADD.W	D0,FACT2(a6)

expr:
*--FPCR, FP2, FP3 ARE SAVED IN ORDER AS SHOWN.
*--ADJFACT CONTAINS 2**(M'), FACT1 + FACT2 = 2**(M) * 2**(J/64).
*--FP0 IS R. THE FOLLOWING CODE COMPUTES
*--	2**(M'+M) * 2**(J/64) * EXP(R)

	FMOVE.X	FP0,FP1
	FMUL.X	FP1,FP1		...FP1 IS S = R*R

	FMOVE.D	EXPT5,FP2	...FP2 IS A5  <1/7/91, JPO>
	FMOVE.D	EXPT4,FP3	...FP3 IS A4  <1/7/91, JPO>

	FMUL.X	FP1,FP2		...FP2 IS S*A5
	FMUL.X	FP1,FP3		...FP3 IS S*A4

	FADD.D	EXPT3,FP2	...FP2 IS A3+S*A5  <1/7/91, JPO>
	FADD.D	EXPT2,FP3	...FP3 IS A2+S*A4  <1/7/91, JPO>

	FMUL.X	FP1,FP2		...FP2 IS S*(A3+S*A5)
	FMUL.X	FP1,FP3		...FP3 IS S*(A2+S*A4)

	FADD.D	EXPT1,FP2	...FP2 IS A1+S*(A3+S*A5)  <1/7/91, JPO>
	FMUL.X	FP0,FP3		...FP3 IS R*S*(A2+S*A4)

	FMUL.X	FP1,FP2		...FP2 IS S*(A1+S*(A3+S*A5))
	FADD.X	FP3,FP0		...FP0 IS R+R*S*(A2+S*A4)
	
	FADD.X	FP2,FP0		...FP0 IS EXP(R) - 1
	

*--FINAL RECONSTRUCTION PROCESS
*--EXP(X) = 2^M*2^(J/64) + 2^M*2^(J/64)*(EXP(R)-1)  -  (1 OR 0)

	FMUL.X	FACT1(a6),FP0
	FADD.X	FACT2(a6),FP0
	FADD.X	FACT1(a6),FP0

	FMOVE.L	d1,FPCR		;restore users exceptions
	clr.w	ADJFACT+2(a6)
	move.l	#$80000000,ADJFACT+4(a6)
	clr.l	ADJFACT+8(a6)
	FMUL.X	ADJFACT(a6),FP0	...FINAL ADJUSTMENT

	bra	t_frcinx