Commit | Line | Data |
---|---|---|
1da177e4 LT |
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 | | | |
646819e8 | 39 | | The program setoxm1 takes approximately ??? / ??? cycles for input |
1da177e4 | 40 | | argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes |
646819e8 | 41 | | approximately ??? / ??? cycles. For the less common arguments, |
1da177e4 LT |
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 | | | |
e00d82d0 MW |
334 | | For details on the license for this file, please see the |
335 | | file, README, in this same directory. | |
1da177e4 LT |
336 | |
337 | |setox idnt 2,1 | Motorola 040 Floating Point Software Package | |
338 | ||
339 | |section 8 | |
340 | ||
341 | #include "fpsp.h" | |
342 | ||
343 | L2: .long 0x3FDC0000,0x82E30865,0x4361C4C6,0x00000000 | |
344 | ||
345 | EXPA3: .long 0x3FA55555,0x55554431 | |
346 | EXPA2: .long 0x3FC55555,0x55554018 | |
347 | ||
348 | HUGE: .long 0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000 | |
349 | TINY: .long 0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000 | |
350 | ||
351 | EM1A4: .long 0x3F811111,0x11174385 | |
352 | EM1A3: .long 0x3FA55555,0x55554F5A | |
353 | ||
354 | EM1A2: .long 0x3FC55555,0x55555555,0x00000000,0x00000000 | |
355 | ||
356 | EM1B8: .long 0x3EC71DE3,0xA5774682 | |
357 | EM1B7: .long 0x3EFA01A0,0x19D7CB68 | |
358 | ||
359 | EM1B6: .long 0x3F2A01A0,0x1A019DF3 | |
360 | EM1B5: .long 0x3F56C16C,0x16C170E2 | |
361 | ||
362 | EM1B4: .long 0x3F811111,0x11111111 | |
363 | EM1B3: .long 0x3FA55555,0x55555555 | |
364 | ||
365 | EM1B2: .long 0x3FFC0000,0xAAAAAAAA,0xAAAAAAAB | |
366 | .long 0x00000000 | |
367 | ||
368 | TWO140: .long 0x48B00000,0x00000000 | |
369 | TWON140: .long 0x37300000,0x00000000 | |
370 | ||
371 | EXPTBL: | |
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 | |
449 | setoxd: | |
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 | |
461 | setox: | |
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 | ||
471 | EXPC1: | |
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 | ||
478 | EXPMAIN: | |
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 | ||
500 | EXPCONT1: | |
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 | |
563 | ADJUST: | |
564 | fmulx ADJSCALE(%a6),%fp0 | |
565 | NORMAL: | |
566 | fmovel %d1,%FPCR | ...restore user FPCR | |
567 | fmulx SCALE(%a6),%fp0 | ...multiply 2^(M) | |
568 | bra t_frcinx | |
569 | ||
570 | EXPSM: | |
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 | ||
577 | EXPBIG: | |
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 | ||
609 | EXP2BIG: | |
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 | |
619 | setoxm1d: | |
620 | |--entry point for EXPM1(X), here X is denormalized | |
621 | |--Step 0. | |
622 | bra t_extdnrm | |
623 | ||
624 | ||
625 | .global setoxm1 | |
626 | setoxm1: | |
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 | ||
637 | EM1CON1: | |
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 | ||
645 | EM1MAIN: | |
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 | |
742 | MLE63: | |
743 | |--Step 6.3 M <= 63 | |
744 | cmpil #-3,%d0 | |
745 | bges MGEN3 | |
746 | MLTN3: | |
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 | |
752 | MGEN3: | |
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 | ||
759 | EM1SCALE: | |
760 | |--Step 6.6 | |
761 | fmovel %d1,%FPCR | |
762 | fmulx SC(%a6),%fp0 | |
763 | ||
764 | bra t_frcinx | |
765 | ||
766 | EM1SM: | |
767 | |--Step 7 |X| < 1/4. | |
768 | cmpil #0x3FBE0000,%d0 | ...2^(-65) | |
769 | bges EM1POLY | |
770 | ||
771 | EM1TINY: | |
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 | ||
785 | EM12TINY: | |
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 | ||
798 | EM1POLY: | |
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 | ||
852 | EM1BIG: | |
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 |