• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1|
2|	setox.sa 3.1 12/10/90
3|
4|	The entry point setox computes the exponential of a value.
5|	setoxd does the same except the input value is a denormalized
6|	number.	setoxm1 computes exp(X)-1, and setoxm1d computes
7|	exp(X)-1 for denormalized X.
8|
9|	INPUT
10|	-----
11|	Double-extended value in memory location pointed to by address
12|	register a0.
13|
14|	OUTPUT
15|	------
16|	exp(X) or exp(X)-1 returned in floating-point register fp0.
17|
18|	ACCURACY and MONOTONICITY
19|	-------------------------
20|	The returned result is within 0.85 ulps in 64 significant bit, i.e.
21|	within 0.5001 ulp to 53 bits if the result is subsequently rounded
22|	to double precision. The result is provably monotonic in double
23|	precision.
24|
25|	SPEED
26|	-----
27|	Two timings are measured, both in the copy-back mode. The
28|	first one is measured when the function is invoked the first time
29|	(so the instructions and data are not in cache), and the
30|	second one is measured when the function is reinvoked at the same
31|	input argument.
32|
33|	The program setox takes approximately 210/190 cycles for input
34|	argument X whose magnitude is less than 16380 log2, which
35|	is the usual situation.	For the less common arguments,
36|	depending on their values, the program may run faster or slower --
37|	but no worse than 10% slower even in the extreme cases.
38|
39|	The program setoxm1 takes approximately ??? / ??? cycles for input
40|	argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes
41|	approximately ??? / ??? cycles. For the less common arguments,
42|	depending on their values, the program may run faster or slower --
43|	but no worse than 10% slower even in the extreme cases.
44|
45|	ALGORITHM and IMPLEMENTATION NOTES
46|	----------------------------------
47|
48|	setoxd
49|	------
50|	Step 1.	Set ans := 1.0
51|
52|	Step 2.	Return	ans := ans + sign(X)*2^(-126). Exit.
53|	Notes:	This will always generate one exception -- inexact.
54|
55|
56|	setox
57|	-----
58|
59|	Step 1.	Filter out extreme cases of input argument.
60|		1.1	If |X| >= 2^(-65), go to Step 1.3.
61|		1.2	Go to Step 7.
62|		1.3	If |X| < 16380 log(2), go to Step 2.
63|		1.4	Go to Step 8.
64|	Notes:	The usual case should take the branches 1.1 -> 1.3 -> 2.
65|		 To avoid the use of floating-point comparisons, a
66|		 compact representation of |X| is used. This format is a
67|		 32-bit integer, the upper (more significant) 16 bits are
68|		 the sign and biased exponent field of |X|; the lower 16
69|		 bits are the 16 most significant fraction (including the
70|		 explicit bit) bits of |X|. Consequently, the comparisons
71|		 in Steps 1.1 and 1.3 can be performed by integer comparison.
72|		 Note also that the constant 16380 log(2) used in Step 1.3
73|		 is also in the compact form. Thus taking the branch
74|		 to Step 2 guarantees |X| < 16380 log(2). There is no harm
75|		 to have a small number of cases where |X| is less than,
76|		 but close to, 16380 log(2) and the branch to Step 9 is
77|		 taken.
78|
79|	Step 2.	Calculate N = round-to-nearest-int( X * 64/log2 ).
80|		2.1	Set AdjFlag := 0 (indicates the branch 1.3 -> 2 was taken)
81|		2.2	N := round-to-nearest-integer( X * 64/log2 ).
82|		2.3	Calculate	J = N mod 64; so J = 0,1,2,..., or 63.
83|		2.4	Calculate	M = (N - J)/64; so N = 64M + J.
84|		2.5	Calculate the address of the stored value of 2^(J/64).
85|		2.6	Create the value Scale = 2^M.
86|	Notes:	The calculation in 2.2 is really performed by
87|
88|			Z := X * constant
89|			N := round-to-nearest-integer(Z)
90|
91|		 where
92|
93|			constant := single-precision( 64/log 2 ).
94|
95|		 Using a single-precision constant avoids memory access.
96|		 Another effect of using a single-precision "constant" is
97|		 that the calculated value Z is
98|
99|			Z = X*(64/log2)*(1+eps), |eps| <= 2^(-24).
100|
101|		 This error has to be considered later in Steps 3 and 4.
102|
103|	Step 3.	Calculate X - N*log2/64.
104|		3.1	R := X + N*L1, where L1 := single-precision(-log2/64).
105|		3.2	R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
106|	Notes:	a) The way L1 and L2 are chosen ensures L1+L2 approximate
107|		 the value	-log2/64	to 88 bits of accuracy.
108|		 b) N*L1 is exact because N is no longer than 22 bits and
109|		 L1 is no longer than 24 bits.
110|		 c) The calculation X+N*L1 is also exact due to cancellation.
111|		 Thus, R is practically X+N(L1+L2) to full 64 bits.
112|		 d) It is important to estimate how large can |R| be after
113|		 Step 3.2.
114|
115|			N = rnd-to-int( X*64/log2 (1+eps) ), |eps|<=2^(-24)
116|			X*64/log2 (1+eps)	=	N + f,	|f| <= 0.5
117|			X*64/log2 - N	=	f - eps*X 64/log2
118|			X - N*log2/64	=	f*log2/64 - eps*X
119|
120|
121|		 Now |X| <= 16446 log2, thus
122|
123|			|X - N*log2/64| <= (0.5 + 16446/2^(18))*log2/64
124|					<= 0.57 log2/64.
125|		 This bound will be used in Step 4.
126|
127|	Step 4.	Approximate exp(R)-1 by a polynomial
128|			p = R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
129|	Notes:	a) In order to reduce memory access, the coefficients are
130|		 made as "short" as possible: A1 (which is 1/2), A4 and A5
131|		 are single precision; A2 and A3 are double precision.
132|		 b) Even with the restrictions above,
133|			|p - (exp(R)-1)| < 2^(-68.8) for all |R| <= 0.0062.
134|		 Note that 0.0062 is slightly bigger than 0.57 log2/64.
135|		 c) To fully utilize the pipeline, p is separated into
136|		 two independent pieces of roughly equal complexities
137|			p = [ R + R*S*(A2 + S*A4) ]	+
138|				[ S*(A1 + S*(A3 + S*A5)) ]
139|		 where S = R*R.
140|
141|	Step 5.	Compute 2^(J/64)*exp(R) = 2^(J/64)*(1+p) by
142|				ans := T + ( T*p + t)
143|		 where T and t are the stored values for 2^(J/64).
144|	Notes:	2^(J/64) is stored as T and t where T+t approximates
145|		 2^(J/64) to roughly 85 bits; T is in extended precision
146|		 and t is in single precision. Note also that T is rounded
147|		 to 62 bits so that the last two bits of T are zero. The
148|		 reason for such a special form is that T-1, T-2, and T-8
149|		 will all be exact --- a property that will give much
150|		 more accurate computation of the function EXPM1.
151|
152|	Step 6.	Reconstruction of exp(X)
153|			exp(X) = 2^M * 2^(J/64) * exp(R).
154|		6.1	If AdjFlag = 0, go to 6.3
155|		6.2	ans := ans * AdjScale
156|		6.3	Restore the user FPCR
157|		6.4	Return ans := ans * Scale. Exit.
158|	Notes:	If AdjFlag = 0, we have X = Mlog2 + Jlog2/64 + R,
159|		 |M| <= 16380, and Scale = 2^M. Moreover, exp(X) will
160|		 neither overflow nor underflow. If AdjFlag = 1, that
161|		 means that
162|			X = (M1+M)log2 + Jlog2/64 + R, |M1+M| >= 16380.
163|		 Hence, exp(X) may overflow or underflow or neither.
164|		 When that is the case, AdjScale = 2^(M1) where M1 is
165|		 approximately M. Thus 6.2 will never cause over/underflow.
166|		 Possible exception in 6.4 is overflow or underflow.
167|		 The inexact exception is not generated in 6.4. Although
168|		 one can argue that the inexact flag should always be
169|		 raised, to simulate that exception cost to much than the
170|		 flag is worth in practical uses.
171|
172|	Step 7.	Return 1 + X.
173|		7.1	ans := X
174|		7.2	Restore user FPCR.
175|		7.3	Return ans := 1 + ans. Exit
176|	Notes:	For non-zero X, the inexact exception will always be
177|		 raised by 7.3. That is the only exception raised by 7.3.
178|		 Note also that we use the FMOVEM instruction to move X
179|		 in Step 7.1 to avoid unnecessary trapping. (Although
180|		 the FMOVEM may not seem relevant since X is normalized,
181|		 the precaution will be useful in the library version of
182|		 this code where the separate entry for denormalized inputs
183|		 will be done away with.)
184|
185|	Step 8.	Handle exp(X) where |X| >= 16380log2.
186|		8.1	If |X| > 16480 log2, go to Step 9.
187|		(mimic 2.2 - 2.6)
188|		8.2	N := round-to-integer( X * 64/log2 )
189|		8.3	Calculate J = N mod 64, J = 0,1,...,63
190|		8.4	K := (N-J)/64, M1 := truncate(K/2), M = K-M1, AdjFlag := 1.
191|		8.5	Calculate the address of the stored value 2^(J/64).
192|		8.6	Create the values Scale = 2^M, AdjScale = 2^M1.
193|		8.7	Go to Step 3.
194|	Notes:	Refer to notes for 2.2 - 2.6.
195|
196|	Step 9.	Handle exp(X), |X| > 16480 log2.
197|		9.1	If X < 0, go to 9.3
198|		9.2	ans := Huge, go to 9.4
199|		9.3	ans := Tiny.
200|		9.4	Restore user FPCR.
201|		9.5	Return ans := ans * ans. Exit.
202|	Notes:	Exp(X) will surely overflow or underflow, depending on
203|		 X's sign. "Huge" and "Tiny" are respectively large/tiny
204|		 extended-precision numbers whose square over/underflow
205|		 with an inexact result. Thus, 9.5 always raises the
206|		 inexact together with either overflow or underflow.
207|
208|
209|	setoxm1d
210|	--------
211|
212|	Step 1.	Set ans := 0
213|
214|	Step 2.	Return	ans := X + ans. Exit.
215|	Notes:	This will return X with the appropriate rounding
216|		 precision prescribed by the user FPCR.
217|
218|	setoxm1
219|	-------
220|
221|	Step 1.	Check |X|
222|		1.1	If |X| >= 1/4, go to Step 1.3.
223|		1.2	Go to Step 7.
224|		1.3	If |X| < 70 log(2), go to Step 2.
225|		1.4	Go to Step 10.
226|	Notes:	The usual case should take the branches 1.1 -> 1.3 -> 2.
227|		 However, it is conceivable |X| can be small very often
228|		 because EXPM1 is intended to evaluate exp(X)-1 accurately
229|		 when |X| is small. For further details on the comparisons,
230|		 see the notes on Step 1 of setox.
231|
232|	Step 2.	Calculate N = round-to-nearest-int( X * 64/log2 ).
233|		2.1	N := round-to-nearest-integer( X * 64/log2 ).
234|		2.2	Calculate	J = N mod 64; so J = 0,1,2,..., or 63.
235|		2.3	Calculate	M = (N - J)/64; so N = 64M + J.
236|		2.4	Calculate the address of the stored value of 2^(J/64).
237|		2.5	Create the values Sc = 2^M and OnebySc := -2^(-M).
238|	Notes:	See the notes on Step 2 of setox.
239|
240|	Step 3.	Calculate X - N*log2/64.
241|		3.1	R := X + N*L1, where L1 := single-precision(-log2/64).
242|		3.2	R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
243|	Notes:	Applying the analysis of Step 3 of setox in this case
244|		 shows that |R| <= 0.0055 (note that |X| <= 70 log2 in
245|		 this case).
246|
247|	Step 4.	Approximate exp(R)-1 by a polynomial
248|			p = R+R*R*(A1+R*(A2+R*(A3+R*(A4+R*(A5+R*A6)))))
249|	Notes:	a) In order to reduce memory access, the coefficients are
250|		 made as "short" as possible: A1 (which is 1/2), A5 and A6
251|		 are single precision; A2, A3 and A4 are double precision.
252|		 b) Even with the restriction above,
253|			|p - (exp(R)-1)| <	|R| * 2^(-72.7)
254|		 for all |R| <= 0.0055.
255|		 c) To fully utilize the pipeline, p is separated into
256|		 two independent pieces of roughly equal complexity
257|			p = [ R*S*(A2 + S*(A4 + S*A6)) ]	+
258|				[ R + S*(A1 + S*(A3 + S*A5)) ]
259|		 where S = R*R.
260|
261|	Step 5.	Compute 2^(J/64)*p by
262|				p := T*p
263|		 where T and t are the stored values for 2^(J/64).
264|	Notes:	2^(J/64) is stored as T and t where T+t approximates
265|		 2^(J/64) to roughly 85 bits; T is in extended precision
266|		 and t is in single precision. Note also that T is rounded
267|		 to 62 bits so that the last two bits of T are zero. The
268|		 reason for such a special form is that T-1, T-2, and T-8
269|		 will all be exact --- a property that will be exploited
270|		 in Step 6 below. The total relative error in p is no
271|		 bigger than 2^(-67.7) compared to the final result.
272|
273|	Step 6.	Reconstruction of exp(X)-1
274|			exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ).
275|		6.1	If M <= 63, go to Step 6.3.
276|		6.2	ans := T + (p + (t + OnebySc)). Go to 6.6
277|		6.3	If M >= -3, go to 6.5.
278|		6.4	ans := (T + (p + t)) + OnebySc. Go to 6.6
279|		6.5	ans := (T + OnebySc) + (p + t).
280|		6.6	Restore user FPCR.
281|		6.7	Return ans := Sc * ans. Exit.
282|	Notes:	The various arrangements of the expressions give accurate
283|		 evaluations.
284|
285|	Step 7.	exp(X)-1 for |X| < 1/4.
286|		7.1	If |X| >= 2^(-65), go to Step 9.
287|		7.2	Go to Step 8.
288|
289|	Step 8.	Calculate exp(X)-1, |X| < 2^(-65).
290|		8.1	If |X| < 2^(-16312), goto 8.3
291|		8.2	Restore FPCR; return ans := X - 2^(-16382). Exit.
292|		8.3	X := X * 2^(140).
293|		8.4	Restore FPCR; ans := ans - 2^(-16382).
294|		 Return ans := ans*2^(140). Exit
295|	Notes:	The idea is to return "X - tiny" under the user
296|		 precision and rounding modes. To avoid unnecessary
297|		 inefficiency, we stay away from denormalized numbers the
298|		 best we can. For |X| >= 2^(-16312), the straightforward
299|		 8.2 generates the inexact exception as the case warrants.
300|
301|	Step 9.	Calculate exp(X)-1, |X| < 1/4, by a polynomial
302|			p = X + X*X*(B1 + X*(B2 + ... + X*B12))
303|	Notes:	a) In order to reduce memory access, the coefficients are
304|		 made as "short" as possible: B1 (which is 1/2), B9 to B12
305|		 are single precision; B3 to B8 are double precision; and
306|		 B2 is double extended.
307|		 b) Even with the restriction above,
308|			|p - (exp(X)-1)| < |X| 2^(-70.6)
309|		 for all |X| <= 0.251.
310|		 Note that 0.251 is slightly bigger than 1/4.
311|		 c) To fully preserve accuracy, the polynomial is computed
312|		 as	X + ( S*B1 +	Q ) where S = X*X and
313|			Q	=	X*S*(B2 + X*(B3 + ... + X*B12))
314|		 d) To fully utilize the pipeline, Q is separated into
315|		 two independent pieces of roughly equal complexity
316|			Q = [ X*S*(B2 + S*(B4 + ... + S*B12)) ] +
317|				[ S*S*(B3 + S*(B5 + ... + S*B11)) ]
318|
319|	Step 10.	Calculate exp(X)-1 for |X| >= 70 log 2.
320|		10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all practical
321|		 purposes. Therefore, go to Step 1 of setox.
322|		10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical purposes.
323|		 ans := -1
324|		 Restore user FPCR
325|		 Return ans := ans + 2^(-126). Exit.
326|	Notes:	10.2 will always create an inexact and return -1 + tiny
327|		 in the user rounding precision and mode.
328|
329|
330
331|		Copyright (C) Motorola, Inc. 1990
332|			All Rights Reserved
333|
334|       For details on the license for this file, please see the
335|       file, README, in this same directory.
336
337|setox	idnt	2,1 | Motorola 040 Floating Point Software Package
338
339	|section	8
340
341#include "fpsp.h"
342
343L2:	.long	0x3FDC0000,0x82E30865,0x4361C4C6,0x00000000
344
345EXPA3:	.long	0x3FA55555,0x55554431
346EXPA2:	.long	0x3FC55555,0x55554018
347
348HUGE:	.long	0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
349TINY:	.long	0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
350
351EM1A4:	.long	0x3F811111,0x11174385
352EM1A3:	.long	0x3FA55555,0x55554F5A
353
354EM1A2:	.long	0x3FC55555,0x55555555,0x00000000,0x00000000
355
356EM1B8:	.long	0x3EC71DE3,0xA5774682
357EM1B7:	.long	0x3EFA01A0,0x19D7CB68
358
359EM1B6:	.long	0x3F2A01A0,0x1A019DF3
360EM1B5:	.long	0x3F56C16C,0x16C170E2
361
362EM1B4:	.long	0x3F811111,0x11111111
363EM1B3:	.long	0x3FA55555,0x55555555
364
365EM1B2:	.long	0x3FFC0000,0xAAAAAAAA,0xAAAAAAAB
366	.long	0x00000000
367
368TWO140:	.long	0x48B00000,0x00000000
369TWON140:	.long	0x37300000,0x00000000
370
371EXPTBL:
372	.long	0x3FFF0000,0x80000000,0x00000000,0x00000000
373	.long	0x3FFF0000,0x8164D1F3,0xBC030774,0x9F841A9B
374	.long	0x3FFF0000,0x82CD8698,0xAC2BA1D8,0x9FC1D5B9
375	.long	0x3FFF0000,0x843A28C3,0xACDE4048,0xA0728369
376	.long	0x3FFF0000,0x85AAC367,0xCC487B14,0x1FC5C95C
377	.long	0x3FFF0000,0x871F6196,0x9E8D1010,0x1EE85C9F
378	.long	0x3FFF0000,0x88980E80,0x92DA8528,0x9FA20729
379	.long	0x3FFF0000,0x8A14D575,0x496EFD9C,0xA07BF9AF
380	.long	0x3FFF0000,0x8B95C1E3,0xEA8BD6E8,0xA0020DCF
381	.long	0x3FFF0000,0x8D1ADF5B,0x7E5BA9E4,0x205A63DA
382	.long	0x3FFF0000,0x8EA4398B,0x45CD53C0,0x1EB70051
383	.long	0x3FFF0000,0x9031DC43,0x1466B1DC,0x1F6EB029
384	.long	0x3FFF0000,0x91C3D373,0xAB11C338,0xA0781494
385	.long	0x3FFF0000,0x935A2B2F,0x13E6E92C,0x9EB319B0
386	.long	0x3FFF0000,0x94F4EFA8,0xFEF70960,0x2017457D
387	.long	0x3FFF0000,0x96942D37,0x20185A00,0x1F11D537
388	.long	0x3FFF0000,0x9837F051,0x8DB8A970,0x9FB952DD
389	.long	0x3FFF0000,0x99E04593,0x20B7FA64,0x1FE43087
390	.long	0x3FFF0000,0x9B8D39B9,0xD54E5538,0x1FA2A818
391	.long	0x3FFF0000,0x9D3ED9A7,0x2CFFB750,0x1FDE494D
392	.long	0x3FFF0000,0x9EF53260,0x91A111AC,0x20504890
393	.long	0x3FFF0000,0xA0B0510F,0xB9714FC4,0xA073691C
394	.long	0x3FFF0000,0xA2704303,0x0C496818,0x1F9B7A05
395	.long	0x3FFF0000,0xA43515AE,0x09E680A0,0xA0797126
396	.long	0x3FFF0000,0xA5FED6A9,0xB15138EC,0xA071A140
397	.long	0x3FFF0000,0xA7CD93B4,0xE9653568,0x204F62DA
398	.long	0x3FFF0000,0xA9A15AB4,0xEA7C0EF8,0x1F283C4A
399	.long	0x3FFF0000,0xAB7A39B5,0xA93ED338,0x9F9A7FDC
400	.long	0x3FFF0000,0xAD583EEA,0x42A14AC8,0xA05B3FAC
401	.long	0x3FFF0000,0xAF3B78AD,0x690A4374,0x1FDF2610
402	.long	0x3FFF0000,0xB123F581,0xD2AC2590,0x9F705F90
403	.long	0x3FFF0000,0xB311C412,0xA9112488,0x201F678A
404	.long	0x3FFF0000,0xB504F333,0xF9DE6484,0x1F32FB13
405	.long	0x3FFF0000,0xB6FD91E3,0x28D17790,0x20038B30
406	.long	0x3FFF0000,0xB8FBAF47,0x62FB9EE8,0x200DC3CC
407	.long	0x3FFF0000,0xBAFF5AB2,0x133E45FC,0x9F8B2AE6
408	.long	0x3FFF0000,0xBD08A39F,0x580C36C0,0xA02BBF70
409	.long	0x3FFF0000,0xBF1799B6,0x7A731084,0xA00BF518
410	.long	0x3FFF0000,0xC12C4CCA,0x66709458,0xA041DD41
411	.long	0x3FFF0000,0xC346CCDA,0x24976408,0x9FDF137B
412	.long	0x3FFF0000,0xC5672A11,0x5506DADC,0x201F1568
413	.long	0x3FFF0000,0xC78D74C8,0xABB9B15C,0x1FC13A2E
414	.long	0x3FFF0000,0xC9B9BD86,0x6E2F27A4,0xA03F8F03
415	.long	0x3FFF0000,0xCBEC14FE,0xF2727C5C,0x1FF4907D
416	.long	0x3FFF0000,0xCE248C15,0x1F8480E4,0x9E6E53E4
417	.long	0x3FFF0000,0xD06333DA,0xEF2B2594,0x1FD6D45C
418	.long	0x3FFF0000,0xD2A81D91,0xF12AE45C,0xA076EDB9
419	.long	0x3FFF0000,0xD4F35AAB,0xCFEDFA20,0x9FA6DE21
420	.long	0x3FFF0000,0xD744FCCA,0xD69D6AF4,0x1EE69A2F
421	.long	0x3FFF0000,0xD99D15C2,0x78AFD7B4,0x207F439F
422	.long	0x3FFF0000,0xDBFBB797,0xDAF23754,0x201EC207
423	.long	0x3FFF0000,0xDE60F482,0x5E0E9124,0x9E8BE175
424	.long	0x3FFF0000,0xE0CCDEEC,0x2A94E110,0x20032C4B
425	.long	0x3FFF0000,0xE33F8972,0xBE8A5A50,0x2004DFF5
426	.long	0x3FFF0000,0xE5B906E7,0x7C8348A8,0x1E72F47A
427	.long	0x3FFF0000,0xE8396A50,0x3C4BDC68,0x1F722F22
428	.long	0x3FFF0000,0xEAC0C6E7,0xDD243930,0xA017E945
429	.long	0x3FFF0000,0xED4F301E,0xD9942B84,0x1F401A5B
430	.long	0x3FFF0000,0xEFE4B99B,0xDCDAF5CC,0x9FB9A9E3
431	.long	0x3FFF0000,0xF281773C,0x59FFB138,0x20744C05
432	.long	0x3FFF0000,0xF5257D15,0x2486CC2C,0x1F773A19
433	.long	0x3FFF0000,0xF7D0DF73,0x0AD13BB8,0x1FFE90D5
434	.long	0x3FFF0000,0xFA83B2DB,0x722A033C,0xA041ED22
435	.long	0x3FFF0000,0xFD3E0C0C,0xF486C174,0x1F853F3A
436
437	.set	ADJFLAG,L_SCR2
438	.set	SCALE,FP_SCR1
439	.set	ADJSCALE,FP_SCR2
440	.set	SC,FP_SCR3
441	.set	ONEBYSC,FP_SCR4
442
443	| xref	t_frcinx
444	|xref	t_extdnrm
445	|xref	t_unfl
446	|xref	t_ovfl
447
448	.global	setoxd
449setoxd:
450|--entry point for EXP(X), X is denormalized
451	movel		(%a0),%d0
452	andil		#0x80000000,%d0
453	oril		#0x00800000,%d0		| ...sign(X)*2^(-126)
454	movel		%d0,-(%sp)
455	fmoves		#0x3F800000,%fp0
456	fmovel		%d1,%fpcr
457	fadds		(%sp)+,%fp0
458	bra		t_frcinx
459
460	.global	setox
461setox:
462|--entry point for EXP(X), here X is finite, non-zero, and not NaN's
463
464|--Step 1.
465	movel		(%a0),%d0	 | ...load part of input X
466	andil		#0x7FFF0000,%d0	| ...biased expo. of X
467	cmpil		#0x3FBE0000,%d0	| ...2^(-65)
468	bges		EXPC1		| ...normal case
469	bra		EXPSM
470
471EXPC1:
472|--The case |X| >= 2^(-65)
473	movew		4(%a0),%d0	| ...expo. and partial sig. of |X|
474	cmpil		#0x400CB167,%d0	| ...16380 log2 trunc. 16 bits
475	blts		EXPMAIN	 | ...normal case
476	bra		EXPBIG
477
478EXPMAIN:
479|--Step 2.
480|--This is the normal branch:	2^(-65) <= |X| < 16380 log2.
481	fmovex		(%a0),%fp0	| ...load input from (a0)
482
483	fmovex		%fp0,%fp1
484	fmuls		#0x42B8AA3B,%fp0	| ...64/log2 * X
485	fmovemx	%fp2-%fp2/%fp3,-(%a7)		| ...save fp2
486	movel		#0,ADJFLAG(%a6)
487	fmovel		%fp0,%d0		| ...N = int( X * 64/log2 )
488	lea		EXPTBL,%a1
489	fmovel		%d0,%fp0		| ...convert to floating-format
490
491	movel		%d0,L_SCR1(%a6)	| ...save N temporarily
492	andil		#0x3F,%d0		| ...D0 is J = N mod 64
493	lsll		#4,%d0
494	addal		%d0,%a1		| ...address of 2^(J/64)
495	movel		L_SCR1(%a6),%d0
496	asrl		#6,%d0		| ...D0 is M
497	addiw		#0x3FFF,%d0	| ...biased expo. of 2^(M)
498	movew		L2,L_SCR1(%a6)	| ...prefetch L2, no need in CB
499
500EXPCONT1:
501|--Step 3.
502|--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
503|--a0 points to 2^(J/64), D0 is biased expo. of 2^(M)
504	fmovex		%fp0,%fp2
505	fmuls		#0xBC317218,%fp0	| ...N * L1, L1 = lead(-log2/64)
506	fmulx		L2,%fp2		| ...N * L2, L1+L2 = -log2/64
507	faddx		%fp1,%fp0		| ...X + N*L1
508	faddx		%fp2,%fp0		| ...fp0 is R, reduced arg.
509|	MOVE.W		#$3FA5,EXPA3	...load EXPA3 in cache
510
511|--Step 4.
512|--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
513|-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
514|--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
515|--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))]
516
517	fmovex		%fp0,%fp1
518	fmulx		%fp1,%fp1		| ...fp1 IS S = R*R
519
520	fmoves		#0x3AB60B70,%fp2	| ...fp2 IS A5
521|	MOVE.W		#0,2(%a1)	...load 2^(J/64) in cache
522
523	fmulx		%fp1,%fp2		| ...fp2 IS S*A5
524	fmovex		%fp1,%fp3
525	fmuls		#0x3C088895,%fp3	| ...fp3 IS S*A4
526
527	faddd		EXPA3,%fp2	| ...fp2 IS A3+S*A5
528	faddd		EXPA2,%fp3	| ...fp3 IS A2+S*A4
529
530	fmulx		%fp1,%fp2		| ...fp2 IS S*(A3+S*A5)
531	movew		%d0,SCALE(%a6)	| ...SCALE is 2^(M) in extended
532	clrw		SCALE+2(%a6)
533	movel		#0x80000000,SCALE+4(%a6)
534	clrl		SCALE+8(%a6)
535
536	fmulx		%fp1,%fp3		| ...fp3 IS S*(A2+S*A4)
537
538	fadds		#0x3F000000,%fp2	| ...fp2 IS A1+S*(A3+S*A5)
539	fmulx		%fp0,%fp3		| ...fp3 IS R*S*(A2+S*A4)
540
541	fmulx		%fp1,%fp2		| ...fp2 IS S*(A1+S*(A3+S*A5))
542	faddx		%fp3,%fp0		| ...fp0 IS R+R*S*(A2+S*A4),
543|					...fp3 released
544
545	fmovex		(%a1)+,%fp1	| ...fp1 is lead. pt. of 2^(J/64)
546	faddx		%fp2,%fp0		| ...fp0 is EXP(R) - 1
547|					...fp2 released
548
549|--Step 5
550|--final reconstruction process
551|--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) )
552
553	fmulx		%fp1,%fp0		| ...2^(J/64)*(Exp(R)-1)
554	fmovemx	(%a7)+,%fp2-%fp2/%fp3	| ...fp2 restored
555	fadds		(%a1),%fp0	| ...accurate 2^(J/64)
556
557	faddx		%fp1,%fp0		| ...2^(J/64) + 2^(J/64)*...
558	movel		ADJFLAG(%a6),%d0
559
560|--Step 6
561	tstl		%d0
562	beqs		NORMAL
563ADJUST:
564	fmulx		ADJSCALE(%a6),%fp0
565NORMAL:
566	fmovel		%d1,%FPCR		| ...restore user FPCR
567	fmulx		SCALE(%a6),%fp0	| ...multiply 2^(M)
568	bra		t_frcinx
569
570EXPSM:
571|--Step 7
572	fmovemx	(%a0),%fp0-%fp0	| ...in case X is denormalized
573	fmovel		%d1,%FPCR
574	fadds		#0x3F800000,%fp0	| ...1+X in user mode
575	bra		t_frcinx
576
577EXPBIG:
578|--Step 8
579	cmpil		#0x400CB27C,%d0	| ...16480 log2
580	bgts		EXP2BIG
581|--Steps 8.2 -- 8.6
582	fmovex		(%a0),%fp0	| ...load input from (a0)
583
584	fmovex		%fp0,%fp1
585	fmuls		#0x42B8AA3B,%fp0	| ...64/log2 * X
586	fmovemx	 %fp2-%fp2/%fp3,-(%a7)		| ...save fp2
587	movel		#1,ADJFLAG(%a6)
588	fmovel		%fp0,%d0		| ...N = int( X * 64/log2 )
589	lea		EXPTBL,%a1
590	fmovel		%d0,%fp0		| ...convert to floating-format
591	movel		%d0,L_SCR1(%a6)			| ...save N temporarily
592	andil		#0x3F,%d0		 | ...D0 is J = N mod 64
593	lsll		#4,%d0
594	addal		%d0,%a1			| ...address of 2^(J/64)
595	movel		L_SCR1(%a6),%d0
596	asrl		#6,%d0			| ...D0 is K
597	movel		%d0,L_SCR1(%a6)			| ...save K temporarily
598	asrl		#1,%d0			| ...D0 is M1
599	subl		%d0,L_SCR1(%a6)			| ...a1 is M
600	addiw		#0x3FFF,%d0		| ...biased expo. of 2^(M1)
601	movew		%d0,ADJSCALE(%a6)		| ...ADJSCALE := 2^(M1)
602	clrw		ADJSCALE+2(%a6)
603	movel		#0x80000000,ADJSCALE+4(%a6)
604	clrl		ADJSCALE+8(%a6)
605	movel		L_SCR1(%a6),%d0			| ...D0 is M
606	addiw		#0x3FFF,%d0		| ...biased expo. of 2^(M)
607	bra		EXPCONT1		| ...go back to Step 3
608
609EXP2BIG:
610|--Step 9
611	fmovel		%d1,%FPCR
612	movel		(%a0),%d0
613	bclrb		#sign_bit,(%a0)		| ...setox always returns positive
614	cmpil		#0,%d0
615	blt		t_unfl
616	bra		t_ovfl
617
618	.global	setoxm1d
619setoxm1d:
620|--entry point for EXPM1(X), here X is denormalized
621|--Step 0.
622	bra		t_extdnrm
623
624
625	.global	setoxm1
626setoxm1:
627|--entry point for EXPM1(X), here X is finite, non-zero, non-NaN
628
629|--Step 1.
630|--Step 1.1
631	movel		(%a0),%d0	 | ...load part of input X
632	andil		#0x7FFF0000,%d0	| ...biased expo. of X
633	cmpil		#0x3FFD0000,%d0	| ...1/4
634	bges		EM1CON1	 | ...|X| >= 1/4
635	bra		EM1SM
636
637EM1CON1:
638|--Step 1.3
639|--The case |X| >= 1/4
640	movew		4(%a0),%d0	| ...expo. and partial sig. of |X|
641	cmpil		#0x4004C215,%d0	| ...70log2 rounded up to 16 bits
642	bles		EM1MAIN	 | ...1/4 <= |X| <= 70log2
643	bra		EM1BIG
644
645EM1MAIN:
646|--Step 2.
647|--This is the case:	1/4 <= |X| <= 70 log2.
648	fmovex		(%a0),%fp0	| ...load input from (a0)
649
650	fmovex		%fp0,%fp1
651	fmuls		#0x42B8AA3B,%fp0	| ...64/log2 * X
652	fmovemx	%fp2-%fp2/%fp3,-(%a7)		| ...save fp2
653|	MOVE.W		#$3F81,EM1A4		...prefetch in CB mode
654	fmovel		%fp0,%d0		| ...N = int( X * 64/log2 )
655	lea		EXPTBL,%a1
656	fmovel		%d0,%fp0		| ...convert to floating-format
657
658	movel		%d0,L_SCR1(%a6)			| ...save N temporarily
659	andil		#0x3F,%d0		 | ...D0 is J = N mod 64
660	lsll		#4,%d0
661	addal		%d0,%a1			| ...address of 2^(J/64)
662	movel		L_SCR1(%a6),%d0
663	asrl		#6,%d0			| ...D0 is M
664	movel		%d0,L_SCR1(%a6)			| ...save a copy of M
665|	MOVE.W		#$3FDC,L2		...prefetch L2 in CB mode
666
667|--Step 3.
668|--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
669|--a0 points to 2^(J/64), D0 and a1 both contain M
670	fmovex		%fp0,%fp2
671	fmuls		#0xBC317218,%fp0	| ...N * L1, L1 = lead(-log2/64)
672	fmulx		L2,%fp2		| ...N * L2, L1+L2 = -log2/64
673	faddx		%fp1,%fp0	 | ...X + N*L1
674	faddx		%fp2,%fp0	 | ...fp0 is R, reduced arg.
675|	MOVE.W		#$3FC5,EM1A2		...load EM1A2 in cache
676	addiw		#0x3FFF,%d0		| ...D0 is biased expo. of 2^M
677
678|--Step 4.
679|--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
680|-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6)))))
681|--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
682|--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))]
683
684	fmovex		%fp0,%fp1
685	fmulx		%fp1,%fp1		| ...fp1 IS S = R*R
686
687	fmoves		#0x3950097B,%fp2	| ...fp2 IS a6
688|	MOVE.W		#0,2(%a1)	...load 2^(J/64) in cache
689
690	fmulx		%fp1,%fp2		| ...fp2 IS S*A6
691	fmovex		%fp1,%fp3
692	fmuls		#0x3AB60B6A,%fp3	| ...fp3 IS S*A5
693
694	faddd		EM1A4,%fp2	| ...fp2 IS A4+S*A6
695	faddd		EM1A3,%fp3	| ...fp3 IS A3+S*A5
696	movew		%d0,SC(%a6)		| ...SC is 2^(M) in extended
697	clrw		SC+2(%a6)
698	movel		#0x80000000,SC+4(%a6)
699	clrl		SC+8(%a6)
700
701	fmulx		%fp1,%fp2		| ...fp2 IS S*(A4+S*A6)
702	movel		L_SCR1(%a6),%d0		| ...D0 is	M
703	negw		%d0		| ...D0 is -M
704	fmulx		%fp1,%fp3		| ...fp3 IS S*(A3+S*A5)
705	addiw		#0x3FFF,%d0	| ...biased expo. of 2^(-M)
706	faddd		EM1A2,%fp2	| ...fp2 IS A2+S*(A4+S*A6)
707	fadds		#0x3F000000,%fp3	| ...fp3 IS A1+S*(A3+S*A5)
708
709	fmulx		%fp1,%fp2		| ...fp2 IS S*(A2+S*(A4+S*A6))
710	oriw		#0x8000,%d0	| ...signed/expo. of -2^(-M)
711	movew		%d0,ONEBYSC(%a6)	| ...OnebySc is -2^(-M)
712	clrw		ONEBYSC+2(%a6)
713	movel		#0x80000000,ONEBYSC+4(%a6)
714	clrl		ONEBYSC+8(%a6)
715	fmulx		%fp3,%fp1		| ...fp1 IS S*(A1+S*(A3+S*A5))
716|					...fp3 released
717
718	fmulx		%fp0,%fp2		| ...fp2 IS R*S*(A2+S*(A4+S*A6))
719	faddx		%fp1,%fp0		| ...fp0 IS R+S*(A1+S*(A3+S*A5))
720|					...fp1 released
721
722	faddx		%fp2,%fp0		| ...fp0 IS EXP(R)-1
723|					...fp2 released
724	fmovemx	(%a7)+,%fp2-%fp2/%fp3	| ...fp2 restored
725
726|--Step 5
727|--Compute 2^(J/64)*p
728
729	fmulx		(%a1),%fp0	| ...2^(J/64)*(Exp(R)-1)
730
731|--Step 6
732|--Step 6.1
733	movel		L_SCR1(%a6),%d0		| ...retrieve M
734	cmpil		#63,%d0
735	bles		MLE63
736|--Step 6.2	M >= 64
737	fmoves		12(%a1),%fp1	| ...fp1 is t
738	faddx		ONEBYSC(%a6),%fp1	| ...fp1 is t+OnebySc
739	faddx		%fp1,%fp0		| ...p+(t+OnebySc), fp1 released
740	faddx		(%a1),%fp0	| ...T+(p+(t+OnebySc))
741	bras		EM1SCALE
742MLE63:
743|--Step 6.3	M <= 63
744	cmpil		#-3,%d0
745	bges		MGEN3
746MLTN3:
747|--Step 6.4	M <= -4
748	fadds		12(%a1),%fp0	| ...p+t
749	faddx		(%a1),%fp0	| ...T+(p+t)
750	faddx		ONEBYSC(%a6),%fp0	| ...OnebySc + (T+(p+t))
751	bras		EM1SCALE
752MGEN3:
753|--Step 6.5	-3 <= M <= 63
754	fmovex		(%a1)+,%fp1	| ...fp1 is T
755	fadds		(%a1),%fp0	| ...fp0 is p+t
756	faddx		ONEBYSC(%a6),%fp1	| ...fp1 is T+OnebySc
757	faddx		%fp1,%fp0		| ...(T+OnebySc)+(p+t)
758
759EM1SCALE:
760|--Step 6.6
761	fmovel		%d1,%FPCR
762	fmulx		SC(%a6),%fp0
763
764	bra		t_frcinx
765
766EM1SM:
767|--Step 7	|X| < 1/4.
768	cmpil		#0x3FBE0000,%d0	| ...2^(-65)
769	bges		EM1POLY
770
771EM1TINY:
772|--Step 8	|X| < 2^(-65)
773	cmpil		#0x00330000,%d0	| ...2^(-16312)
774	blts		EM12TINY
775|--Step 8.2
776	movel		#0x80010000,SC(%a6)	| ...SC is -2^(-16382)
777	movel		#0x80000000,SC+4(%a6)
778	clrl		SC+8(%a6)
779	fmovex		(%a0),%fp0
780	fmovel		%d1,%FPCR
781	faddx		SC(%a6),%fp0
782
783	bra		t_frcinx
784
785EM12TINY:
786|--Step 8.3
787	fmovex		(%a0),%fp0
788	fmuld		TWO140,%fp0
789	movel		#0x80010000,SC(%a6)
790	movel		#0x80000000,SC+4(%a6)
791	clrl		SC+8(%a6)
792	faddx		SC(%a6),%fp0
793	fmovel		%d1,%FPCR
794	fmuld		TWON140,%fp0
795
796	bra		t_frcinx
797
798EM1POLY:
799|--Step 9	exp(X)-1 by a simple polynomial
800	fmovex		(%a0),%fp0	| ...fp0 is X
801	fmulx		%fp0,%fp0		| ...fp0 is S := X*X
802	fmovemx	%fp2-%fp2/%fp3,-(%a7)	| ...save fp2
803	fmoves		#0x2F30CAA8,%fp1	| ...fp1 is B12
804	fmulx		%fp0,%fp1		| ...fp1 is S*B12
805	fmoves		#0x310F8290,%fp2	| ...fp2 is B11
806	fadds		#0x32D73220,%fp1	| ...fp1 is B10+S*B12
807
808	fmulx		%fp0,%fp2		| ...fp2 is S*B11
809	fmulx		%fp0,%fp1		| ...fp1 is S*(B10 + ...
810
811	fadds		#0x3493F281,%fp2	| ...fp2 is B9+S*...
812	faddd		EM1B8,%fp1	| ...fp1 is B8+S*...
813
814	fmulx		%fp0,%fp2		| ...fp2 is S*(B9+...
815	fmulx		%fp0,%fp1		| ...fp1 is S*(B8+...
816
817	faddd		EM1B7,%fp2	| ...fp2 is B7+S*...
818	faddd		EM1B6,%fp1	| ...fp1 is B6+S*...
819
820	fmulx		%fp0,%fp2		| ...fp2 is S*(B7+...
821	fmulx		%fp0,%fp1		| ...fp1 is S*(B6+...
822
823	faddd		EM1B5,%fp2	| ...fp2 is B5+S*...
824	faddd		EM1B4,%fp1	| ...fp1 is B4+S*...
825
826	fmulx		%fp0,%fp2		| ...fp2 is S*(B5+...
827	fmulx		%fp0,%fp1		| ...fp1 is S*(B4+...
828
829	faddd		EM1B3,%fp2	| ...fp2 is B3+S*...
830	faddx		EM1B2,%fp1	| ...fp1 is B2+S*...
831
832	fmulx		%fp0,%fp2		| ...fp2 is S*(B3+...
833	fmulx		%fp0,%fp1		| ...fp1 is S*(B2+...
834
835	fmulx		%fp0,%fp2		| ...fp2 is S*S*(B3+...)
836	fmulx		(%a0),%fp1	| ...fp1 is X*S*(B2...
837
838	fmuls		#0x3F000000,%fp0	| ...fp0 is S*B1
839	faddx		%fp2,%fp1		| ...fp1 is Q
840|					...fp2 released
841
842	fmovemx	(%a7)+,%fp2-%fp2/%fp3	| ...fp2 restored
843
844	faddx		%fp1,%fp0		| ...fp0 is S*B1+Q
845|					...fp1 released
846
847	fmovel		%d1,%FPCR
848	faddx		(%a0),%fp0
849
850	bra		t_frcinx
851
852EM1BIG:
853|--Step 10	|X| > 70 log2
854	movel		(%a0),%d0
855	cmpil		#0,%d0
856	bgt		EXPC1
857|--Step 10.2
858	fmoves		#0xBF800000,%fp0	| ...fp0 is -1
859	fmovel		%d1,%FPCR
860	fadds		#0x00800000,%fp0	| ...-1 + 2^(-126)
861
862	bra		t_frcinx
863
864	|end
865