Home | History | Annotate | Line # | Download | only in fpsp
setox.sa revision 1.4
      1 *	$NetBSD: setox.sa,v 1.4 2003/02/05 00:02:35 perry Exp $
      2 
      3 *	MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
      4 *	M68000 Hi-Performance Microprocessor Division
      5 *	M68040 Software Package 
      6 *
      7 *	M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
      8 *	All rights reserved.
      9 *
     10 *	THE SOFTWARE is provided on an "AS IS" basis and without warranty.
     11 *	To the maximum extent permitted by applicable law,
     12 *	MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
     13 *	INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
     14 *	PARTICULAR PURPOSE and any warranty against infringement with
     15 *	regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
     16 *	and any accompanying written materials. 
     17 *
     18 *	To the maximum extent permitted by applicable law,
     19 *	IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
     20 *	(INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
     21 *	PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
     22 *	OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
     23 *	SOFTWARE.  Motorola assumes no responsibility for the maintenance
     24 *	and support of the SOFTWARE.  
     25 *
     26 *	You are hereby granted a copyright license to use, modify, and
     27 *	distribute the SOFTWARE so long as this entire notice is retained
     28 *	without alteration in any modified and/or redistributed versions,
     29 *	and that such modified versions are clearly identified as such.
     30 *	No licenses are granted by implication, estoppel or otherwise
     31 *	under any patents or trademarks of Motorola, Inc.
     32 
     33 *
     34 *	setox.sa 3.1 12/10/90
     35 *
     36 *	The entry point setox computes the exponential of a value.
     37 *	setoxd does the same except the input value is a denormalized
     38 *	number.	setoxm1 computes exp(X)-1, and setoxm1d computes
     39 *	exp(X)-1 for denormalized X.
     40 *
     41 *	INPUT
     42 *	-----
     43 *	Double-extended value in memory location pointed to by address
     44 *	register a0.
     45 *
     46 *	OUTPUT
     47 *	------
     48 *	exp(X) or exp(X)-1 returned in floating-point register fp0.
     49 *
     50 *	ACCURACY and MONOTONICITY
     51 *	-------------------------
     52 *	The returned result is within 0.85 ulps in 64 significant bit, i.e.
     53 *	within 0.5001 ulp to 53 bits if the result is subsequently rounded
     54 *	to double precision. The result is provably monotonic in double
     55 *	precision.
     56 *
     57 *	SPEED
     58 *	-----
     59 *	Two timings are measured, both in the copy-back mode. The
     60 *	first one is measured when the function is invoked the first time
     61 *	(so the instructions and data are not in cache), and the
     62 *	second one is measured when the function is reinvoked at the same
     63 *	input argument.
     64 *
     65 *	The program setox takes approximately 210/190 cycles for input
     66 *	argument X whose magnitude is less than 16380 log2, which
     67 *	is the usual situation.	For the less common arguments,
     68 *	depending on their values, the program may run faster or slower --
     69 *	but no worse than 10% slower even in the extreme cases.
     70 *
     71 *	The program setoxm1 takes approximately ???/??? cycles for input
     72 *	argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes
     73 *	approximately ???/??? cycles. For the less common arguments,
     74 *	depending on their values, the program may run faster or slower --
     75 *	but no worse than 10% slower even in the extreme cases.
     76 *
     77 *	ALGORITHM and IMPLEMENTATION NOTES
     78 *	----------------------------------
     79 *
     80 *	setoxd
     81 *	------
     82 *	Step 1.	Set ans := 1.0
     83 *
     84 *	Step 2.	Return	ans := ans + sign(X)*2^(-126). Exit.
     85 *	Notes:	This will always generate one exception -- inexact.
     86 *
     87 *
     88 *	setox
     89 *	-----
     90 *
     91 *	Step 1.	Filter out extreme cases of input argument.
     92 *		1.1	If |X| >= 2^(-65), go to Step 1.3.
     93 *		1.2	Go to Step 7.
     94 *		1.3	If |X| < 16380 log(2), go to Step 2.
     95 *		1.4	Go to Step 8.
     96 *	Notes:	The usual case should take the branches 1.1 -> 1.3 -> 2.
     97 *		 To avoid the use of floating-point comparisons, a
     98 *		 compact representation of |X| is used. This format is a
     99 *		 32-bit integer, the upper (more significant) 16 bits are
    100 *		 the sign and biased exponent field of |X|; the lower 16
    101 *		 bits are the 16 most significant fraction (including the
    102 *		 explicit bit) bits of |X|. Consequently, the comparisons
    103 *		 in Steps 1.1 and 1.3 can be performed by integer comparison.
    104 *		 Note also that the constant 16380 log(2) used in Step 1.3
    105 *		 is also in the compact form. Thus taking the branch
    106 *		 to Step 2 guarantees |X| < 16380 log(2). There is no harm
    107 *		 to have a small number of cases where |X| is less than,
    108 *		 but close to, 16380 log(2) and the branch to Step 9 is
    109 *		 taken.
    110 *
    111 *	Step 2.	Calculate N = round-to-nearest-int( X * 64/log2 ).
    112 *		2.1	Set AdjFlag := 0 (indicates the branch 1.3 -> 2 was taken)
    113 *		2.2	N := round-to-nearest-integer( X * 64/log2 ).
    114 *		2.3	Calculate	J = N mod 64; so J = 0,1,2,..., or 63.
    115 *		2.4	Calculate	M = (N - J)/64; so N = 64M + J.
    116 *		2.5	Calculate the address of the stored value of 2^(J/64).
    117 *		2.6	Create the value Scale = 2^M.
    118 *	Notes:	The calculation in 2.2 is really performed by
    119 *
    120 *			Z := X * constant
    121 *			N := round-to-nearest-integer(Z)
    122 *
    123 *		 where
    124 *
    125 *			constant := single-precision( 64/log 2 ).
    126 *
    127 *		 Using a single-precision constant avoids memory access.
    128 *		 Another effect of using a single-precision "constant" is
    129 *		 that the calculated value Z is
    130 *
    131 *			Z = X*(64/log2)*(1+eps), |eps| <= 2^(-24).
    132 *
    133 *		 This error has to be considered later in Steps 3 and 4.
    134 *
    135 *	Step 3.	Calculate X - N*log2/64.
    136 *		3.1	R := X + N*L1, where L1 := single-precision(-log2/64).
    137 *		3.2	R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
    138 *	Notes:	a) The way L1 and L2 are chosen ensures L1+L2 approximate
    139 *		 the value	-log2/64	to 88 bits of accuracy.
    140 *		 b) N*L1 is exact because N is no longer than 22 bits and
    141 *		 L1 is no longer than 24 bits.
    142 *		 c) The calculation X+N*L1 is also exact due to cancellation.
    143 *		 Thus, R is practically X+N(L1+L2) to full 64 bits.
    144 *		 d) It is important to estimate how large can |R| be after
    145 *		 Step 3.2.
    146 *
    147 *			N = rnd-to-int( X*64/log2 (1+eps) ), |eps|<=2^(-24)
    148 *			X*64/log2 (1+eps)	=	N + f,	|f| <= 0.5
    149 *			X*64/log2 - N	=	f - eps*X 64/log2
    150 *			X - N*log2/64	=	f*log2/64 - eps*X
    151 *
    152 *
    153 *		 Now |X| <= 16446 log2, thus
    154 *
    155 *			|X - N*log2/64| <= (0.5 + 16446/2^(18))*log2/64
    156 *					<= 0.57 log2/64.
    157 *		 This bound will be used in Step 4.
    158 *
    159 *	Step 4.	Approximate exp(R)-1 by a polynomial
    160 *			p = R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
    161 *	Notes:	a) In order to reduce memory access, the coefficients are
    162 *		 made as "short" as possible: A1 (which is 1/2), A4 and A5
    163 *		 are single precision; A2 and A3 are double precision.
    164 *		 b) Even with the restrictions above,
    165 *			|p - (exp(R)-1)| < 2^(-68.8) for all |R| <= 0.0062.
    166 *		 Note that 0.0062 is slightly bigger than 0.57 log2/64.
    167 *		 c) To fully use the pipeline, p is separated into
    168 *		 two independent pieces of roughly equal complexities
    169 *			p = [ R + R*S*(A2 + S*A4) ]	+
    170 *				[ S*(A1 + S*(A3 + S*A5)) ]
    171 *		 where S = R*R.
    172 *
    173 *	Step 5.	Compute 2^(J/64)*exp(R) = 2^(J/64)*(1+p) by
    174 *				ans := T + ( T*p + t)
    175 *		 where T and t are the stored values for 2^(J/64).
    176 *	Notes:	2^(J/64) is stored as T and t where T+t approximates
    177 *		 2^(J/64) to roughly 85 bits; T is in extended precision
    178 *		 and t is in single precision. Note also that T is rounded
    179 *		 to 62 bits so that the last two bits of T are zero. The
    180 *		 reason for such a special form is that T-1, T-2, and T-8
    181 *		 will all be exact --- a property that will give much
    182 *		 more accurate computation of the function EXPM1.
    183 *
    184 *	Step 6.	Reconstruction of exp(X)
    185 *			exp(X) = 2^M * 2^(J/64) * exp(R).
    186 *		6.1	If AdjFlag = 0, go to 6.3
    187 *		6.2	ans := ans * AdjScale
    188 *		6.3	Restore the user FPCR
    189 *		6.4	Return ans := ans * Scale. Exit.
    190 *	Notes:	If AdjFlag = 0, we have X = Mlog2 + Jlog2/64 + R,
    191 *		 |M| <= 16380, and Scale = 2^M. Moreover, exp(X) will
    192 *		 neither overflow nor underflow. If AdjFlag = 1, that
    193 *		 means that
    194 *			X = (M1+M)log2 + Jlog2/64 + R, |M1+M| >= 16380.
    195 *		 Hence, exp(X) may overflow or underflow or neither.
    196 *		 When that is the case, AdjScale = 2^(M1) where M1 is
    197 *		 approximately M. Thus 6.2 will never cause over/underflow.
    198 *		 Possible exception in 6.4 is overflow or underflow.
    199 *		 The inexact exception is not generated in 6.4. Although
    200 *		 one can argue that the inexact flag should always be
    201 *		 raised, to simulate that exception cost to much than the
    202 *		 flag is worth in practical uses.
    203 *
    204 *	Step 7.	Return 1 + X.
    205 *		7.1	ans := X
    206 *		7.2	Restore user FPCR.
    207 *		7.3	Return ans := 1 + ans. Exit
    208 *	Notes:	For non-zero X, the inexact exception will always be
    209 *		 raised by 7.3. That is the only exception raised by 7.3.
    210 *		 Note also that we use the FMOVEM instruction to move X
    211 *		 in Step 7.1 to avoid unnecessary trapping. (Although
    212 *		 the FMOVEM may not seem relevant since X is normalized,
    213 *		 the precaution will be useful in the library version of
    214 *		 this code where the separate entry for denormalized inputs
    215 *		 will be done away with.)
    216 *
    217 *	Step 8.	Handle exp(X) where |X| >= 16380log2.
    218 *		8.1	If |X| > 16480 log2, go to Step 9.
    219 *		(mimic 2.2 - 2.6)
    220 *		8.2	N := round-to-integer( X * 64/log2 )
    221 *		8.3	Calculate J = N mod 64, J = 0,1,...,63
    222 *		8.4	K := (N-J)/64, M1 := truncate(K/2), M = K-M1, AdjFlag := 1.
    223 *		8.5	Calculate the address of the stored value 2^(J/64).
    224 *		8.6	Create the values Scale = 2^M, AdjScale = 2^M1.
    225 *		8.7	Go to Step 3.
    226 *	Notes:	Refer to notes for 2.2 - 2.6.
    227 *
    228 *	Step 9.	Handle exp(X), |X| > 16480 log2.
    229 *		9.1	If X < 0, go to 9.3
    230 *		9.2	ans := Huge, go to 9.4
    231 *		9.3	ans := Tiny.
    232 *		9.4	Restore user FPCR.
    233 *		9.5	Return ans := ans * ans. Exit.
    234 *	Notes:	Exp(X) will surely overflow or underflow, depending on
    235 *		 X's sign. "Huge" and "Tiny" are respectively large/tiny
    236 *		 extended-precision numbers whose square over/underflow
    237 *		 with an inexact result. Thus, 9.5 always raises the
    238 *		 inexact together with either overflow or underflow.
    239 *
    240 *
    241 *	setoxm1d
    242 *	--------
    243 *
    244 *	Step 1.	Set ans := 0
    245 *
    246 *	Step 2.	Return	ans := X + ans. Exit.
    247 *	Notes:	This will return X with the appropriate rounding
    248 *		 precision prescribed by the user FPCR.
    249 *
    250 *	setoxm1
    251 *	-------
    252 *
    253 *	Step 1.	Check |X|
    254 *		1.1	If |X| >= 1/4, go to Step 1.3.
    255 *		1.2	Go to Step 7.
    256 *		1.3	If |X| < 70 log(2), go to Step 2.
    257 *		1.4	Go to Step 10.
    258 *	Notes:	The usual case should take the branches 1.1 -> 1.3 -> 2.
    259 *		 However, it is conceivable |X| can be small very often
    260 *		 because EXPM1 is intended to evaluate exp(X)-1 accurately
    261 *		 when |X| is small. For further details on the comparisons,
    262 *		 see the notes on Step 1 of setox.
    263 *
    264 *	Step 2.	Calculate N = round-to-nearest-int( X * 64/log2 ).
    265 *		2.1	N := round-to-nearest-integer( X * 64/log2 ).
    266 *		2.2	Calculate	J = N mod 64; so J = 0,1,2,..., or 63.
    267 *		2.3	Calculate	M = (N - J)/64; so N = 64M + J.
    268 *		2.4	Calculate the address of the stored value of 2^(J/64).
    269 *		2.5	Create the values Sc = 2^M and OnebySc := -2^(-M).
    270 *	Notes:	See the notes on Step 2 of setox.
    271 *
    272 *	Step 3.	Calculate X - N*log2/64.
    273 *		3.1	R := X + N*L1, where L1 := single-precision(-log2/64).
    274 *		3.2	R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
    275 *	Notes:	Applying the analysis of Step 3 of setox in this case
    276 *		 shows that |R| <= 0.0055 (note that |X| <= 70 log2 in
    277 *		 this case).
    278 *
    279 *	Step 4.	Approximate exp(R)-1 by a polynomial
    280 *			p = R+R*R*(A1+R*(A2+R*(A3+R*(A4+R*(A5+R*A6)))))
    281 *	Notes:	a) In order to reduce memory access, the coefficients are
    282 *		 made as "short" as possible: A1 (which is 1/2), A5 and A6
    283 *		 are single precision; A2, A3 and A4 are double precision.
    284 *		 b) Even with the restriction above,
    285 *			|p - (exp(R)-1)| <	|R| * 2^(-72.7)
    286 *		 for all |R| <= 0.0055.
    287 *		 c) To fully use the pipeline, p is separated into
    288 *		 two independent pieces of roughly equal complexity
    289 *			p = [ R*S*(A2 + S*(A4 + S*A6)) ]	+
    290 *				[ R + S*(A1 + S*(A3 + S*A5)) ]
    291 *		 where S = R*R.
    292 *
    293 *	Step 5.	Compute 2^(J/64)*p by
    294 *				p := T*p
    295 *		 where T and t are the stored values for 2^(J/64).
    296 *	Notes:	2^(J/64) is stored as T and t where T+t approximates
    297 *		 2^(J/64) to roughly 85 bits; T is in extended precision
    298 *		 and t is in single precision. Note also that T is rounded
    299 *		 to 62 bits so that the last two bits of T are zero. The
    300 *		 reason for such a special form is that T-1, T-2, and T-8
    301 *		 will all be exact --- a property that will be exploited
    302 *		 in Step 6 below. The total relative error in p is no
    303 *		 bigger than 2^(-67.7) compared to the final result.
    304 *
    305 *	Step 6.	Reconstruction of exp(X)-1
    306 *			exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ).
    307 *		6.1	If M <= 63, go to Step 6.3.
    308 *		6.2	ans := T + (p + (t + OnebySc)). Go to 6.6
    309 *		6.3	If M >= -3, go to 6.5.
    310 *		6.4	ans := (T + (p + t)) + OnebySc. Go to 6.6
    311 *		6.5	ans := (T + OnebySc) + (p + t).
    312 *		6.6	Restore user FPCR.
    313 *		6.7	Return ans := Sc * ans. Exit.
    314 *	Notes:	The various arrangements of the expressions give accurate
    315 *		 evaluations.
    316 *
    317 *	Step 7.	exp(X)-1 for |X| < 1/4.
    318 *		7.1	If |X| >= 2^(-65), go to Step 9.
    319 *		7.2	Go to Step 8.
    320 *
    321 *	Step 8.	Calculate exp(X)-1, |X| < 2^(-65).
    322 *		8.1	If |X| < 2^(-16312), goto 8.3
    323 *		8.2	Restore FPCR; return ans := X - 2^(-16382). Exit.
    324 *		8.3	X := X * 2^(140).
    325 *		8.4	Restore FPCR; ans := ans - 2^(-16382).
    326 *		 Return ans := ans*2^(140). Exit
    327 *	Notes:	The idea is to return "X - tiny" under the user
    328 *		 precision and rounding modes. To avoid unnecessary
    329 *		 inefficiency, we stay away from denormalized numbers the
    330 *		 best we can. For |X| >= 2^(-16312), the straightforward
    331 *		 8.2 generates the inexact exception as the case warrants.
    332 *
    333 *	Step 9.	Calculate exp(X)-1, |X| < 1/4, by a polynomial
    334 *			p = X + X*X*(B1 + X*(B2 + ... + X*B12))
    335 *	Notes:	a) In order to reduce memory access, the coefficients are
    336 *		 made as "short" as possible: B1 (which is 1/2), B9 to B12
    337 *		 are single precision; B3 to B8 are double precision; and
    338 *		 B2 is double extended.
    339 *		 b) Even with the restriction above,
    340 *			|p - (exp(X)-1)| < |X| 2^(-70.6)
    341 *		 for all |X| <= 0.251.
    342 *		 Note that 0.251 is slightly bigger than 1/4.
    343 *		 c) To fully preserve accuracy, the polynomial is computed
    344 *		 as	X + ( S*B1 +	Q ) where S = X*X and
    345 *			Q	=	X*S*(B2 + X*(B3 + ... + X*B12))
    346 *		 d) To fully use the pipeline, Q is separated into
    347 *		 two independent pieces of roughly equal complexity
    348 *			Q = [ X*S*(B2 + S*(B4 + ... + S*B12)) ] +
    349 *				[ S*S*(B3 + S*(B5 + ... + S*B11)) ]
    350 *
    351 *	Step 10.	Calculate exp(X)-1 for |X| >= 70 log 2.
    352 *		10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all practical
    353 *		 purposes. Therefore, go to Step 1 of setox.
    354 *		10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical purposes.
    355 *		 ans := -1
    356 *		 Restore user FPCR
    357 *		 Return ans := ans + 2^(-126). Exit.
    358 *	Notes:	10.2 will always create an inexact and return -1 + tiny
    359 *		 in the user rounding precision and mode.
    360 *
    361 
    362 setox	IDNT	2,1 Motorola 040 Floating Point Software Package
    363 
    364 	section	8
    365 
    366 	include	fpsp.h
    367 
    368 L2	DC.L	$3FDC0000,$82E30865,$4361C4C6,$00000000
    369 
    370 EXPA3	DC.L	$3FA55555,$55554431
    371 EXPA2	DC.L	$3FC55555,$55554018
    372 
    373 HUGE	DC.L	$7FFE0000,$FFFFFFFF,$FFFFFFFF,$00000000
    374 TINY	DC.L	$00010000,$FFFFFFFF,$FFFFFFFF,$00000000
    375 
    376 EM1A4	DC.L	$3F811111,$11174385
    377 EM1A3	DC.L	$3FA55555,$55554F5A
    378 
    379 EM1A2	DC.L	$3FC55555,$55555555,$00000000,$00000000
    380 
    381 EM1B8	DC.L	$3EC71DE3,$A5774682
    382 EM1B7	DC.L	$3EFA01A0,$19D7CB68
    383 
    384 EM1B6	DC.L	$3F2A01A0,$1A019DF3
    385 EM1B5	DC.L	$3F56C16C,$16C170E2
    386 
    387 EM1B4	DC.L	$3F811111,$11111111
    388 EM1B3	DC.L	$3FA55555,$55555555
    389 
    390 EM1B2	DC.L	$3FFC0000,$AAAAAAAA,$AAAAAAAB
    391 	DC.L	$00000000
    392 
    393 TWO140	DC.L	$48B00000,$00000000
    394 TWON140	DC.L	$37300000,$00000000
    395 
    396 EXPTBL
    397 	DC.L	$3FFF0000,$80000000,$00000000,$00000000
    398 	DC.L	$3FFF0000,$8164D1F3,$BC030774,$9F841A9B
    399 	DC.L	$3FFF0000,$82CD8698,$AC2BA1D8,$9FC1D5B9
    400 	DC.L	$3FFF0000,$843A28C3,$ACDE4048,$A0728369
    401 	DC.L	$3FFF0000,$85AAC367,$CC487B14,$1FC5C95C
    402 	DC.L	$3FFF0000,$871F6196,$9E8D1010,$1EE85C9F
    403 	DC.L	$3FFF0000,$88980E80,$92DA8528,$9FA20729
    404 	DC.L	$3FFF0000,$8A14D575,$496EFD9C,$A07BF9AF
    405 	DC.L	$3FFF0000,$8B95C1E3,$EA8BD6E8,$A0020DCF
    406 	DC.L	$3FFF0000,$8D1ADF5B,$7E5BA9E4,$205A63DA
    407 	DC.L	$3FFF0000,$8EA4398B,$45CD53C0,$1EB70051
    408 	DC.L	$3FFF0000,$9031DC43,$1466B1DC,$1F6EB029
    409 	DC.L	$3FFF0000,$91C3D373,$AB11C338,$A0781494
    410 	DC.L	$3FFF0000,$935A2B2F,$13E6E92C,$9EB319B0
    411 	DC.L	$3FFF0000,$94F4EFA8,$FEF70960,$2017457D
    412 	DC.L	$3FFF0000,$96942D37,$20185A00,$1F11D537
    413 	DC.L	$3FFF0000,$9837F051,$8DB8A970,$9FB952DD
    414 	DC.L	$3FFF0000,$99E04593,$20B7FA64,$1FE43087
    415 	DC.L	$3FFF0000,$9B8D39B9,$D54E5538,$1FA2A818
    416 	DC.L	$3FFF0000,$9D3ED9A7,$2CFFB750,$1FDE494D
    417 	DC.L	$3FFF0000,$9EF53260,$91A111AC,$20504890
    418 	DC.L	$3FFF0000,$A0B0510F,$B9714FC4,$A073691C
    419 	DC.L	$3FFF0000,$A2704303,$0C496818,$1F9B7A05
    420 	DC.L	$3FFF0000,$A43515AE,$09E680A0,$A0797126
    421 	DC.L	$3FFF0000,$A5FED6A9,$B15138EC,$A071A140
    422 	DC.L	$3FFF0000,$A7CD93B4,$E9653568,$204F62DA
    423 	DC.L	$3FFF0000,$A9A15AB4,$EA7C0EF8,$1F283C4A
    424 	DC.L	$3FFF0000,$AB7A39B5,$A93ED338,$9F9A7FDC
    425 	DC.L	$3FFF0000,$AD583EEA,$42A14AC8,$A05B3FAC
    426 	DC.L	$3FFF0000,$AF3B78AD,$690A4374,$1FDF2610
    427 	DC.L	$3FFF0000,$B123F581,$D2AC2590,$9F705F90
    428 	DC.L	$3FFF0000,$B311C412,$A9112488,$201F678A
    429 	DC.L	$3FFF0000,$B504F333,$F9DE6484,$1F32FB13
    430 	DC.L	$3FFF0000,$B6FD91E3,$28D17790,$20038B30
    431 	DC.L	$3FFF0000,$B8FBAF47,$62FB9EE8,$200DC3CC
    432 	DC.L	$3FFF0000,$BAFF5AB2,$133E45FC,$9F8B2AE6
    433 	DC.L	$3FFF0000,$BD08A39F,$580C36C0,$A02BBF70
    434 	DC.L	$3FFF0000,$BF1799B6,$7A731084,$A00BF518
    435 	DC.L	$3FFF0000,$C12C4CCA,$66709458,$A041DD41
    436 	DC.L	$3FFF0000,$C346CCDA,$24976408,$9FDF137B
    437 	DC.L	$3FFF0000,$C5672A11,$5506DADC,$201F1568
    438 	DC.L	$3FFF0000,$C78D74C8,$ABB9B15C,$1FC13A2E
    439 	DC.L	$3FFF0000,$C9B9BD86,$6E2F27A4,$A03F8F03
    440 	DC.L	$3FFF0000,$CBEC14FE,$F2727C5C,$1FF4907D
    441 	DC.L	$3FFF0000,$CE248C15,$1F8480E4,$9E6E53E4
    442 	DC.L	$3FFF0000,$D06333DA,$EF2B2594,$1FD6D45C
    443 	DC.L	$3FFF0000,$D2A81D91,$F12AE45C,$A076EDB9
    444 	DC.L	$3FFF0000,$D4F35AAB,$CFEDFA20,$9FA6DE21
    445 	DC.L	$3FFF0000,$D744FCCA,$D69D6AF4,$1EE69A2F
    446 	DC.L	$3FFF0000,$D99D15C2,$78AFD7B4,$207F439F
    447 	DC.L	$3FFF0000,$DBFBB797,$DAF23754,$201EC207
    448 	DC.L	$3FFF0000,$DE60F482,$5E0E9124,$9E8BE175
    449 	DC.L	$3FFF0000,$E0CCDEEC,$2A94E110,$20032C4B
    450 	DC.L	$3FFF0000,$E33F8972,$BE8A5A50,$2004DFF5
    451 	DC.L	$3FFF0000,$E5B906E7,$7C8348A8,$1E72F47A
    452 	DC.L	$3FFF0000,$E8396A50,$3C4BDC68,$1F722F22
    453 	DC.L	$3FFF0000,$EAC0C6E7,$DD243930,$A017E945
    454 	DC.L	$3FFF0000,$ED4F301E,$D9942B84,$1F401A5B
    455 	DC.L	$3FFF0000,$EFE4B99B,$DCDAF5CC,$9FB9A9E3
    456 	DC.L	$3FFF0000,$F281773C,$59FFB138,$20744C05
    457 	DC.L	$3FFF0000,$F5257D15,$2486CC2C,$1F773A19
    458 	DC.L	$3FFF0000,$F7D0DF73,$0AD13BB8,$1FFE90D5
    459 	DC.L	$3FFF0000,$FA83B2DB,$722A033C,$A041ED22
    460 	DC.L	$3FFF0000,$FD3E0C0C,$F486C174,$1F853F3A
    461 
    462 ADJFLAG	equ L_SCR2
    463 SCALE	equ FP_SCR1
    464 ADJSCALE equ FP_SCR2
    465 SC	equ FP_SCR3
    466 ONEBYSC	equ FP_SCR4
    467 
    468 	xref	t_frcinx
    469 	xref	t_extdnrm
    470 	xref	t_unfl
    471 	xref	t_ovfl
    472 
    473 	xdef	setoxd
    474 setoxd:
    475 *--entry point for EXP(X), X is denormalized
    476 	MOVE.L		(a0),d0
    477 	ANDI.L		#$80000000,d0
    478 	ORI.L		#$00800000,d0		...sign(X)*2^(-126)
    479 	MOVE.L		d0,-(sp)
    480 	FMOVE.S		#:3F800000,fp0
    481 	fmove.l		d1,fpcr
    482 	FADD.S		(sp)+,fp0
    483 	bra		t_frcinx
    484 
    485 	xdef	setox
    486 setox:
    487 *--entry point for EXP(X), here X is finite, non-zero, and not NaN's
    488 
    489 *--Step 1.
    490 	MOVE.L		(a0),d0	 ...load part of input X
    491 	ANDI.L		#$7FFF0000,d0	...biased expo. of X
    492 	CMPI.L		#$3FBE0000,d0	...2^(-65)
    493 	BGE.B		EXPC1		...normal case
    494 	BRA.W		EXPSM
    495 
    496 EXPC1:
    497 *--The case |X| >= 2^(-65)
    498 	MOVE.W		4(a0),d0	...expo. and partial sig. of |X|
    499 	CMPI.L		#$400CB167,d0	...16380 log2 trunc. 16 bits
    500 	BLT.B		EXPMAIN	 ...normal case
    501 	BRA.W		EXPBIG
    502 
    503 EXPMAIN:
    504 *--Step 2.
    505 *--This is the normal branch:	2^(-65) <= |X| < 16380 log2.
    506 	FMOVE.X		(a0),fp0	...load input from (a0)
    507 
    508 	FMOVE.X		fp0,fp1
    509 	FMUL.S		#:42B8AA3B,fp0	...64/log2 * X
    510 	fmovem.x	fp2/fp3,-(a7)		...save fp2
    511 	CLR.L		ADJFLAG(a6)
    512 	FMOVE.L		fp0,d0		...N = int( X * 64/log2 )
    513 	LEA		EXPTBL,a1
    514 	FMOVE.L		d0,fp0		...convert to floating-format
    515 
    516 	MOVE.L		d0,L_SCR1(a6)	...save N temporarily
    517 	ANDI.L		#$3F,d0		...D0 is J = N mod 64
    518 	LSL.L		#4,d0
    519 	ADDA.L		d0,a1		...address of 2^(J/64)
    520 	MOVE.L		L_SCR1(a6),d0
    521 	ASR.L		#6,d0		...D0 is M
    522 	ADDI.W		#$3FFF,d0	...biased expo. of 2^(M)
    523 	MOVE.W		L2,L_SCR1(a6)	...prefetch L2, no need in CB
    524 
    525 EXPCONT1:
    526 *--Step 3.
    527 *--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
    528 *--a0 points to 2^(J/64), D0 is biased expo. of 2^(M)
    529 	FMOVE.X		fp0,fp2
    530 	FMUL.S		#:BC317218,fp0	...N * L1, L1 = lead(-log2/64)
    531 	FMUL.X		L2,fp2		...N * L2, L1+L2 = -log2/64
    532 	FADD.X		fp1,fp0	 	...X + N*L1
    533 	FADD.X		fp2,fp0		...fp0 is R, reduced arg.
    534 *	MOVE.W		#$3FA5,EXPA3	...load EXPA3 in cache
    535 
    536 *--Step 4.
    537 *--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
    538 *-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
    539 *--TO FULLY USE THE PIPELINE, WE COMPUTE S = R*R
    540 *--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))]
    541 
    542 	FMOVE.X		fp0,fp1
    543 	FMUL.X		fp1,fp1	 	...fp1 IS S = R*R
    544 
    545 	FMOVE.S		#:3AB60B70,fp2	...fp2 IS A5
    546 *	CLR.W		2(a1)		...load 2^(J/64) in cache
    547 
    548 	FMUL.X		fp1,fp2	 	...fp2 IS S*A5
    549 	FMOVE.X		fp1,fp3
    550 	FMUL.S		#:3C088895,fp3	...fp3 IS S*A4
    551 
    552 	FADD.D		EXPA3,fp2	...fp2 IS A3+S*A5
    553 	FADD.D		EXPA2,fp3	...fp3 IS A2+S*A4
    554 
    555 	FMUL.X		fp1,fp2	 	...fp2 IS S*(A3+S*A5)
    556 	MOVE.W		d0,SCALE(a6)	...SCALE is 2^(M) in extended
    557 	clr.w		SCALE+2(a6)
    558 	move.l		#$80000000,SCALE+4(a6)
    559 	clr.l		SCALE+8(a6)
    560 
    561 	FMUL.X		fp1,fp3	 	...fp3 IS S*(A2+S*A4)
    562 
    563 	FADD.S		#:3F000000,fp2	...fp2 IS A1+S*(A3+S*A5)
    564 	FMUL.X		fp0,fp3	 	...fp3 IS R*S*(A2+S*A4)
    565 
    566 	FMUL.X		fp1,fp2	 	...fp2 IS S*(A1+S*(A3+S*A5))
    567 	FADD.X		fp3,fp0	 	...fp0 IS R+R*S*(A2+S*A4),
    568 *					...fp3 released
    569 
    570 	FMOVE.X		(a1)+,fp1	...fp1 is lead. pt. of 2^(J/64)
    571 	FADD.X		fp2,fp0	 	...fp0 is EXP(R) - 1
    572 *					...fp2 released
    573 
    574 *--Step 5
    575 *--final reconstruction process
    576 *--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) )
    577 
    578 	FMUL.X		fp1,fp0	 	...2^(J/64)*(Exp(R)-1)
    579 	fmovem.x	(a7)+,fp2/fp3	...fp2 restored
    580 	FADD.S		(a1),fp0	...accurate 2^(J/64)
    581 
    582 	FADD.X		fp1,fp0	 	...2^(J/64) + 2^(J/64)*...
    583 	MOVE.L		ADJFLAG(a6),d0
    584 
    585 *--Step 6
    586 	TST.L		D0
    587 	BEQ.B		NORMAL
    588 ADJUST:
    589 	FMUL.X		ADJSCALE(a6),fp0
    590 NORMAL:
    591 	FMOVE.L		d1,FPCR	 	...restore user FPCR
    592 	FMUL.X		SCALE(a6),fp0	...multiply 2^(M)
    593 	bra		t_frcinx
    594 
    595 EXPSM:
    596 *--Step 7
    597 	FMOVEM.X	(a0),fp0	...in case X is denormalized
    598 	FMOVE.L		d1,FPCR
    599 	FADD.S		#:3F800000,fp0	...1+X in user mode
    600 	bra		t_frcinx
    601 
    602 EXPBIG:
    603 *--Step 8
    604 	CMPI.L		#$400CB27C,d0	...16480 log2
    605 	BGT.B		EXP2BIG
    606 *--Steps 8.2 -- 8.6
    607 	FMOVE.X		(a0),fp0	...load input from (a0)
    608 
    609 	FMOVE.X		fp0,fp1
    610 	FMUL.S		#:42B8AA3B,fp0	...64/log2 * X
    611 	fmovem.x	 fp2/fp3,-(a7)		...save fp2
    612 	MOVE.L		#1,ADJFLAG(a6)
    613 	FMOVE.L		fp0,d0		...N = int( X * 64/log2 )
    614 	LEA		EXPTBL,a1
    615 	FMOVE.L		d0,fp0		...convert to floating-format
    616 	MOVE.L		d0,L_SCR1(a6)			...save N temporarily
    617 	ANDI.L		#$3F,d0		 ...D0 is J = N mod 64
    618 	LSL.L		#4,d0
    619 	ADDA.L		d0,a1			...address of 2^(J/64)
    620 	MOVE.L		L_SCR1(a6),d0
    621 	ASR.L		#6,d0			...D0 is K
    622 	MOVE.L		d0,L_SCR1(a6)			...save K temporarily
    623 	ASR.L		#1,d0			...D0 is M1
    624 	SUB.L		d0,L_SCR1(a6)			...a1 is M
    625 	ADDI.W		#$3FFF,d0		...biased expo. of 2^(M1)
    626 	MOVE.W		d0,ADJSCALE(a6)		...ADJSCALE := 2^(M1)
    627 	clr.w		ADJSCALE+2(a6)
    628 	move.l		#$80000000,ADJSCALE+4(a6)
    629 	clr.l		ADJSCALE+8(a6)
    630 	MOVE.L		L_SCR1(a6),d0			...D0 is M
    631 	ADDI.W		#$3FFF,d0		...biased expo. of 2^(M)
    632 	BRA.W		EXPCONT1		...go back to Step 3
    633 
    634 EXP2BIG:
    635 *--Step 9
    636 	FMOVE.L		d1,FPCR
    637 	MOVE.L		(a0),d0
    638 	bclr.b		#sign_bit,(a0)		...setox always returns positive
    639 	TST.L		d0
    640 	BLT		t_unfl
    641 	BRA		t_ovfl
    642 
    643 	xdef	setoxm1d
    644 setoxm1d:
    645 *--entry point for EXPM1(X), here X is denormalized
    646 *--Step 0.
    647 	bra		t_extdnrm
    648 
    649 
    650 	xdef	setoxm1
    651 setoxm1:
    652 *--entry point for EXPM1(X), here X is finite, non-zero, non-NaN
    653 
    654 *--Step 1.
    655 *--Step 1.1
    656 	MOVE.L		(a0),d0	 ...load part of input X
    657 	ANDI.L		#$7FFF0000,d0	...biased expo. of X
    658 	CMPI.L		#$3FFD0000,d0	...1/4
    659 	BGE.B		EM1CON1	 ...|X| >= 1/4
    660 	BRA.W		EM1SM
    661 
    662 EM1CON1:
    663 *--Step 1.3
    664 *--The case |X| >= 1/4
    665 	MOVE.W		4(a0),d0	...expo. and partial sig. of |X|
    666 	CMPI.L		#$4004C215,d0	...70log2 rounded up to 16 bits
    667 	BLE.B		EM1MAIN	 ...1/4 <= |X| <= 70log2
    668 	BRA.W		EM1BIG
    669 
    670 EM1MAIN:
    671 *--Step 2.
    672 *--This is the case:	1/4 <= |X| <= 70 log2.
    673 	FMOVE.X		(a0),fp0	...load input from (a0)
    674 
    675 	FMOVE.X		fp0,fp1
    676 	FMUL.S		#:42B8AA3B,fp0	...64/log2 * X
    677 	fmovem.x	fp2/fp3,-(a7)		...save fp2
    678 *	MOVE.W		#$3F81,EM1A4		...prefetch in CB mode
    679 	FMOVE.L		fp0,d0		...N = int( X * 64/log2 )
    680 	LEA		EXPTBL,a1
    681 	FMOVE.L		d0,fp0		...convert to floating-format
    682 
    683 	MOVE.L		d0,L_SCR1(a6)			...save N temporarily
    684 	ANDI.L		#$3F,d0		 ...D0 is J = N mod 64
    685 	LSL.L		#4,d0
    686 	ADDA.L		d0,a1			...address of 2^(J/64)
    687 	MOVE.L		L_SCR1(a6),d0
    688 	ASR.L		#6,d0			...D0 is M
    689 	MOVE.L		d0,L_SCR1(a6)			...save a copy of M
    690 *	MOVE.W		#$3FDC,L2		...prefetch L2 in CB mode
    691 
    692 *--Step 3.
    693 *--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
    694 *--a0 points to 2^(J/64), D0 and a1 both contain M
    695 	FMOVE.X		fp0,fp2
    696 	FMUL.S		#:BC317218,fp0	...N * L1, L1 = lead(-log2/64)
    697 	FMUL.X		L2,fp2		...N * L2, L1+L2 = -log2/64
    698 	FADD.X		fp1,fp0	 ...X + N*L1
    699 	FADD.X		fp2,fp0	 ...fp0 is R, reduced arg.
    700 *	MOVE.W		#$3FC5,EM1A2		...load EM1A2 in cache
    701 	ADDI.W		#$3FFF,d0		...D0 is biased expo. of 2^M
    702 
    703 *--Step 4.
    704 *--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
    705 *-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6)))))
    706 *--TO FULLY USE THE PIPELINE, WE COMPUTE S = R*R
    707 *--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))]
    708 
    709 	FMOVE.X		fp0,fp1
    710 	FMUL.X		fp1,fp1		...fp1 IS S = R*R
    711 
    712 	FMOVE.S		#:3950097B,fp2	...fp2 IS a6
    713 *	CLR.W		2(a1)		...load 2^(J/64) in cache
    714 
    715 	FMUL.X		fp1,fp2		...fp2 IS S*A6
    716 	FMOVE.X		fp1,fp3
    717 	FMUL.S		#:3AB60B6A,fp3	...fp3 IS S*A5
    718 
    719 	FADD.D		EM1A4,fp2	...fp2 IS A4+S*A6
    720 	FADD.D		EM1A3,fp3	...fp3 IS A3+S*A5
    721 	MOVE.W		d0,SC(a6)		...SC is 2^(M) in extended
    722 	clr.w		SC+2(a6)
    723 	move.l		#$80000000,SC+4(a6)
    724 	clr.l		SC+8(a6)
    725 
    726 	FMUL.X		fp1,fp2		...fp2 IS S*(A4+S*A6)
    727 	MOVE.L		L_SCR1(a6),d0		...D0 is	M
    728 	NEG.W		D0		...D0 is -M
    729 	FMUL.X		fp1,fp3		...fp3 IS S*(A3+S*A5)
    730 	ADDI.W		#$3FFF,d0	...biased expo. of 2^(-M)
    731 	FADD.D		EM1A2,fp2	...fp2 IS A2+S*(A4+S*A6)
    732 	FADD.S		#:3F000000,fp3	...fp3 IS A1+S*(A3+S*A5)
    733 
    734 	FMUL.X		fp1,fp2		...fp2 IS S*(A2+S*(A4+S*A6))
    735 	ORI.W		#$8000,d0	...signed/expo. of -2^(-M)
    736 	MOVE.W		d0,ONEBYSC(a6)	...OnebySc is -2^(-M)
    737 	clr.w		ONEBYSC+2(a6)
    738 	move.l		#$80000000,ONEBYSC+4(a6)
    739 	clr.l		ONEBYSC+8(a6)
    740 	FMUL.X		fp3,fp1		...fp1 IS S*(A1+S*(A3+S*A5))
    741 *					...fp3 released
    742 
    743 	FMUL.X		fp0,fp2		...fp2 IS R*S*(A2+S*(A4+S*A6))
    744 	FADD.X		fp1,fp0		...fp0 IS R+S*(A1+S*(A3+S*A5))
    745 *					...fp1 released
    746 
    747 	FADD.X		fp2,fp0		...fp0 IS EXP(R)-1
    748 *					...fp2 released
    749 	fmovem.x	(a7)+,fp2/fp3	...fp2 restored
    750 
    751 *--Step 5
    752 *--Compute 2^(J/64)*p
    753 
    754 	FMUL.X		(a1),fp0	...2^(J/64)*(Exp(R)-1)
    755 
    756 *--Step 6
    757 *--Step 6.1
    758 	MOVE.L		L_SCR1(a6),d0		...retrieve M
    759 	CMPI.L		#63,d0
    760 	BLE.B		MLE63
    761 *--Step 6.2	M >= 64
    762 	FMOVE.S		12(a1),fp1	...fp1 is t
    763 	FADD.X		ONEBYSC(a6),fp1	...fp1 is t+OnebySc
    764 	FADD.X		fp1,fp0		...p+(t+OnebySc), fp1 released
    765 	FADD.X		(a1),fp0	...T+(p+(t+OnebySc))
    766 	BRA.B		EM1SCALE
    767 MLE63:
    768 *--Step 6.3	M <= 63
    769 	CMPI.L		#-3,d0
    770 	BGE.B		MGEN3
    771 MLTN3:
    772 *--Step 6.4	M <= -4
    773 	FADD.S		12(a1),fp0	...p+t
    774 	FADD.X		(a1),fp0	...T+(p+t)
    775 	FADD.X		ONEBYSC(a6),fp0	...OnebySc + (T+(p+t))
    776 	BRA.B		EM1SCALE
    777 MGEN3:
    778 *--Step 6.5	-3 <= M <= 63
    779 	FMOVE.X		(a1)+,fp1	...fp1 is T
    780 	FADD.S		(a1),fp0	...fp0 is p+t
    781 	FADD.X		ONEBYSC(a6),fp1	...fp1 is T+OnebySc
    782 	FADD.X		fp1,fp0		...(T+OnebySc)+(p+t)
    783 
    784 EM1SCALE:
    785 *--Step 6.6
    786 	FMOVE.L		d1,FPCR
    787 	FMUL.X		SC(a6),fp0
    788 
    789 	bra		t_frcinx
    790 
    791 EM1SM:
    792 *--Step 7	|X| < 1/4.
    793 	CMPI.L		#$3FBE0000,d0	...2^(-65)
    794 	BGE.B		EM1POLY
    795 
    796 EM1TINY:
    797 *--Step 8	|X| < 2^(-65)
    798 	CMPI.L		#$00330000,d0	...2^(-16312)
    799 	BLT.B		EM12TINY
    800 *--Step 8.2
    801 	MOVE.L		#$80010000,SC(a6)	...SC is -2^(-16382)
    802 	move.l		#$80000000,SC+4(a6)
    803 	clr.l		SC+8(a6)
    804 	FMOVE.X		(a0),fp0
    805 	FMOVE.L		d1,FPCR
    806 	FADD.X		SC(a6),fp0
    807 
    808 	bra		t_frcinx
    809 
    810 EM12TINY:
    811 *--Step 8.3
    812 	FMOVE.X		(a0),fp0
    813 	FMUL.D		TWO140,fp0
    814 	MOVE.L		#$80010000,SC(a6)
    815 	move.l		#$80000000,SC+4(a6)
    816 	clr.l		SC+8(a6)
    817 	FADD.X		SC(a6),fp0
    818 	FMOVE.L		d1,FPCR
    819 	FMUL.D		TWON140,fp0
    820 
    821 	bra		t_frcinx
    822 
    823 EM1POLY:
    824 *--Step 9	exp(X)-1 by a simple polynomial
    825 	FMOVE.X		(a0),fp0	...fp0 is X
    826 	FMUL.X		fp0,fp0		...fp0 is S := X*X
    827 	fmovem.x	fp2/fp3,-(a7)	...save fp2
    828 	FMOVE.S		#:2F30CAA8,fp1	...fp1 is B12
    829 	FMUL.X		fp0,fp1		...fp1 is S*B12
    830 	FMOVE.S		#:310F8290,fp2	...fp2 is B11
    831 	FADD.S		#:32D73220,fp1	...fp1 is B10+S*B12
    832 
    833 	FMUL.X		fp0,fp2		...fp2 is S*B11
    834 	FMUL.X		fp0,fp1		...fp1 is S*(B10 + ...
    835 
    836 	FADD.S		#:3493F281,fp2	...fp2 is B9+S*...
    837 	FADD.D		EM1B8,fp1	...fp1 is B8+S*...
    838 
    839 	FMUL.X		fp0,fp2		...fp2 is S*(B9+...
    840 	FMUL.X		fp0,fp1		...fp1 is S*(B8+...
    841 
    842 	FADD.D		EM1B7,fp2	...fp2 is B7+S*...
    843 	FADD.D		EM1B6,fp1	...fp1 is B6+S*...
    844 
    845 	FMUL.X		fp0,fp2		...fp2 is S*(B7+...
    846 	FMUL.X		fp0,fp1		...fp1 is S*(B6+...
    847 
    848 	FADD.D		EM1B5,fp2	...fp2 is B5+S*...
    849 	FADD.D		EM1B4,fp1	...fp1 is B4+S*...
    850 
    851 	FMUL.X		fp0,fp2		...fp2 is S*(B5+...
    852 	FMUL.X		fp0,fp1		...fp1 is S*(B4+...
    853 
    854 	FADD.D		EM1B3,fp2	...fp2 is B3+S*...
    855 	FADD.X		EM1B2,fp1	...fp1 is B2+S*...
    856 
    857 	FMUL.X		fp0,fp2		...fp2 is S*(B3+...
    858 	FMUL.X		fp0,fp1		...fp1 is S*(B2+...
    859 
    860 	FMUL.X		fp0,fp2		...fp2 is S*S*(B3+...)
    861 	FMUL.X		(a0),fp1	...fp1 is X*S*(B2...
    862 
    863 	FMUL.S		#:3F000000,fp0	...fp0 is S*B1
    864 	FADD.X		fp2,fp1		...fp1 is Q
    865 *					...fp2 released
    866 
    867 	fmovem.x	(a7)+,fp2/fp3	...fp2 restored
    868 
    869 	FADD.X		fp1,fp0		...fp0 is S*B1+Q
    870 *					...fp1 released
    871 
    872 	FMOVE.L		d1,FPCR
    873 	FADD.X		(a0),fp0
    874 
    875 	bra		t_frcinx
    876 
    877 EM1BIG:
    878 *--Step 10	|X| > 70 log2
    879 	MOVE.L		(a0),d0
    880 	TST.L		d0
    881 	BGT.W		EXPC1
    882 *--Step 10.2
    883 	FMOVE.S		#:BF800000,fp0	...fp0 is -1
    884 	FMOVE.L		d1,FPCR
    885 	FADD.S		#:00800000,fp0	...-1 + 2^(-126)
    886 
    887 	bra		t_frcinx
    888 
    889 	end
    890