xref: /src/sys/arch/m68k/fpsp/setox.sa

*	$NetBSD: setox.sa,v 1.4.152.1 2017/12/03 11:36:23 jdolecek Exp $

*	MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
*	M68000 Hi-Performance Microprocessor Division
*	M68040 Software Package 
*
*	M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
*	All rights reserved.
*
*	THE SOFTWARE is provided on an "AS IS" basis and without warranty.
*	To the maximum extent permitted by applicable law,
*	MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
*	INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
*	PARTICULAR PURPOSE and any warranty against infringement with
*	regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
*	and any accompanying written materials. 
*
*	To the maximum extent permitted by applicable law,
*	IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
*	(INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
*	PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
*	OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
*	SOFTWARE.  Motorola assumes no responsibility for the maintenance
*	and support of the SOFTWARE.  
*
*	You are hereby granted a copyright license to use, modify, and
*	distribute the SOFTWARE so long as this entire notice is retained
*	without alteration in any modified and/or redistributed versions,
*	and that such modified versions are clearly identified as such.
*	No licenses are granted by implication, estoppel or otherwise
*	under any patents or trademarks of Motorola, Inc.

*
*	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 use 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 use 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 use 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.
*

setox	IDNT	2,1 Motorola 040 Floating Point Software Package

	section	8

	include	fpsp.h

L2	DC.L	$3FDC0000,$82E30865,$4361C4C6,$00000000

EXPA3	DC.L	$3FA55555,$55554431
EXPA2	DC.L	$3FC55555,$55554018

HUGE	DC.L	$7FFE0000,$FFFFFFFF,$FFFFFFFF,$00000000
TINY	DC.L	$00010000,$FFFFFFFF,$FFFFFFFF,$00000000

EM1A4	DC.L	$3F811111,$11174385
EM1A3	DC.L	$3FA55555,$55554F5A

EM1A2	DC.L	$3FC55555,$55555555,$00000000,$00000000

EM1B8	DC.L	$3EC71DE3,$A5774682
EM1B7	DC.L	$3EFA01A0,$19D7CB68

EM1B6	DC.L	$3F2A01A0,$1A019DF3
EM1B5	DC.L	$3F56C16C,$16C170E2

EM1B4	DC.L	$3F811111,$11111111
EM1B3	DC.L	$3FA55555,$55555555

EM1B2	DC.L	$3FFC0000,$AAAAAAAA,$AAAAAAAB
	DC.L	$00000000

TWO140	DC.L	$48B00000,$00000000
TWON140	DC.L	$37300000,$00000000

EXPTBL
	DC.L	$3FFF0000,$80000000,$00000000,$00000000
	DC.L	$3FFF0000,$8164D1F3,$BC030774,$9F841A9B
	DC.L	$3FFF0000,$82CD8698,$AC2BA1D8,$9FC1D5B9
	DC.L	$3FFF0000,$843A28C3,$ACDE4048,$A0728369
	DC.L	$3FFF0000,$85AAC367,$CC487B14,$1FC5C95C
	DC.L	$3FFF0000,$871F6196,$9E8D1010,$1EE85C9F
	DC.L	$3FFF0000,$88980E80,$92DA8528,$9FA20729
	DC.L	$3FFF0000,$8A14D575,$496EFD9C,$A07BF9AF
	DC.L	$3FFF0000,$8B95C1E3,$EA8BD6E8,$A0020DCF
	DC.L	$3FFF0000,$8D1ADF5B,$7E5BA9E4,$205A63DA
	DC.L	$3FFF0000,$8EA4398B,$45CD53C0,$1EB70051
	DC.L	$3FFF0000,$9031DC43,$1466B1DC,$1F6EB029
	DC.L	$3FFF0000,$91C3D373,$AB11C338,$A0781494
	DC.L	$3FFF0000,$935A2B2F,$13E6E92C,$9EB319B0
	DC.L	$3FFF0000,$94F4EFA8,$FEF70960,$2017457D
	DC.L	$3FFF0000,$96942D37,$20185A00,$1F11D537
	DC.L	$3FFF0000,$9837F051,$8DB8A970,$9FB952DD
	DC.L	$3FFF0000,$99E04593,$20B7FA64,$1FE43087
	DC.L	$3FFF0000,$9B8D39B9,$D54E5538,$1FA2A818
	DC.L	$3FFF0000,$9D3ED9A7,$2CFFB750,$1FDE494D
	DC.L	$3FFF0000,$9EF53260,$91A111AC,$20504890
	DC.L	$3FFF0000,$A0B0510F,$B9714FC4,$A073691C
	DC.L	$3FFF0000,$A2704303,$0C496818,$1F9B7A05
	DC.L	$3FFF0000,$A43515AE,$09E680A0,$A0797126
	DC.L	$3FFF0000,$A5FED6A9,$B15138EC,$A071A140
	DC.L	$3FFF0000,$A7CD93B4,$E9653568,$204F62DA
	DC.L	$3FFF0000,$A9A15AB4,$EA7C0EF8,$1F283C4A
	DC.L	$3FFF0000,$AB7A39B5,$A93ED338,$9F9A7FDC
	DC.L	$3FFF0000,$AD583EEA,$42A14AC8,$A05B3FAC
	DC.L	$3FFF0000,$AF3B78AD,$690A4374,$1FDF2610
	DC.L	$3FFF0000,$B123F581,$D2AC2590,$9F705F90
	DC.L	$3FFF0000,$B311C412,$A9112488,$201F678A
	DC.L	$3FFF0000,$B504F333,$F9DE6484,$1F32FB13
	DC.L	$3FFF0000,$B6FD91E3,$28D17790,$20038B30
	DC.L	$3FFF0000,$B8FBAF47,$62FB9EE8,$200DC3CC
	DC.L	$3FFF0000,$BAFF5AB2,$133E45FC,$9F8B2AE6
	DC.L	$3FFF0000,$BD08A39F,$580C36C0,$A02BBF70
	DC.L	$3FFF0000,$BF1799B6,$7A731084,$A00BF518
	DC.L	$3FFF0000,$C12C4CCA,$66709458,$A041DD41
	DC.L	$3FFF0000,$C346CCDA,$24976408,$9FDF137B
	DC.L	$3FFF0000,$C5672A11,$5506DADC,$201F1568
	DC.L	$3FFF0000,$C78D74C8,$ABB9B15C,$1FC13A2E
	DC.L	$3FFF0000,$C9B9BD86,$6E2F27A4,$A03F8F03
	DC.L	$3FFF0000,$CBEC14FE,$F2727C5C,$1FF4907D
	DC.L	$3FFF0000,$CE248C15,$1F8480E4,$9E6E53E4
	DC.L	$3FFF0000,$D06333DA,$EF2B2594,$1FD6D45C
	DC.L	$3FFF0000,$D2A81D91,$F12AE45C,$A076EDB9
	DC.L	$3FFF0000,$D4F35AAB,$CFEDFA20,$9FA6DE21
	DC.L	$3FFF0000,$D744FCCA,$D69D6AF4,$1EE69A2F
	DC.L	$3FFF0000,$D99D15C2,$78AFD7B4,$207F439F
	DC.L	$3FFF0000,$DBFBB797,$DAF23754,$201EC207
	DC.L	$3FFF0000,$DE60F482,$5E0E9124,$9E8BE175
	DC.L	$3FFF0000,$E0CCDEEC,$2A94E110,$20032C4B
	DC.L	$3FFF0000,$E33F8972,$BE8A5A50,$2004DFF5
	DC.L	$3FFF0000,$E5B906E7,$7C8348A8,$1E72F47A
	DC.L	$3FFF0000,$E8396A50,$3C4BDC68,$1F722F22
	DC.L	$3FFF0000,$EAC0C6E7,$DD243930,$A017E945
	DC.L	$3FFF0000,$ED4F301E,$D9942B84,$1F401A5B
	DC.L	$3FFF0000,$EFE4B99B,$DCDAF5CC,$9FB9A9E3
	DC.L	$3FFF0000,$F281773C,$59FFB138,$20744C05
	DC.L	$3FFF0000,$F5257D15,$2486CC2C,$1F773A19
	DC.L	$3FFF0000,$F7D0DF73,$0AD13BB8,$1FFE90D5
	DC.L	$3FFF0000,$FA83B2DB,$722A033C,$A041ED22
	DC.L	$3FFF0000,$FD3E0C0C,$F486C174,$1F853F3A

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

	xref	t_frcinx
	xref	t_extdnrm
	xref	t_unfl
	xref	t_ovfl

	xdef	setoxd
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

	xdef	setox
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
	CLR.L		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 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 USE 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
*	CLR.W		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
	TST.L		d0
	BLT		t_unfl
	BRA		t_ovfl

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


	xdef	setoxm1
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 USE 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
*	CLR.W		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)
	BLT.B		EM12TINY
*--Step 8.2
	MOVE.L		#$80010000,SC(a6)	...SC is -2^(-16382)
	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)
	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
	TST.L		d0
	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

	end