| blob | e53d3c8d7ccbf8d7164c599c7f834bb97db0e082 |
1 .file "exp_m1l.s"
3 // Copyright (C) 2000, 2001, Intel Corporation
4 // All rights reserved.
5 //
6 // Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
7 // and Ping Tak Peter Tang of the Computational Software Lab, Intel Corporation.
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are
11 // met:
12 //
13 // * Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
15 //
16 // * Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
19 //
20 // * The name of Intel Corporation may not be used to endorse or promote
21 // products derived from this software without specific prior written
22 // permission.
23 //
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35 //
36 // Intel Corporation is the author of this code, and requests that all
37 // problem reports or change requests be submitted to it directly at
38 // http://developer.intel.com/opensource.
39 //
40 // History
41 //==============================================================
42 // 4/04/00 Unwind support added
43 // 8/15/00 Bundle added after call to __libm_error_support to properly
44 // set [the previously overwritten] GR_Parameter_RESULT.
45 //
46 // *********************************************************************
47 //
48 // Function: Combined expl(x) and expm1l(x), where
49 // x
50 // expl(x) = e , for double-extended precision x values
51 // x
52 // expm1l(x) = e - 1 for double-extended precision x values
53 //
54 // *********************************************************************
55 //
56 // Resources Used:
57 //
58 // Floating-Point Registers: f8 (Input and Return Value)
59 // f9,f32-f61, f99-f102
60 //
61 // General Purpose Registers:
62 // r32-r61
63 // r62-r65 (Used to pass arguments to error handling routine)
64 //
65 // Predicate Registers: p6-p15
66 //
67 // *********************************************************************
68 //
69 // IEEE Special Conditions:
70 //
71 // Denormal fault raised on denormal inputs
72 // Overflow exceptions raised when appropriate for exp and expm1
73 // Underflow exceptions raised when appropriate for exp and expm1
74 // (Error Handling Routine called for overflow and Underflow)
75 // Inexact raised when appropriate by algorithm
76 //
77 // expl(inf) = inf
78 // expl(-inf) = +0
79 // expl(SNaN) = QNaN
80 // expl(QNaN) = QNaN
81 // expl(0) = 1
82 // expl(EM_special Values) = QNaN
83 // expl(inf) = inf
84 // expm1l(-inf) = -1
85 // expm1l(SNaN) = QNaN
86 // expm1l(QNaN) = QNaN
87 // expm1l(0) = 0
88 // expm1l(EM_special Values) = QNaN
89 //
90 // *********************************************************************
91 //
92 // Implementation and Algorithm Notes:
93 //
94 // ker_exp_64( in_FR : X,
95 // in_GR : Flag,
96 // in_GR : Expo_Range
97 // out_FR : Y_hi,
98 // out_FR : Y_lo,
99 // out_FR : scale,
100 // out_PR : Safe )
101 //
102 // On input, X is in register format and
103 // Flag = 0 for exp,
104 // Flag = 1 for expm1,
105 //
106 // On output, provided X and X_cor are real numbers, then
107 //
108 // scale*(Y_hi + Y_lo) approximates expl(X) if Flag is 0
109 // scale*(Y_hi + Y_lo) approximates expl(X)-1 if Flag is 1
110 //
111 // The accuracy is sufficient for a highly accurate 64 sig.
112 // bit implementation. Safe is set if there is no danger of
113 // overflow/underflow when the result is composed from scale,
114 // Y_hi and Y_lo. Thus, we can have a fast return if Safe is set.
115 // Otherwise, one must prepare to handle the possible exception
116 // appropriately. Note that SAFE not set (false) does not mean
117 // that overflow/underflow will occur; only the setting of SAFE
118 // guarantees the opposite.
119 //
120 // **** High Level Overview ****
121 //
122 // The method consists of three cases.
123 //
124 // If |X| < Tiny use case exp_tiny;
125 // else if |X| < 2^(-6) use case exp_small;
126 // else use case exp_regular;
127 //
128 // Case exp_tiny:
129 //
130 // 1 + X can be used to approximate expl(X) or expl(X+X_cor);
131 // X + X^2/2 can be used to approximate expl(X) - 1
132 //
133 // Case exp_small:
134 //
135 // Here, expl(X), expl(X+X_cor), and expl(X) - 1 can all be
136 // appproximated by a relatively simple polynomial.
137 //
138 // This polynomial resembles the truncated Taylor series
139 //
140 // expl(w) = 1 + w + w^2/2! + w^3/3! + ... + w^n/n!
141 //
142 // Case exp_regular:
143 //
144 // Here we use a table lookup method. The basic idea is that in
145 // order to compute expl(X), we accurately decompose X into
146 //
147 // X = N * log(2)/(2^12) + r, |r| <= log(2)/2^13.
148 //
149 // Hence
150 //
151 // expl(X) = 2^( N / 2^12 ) * expl(r).
152 //
153 // The value 2^( N / 2^12 ) is obtained by simple combinations
154 // of values calculated beforehand and stored in table; expl(r)
155 // is approximated by a short polynomial because |r| is small.
156 //
157 // We elaborate this method in 4 steps.
158 //
159 // Step 1: Reduction
160 //
161 // The value 2^12/log(2) is stored as a double-extended number
162 // L_Inv.
163 //
164 // N := round_to_nearest_integer( X * L_Inv )
165 //
166 // The value log(2)/2^12 is stored as two numbers L_hi and L_lo so
167 // that r can be computed accurately via
168 //
169 // r := (X - N*L_hi) - N*L_lo
170 //
171 // We pick L_hi such that N*L_hi is representable in 64 sig. bits
172 // and thus the FMA X - N*L_hi is error free. So r is the
173 // 1 rounding error from an exact reduction with respect to
174 //
175 // L_hi + L_lo.
176 //
177 // In particular, L_hi has 30 significant bit and can be stored
178 // as a double-precision number; L_lo has 64 significant bits and
179 // stored as a double-extended number.
180 //
181 // In the case Flag = 2, we further modify r by
182 //
183 // r := r + X_cor.
184 //
185 // Step 2: Approximation
186 //
187 // expl(r) - 1 is approximated by a short polynomial of the form
188 //
189 // r + A_1 r^2 + A_2 r^3 + A_3 r^4 .
190 //
191 // Step 3: Composition from Table Values
192 //
193 // The value 2^( N / 2^12 ) can be composed from a couple of tables
194 // of precalculated values. First, express N as three integers
195 // K, M_1, and M_2 as
196 //
197 // N = K * 2^12 + M_1 * 2^6 + M_2
198 //
199 // Where 0 <= M_1, M_2 < 2^6; and K can be positive or negative.
200 // When N is represented in 2's complement, M_2 is simply the 6
201 // lsb's, M_1 is the next 6, and K is simply N shifted right
202 // arithmetically (sign extended) by 12 bits.
203 //
204 // Now, 2^( N / 2^12 ) is simply
205 //
206 // 2^K * 2^( M_1 / 2^6 ) * 2^( M_2 / 2^12 )
207 //
208 // Clearly, 2^K needs no tabulation. The other two values are less
209 // trivial because if we store each accurately to more than working
210 // precision, than its product is too expensive to calculate. We
211 // use the following method.
212 //
213 // Define two mathematical values, delta_1 and delta_2, implicitly
214 // such that
215 //
216 // T_1 = expl( [M_1 log(2)/2^6] - delta_1 )
217 // T_2 = expl( [M_2 log(2)/2^12] - delta_2 )
218 //
219 // are representable as 24 significant bits. To illustrate the idea,
220 // we show how we define delta_1:
221 //
222 // T_1 := round_to_24_bits( expl( M_1 log(2)/2^6 ) )
223 // delta_1 = (M_1 log(2)/2^6) - log( T_1 )
224 //
225 // The last equality means mathematical equality. We then tabulate
226 //
227 // W_1 := expl(delta_1) - 1
228 // W_2 := expl(delta_2) - 1
229 //
230 // Both in double precision.
231 //
232 // From the tabulated values T_1, T_2, W_1, W_2, we compose the values
233 // T and W via
234 //
235 // T := T_1 * T_2 ...exactly
236 // W := W_1 + (1 + W_1)*W_2
237 //
238 // W approximates expl( delta ) - 1 where delta = delta_1 + delta_2.
239 // The mathematical product of T and (W+1) is an accurate representation
240 // of 2^(M_1/2^6) * 2^(M_2/2^12).
241 //
242 // Step 4. Reconstruction
243 //
244 // Finally, we can reconstruct expl(X), expl(X) - 1.
245 // Because
246 //
247 // X = K * log(2) + (M_1*log(2)/2^6 - delta_1)
248 // + (M_2*log(2)/2^12 - delta_2)
249 // + delta_1 + delta_2 + r ...accurately
250 // We have
251 //
252 // expl(X) ~=~ 2^K * ( T + T*[expl(delta_1+delta_2+r) - 1] )
253 // ~=~ 2^K * ( T + T*[expl(delta + r) - 1] )
254 // ~=~ 2^K * ( T + T*[(expl(delta)-1)
255 // + expl(delta)*(expl(r)-1)] )
256 // ~=~ 2^K * ( T + T*( W + (1+W)*poly(r) ) )
257 // ~=~ 2^K * ( Y_hi + Y_lo )
258 //
259 // where Y_hi = T and Y_lo = T*(W + (1+W)*poly(r))
260 //
261 // For expl(X)-1, we have
262 //
263 // expl(X)-1 ~=~ 2^K * ( Y_hi + Y_lo ) - 1
264 // ~=~ 2^K * ( Y_hi + Y_lo - 2^(-K) )
265 //
266 // and we combine Y_hi + Y_lo - 2^(-N) into the form of two
267 // numbers Y_hi + Y_lo carefully.
268 //
269 // **** Algorithm Details ****
270 //
271 // A careful algorithm must be used to realize the mathematical ideas
272 // accurately. We describe each of the three cases. We assume SAFE
273 // is preset to be TRUE.
274 //
275 // Case exp_tiny:
276 //
277 // The important points are to ensure an accurate result under
278 // different rounding directions and a correct setting of the SAFE
279 // flag.
280 //
281 // If Flag is 1, then
282 // SAFE := False ...possibility of underflow
283 // Scale := 1.0
284 // Y_hi := X
285 // Y_lo := 2^(-17000)
286 // Else
287 // Scale := 1.0
288 // Y_hi := 1.0
289 // Y_lo := X ...for different rounding modes
290 // Endif
291 //
292 // Case exp_small:
293 //
294 // Here we compute a simple polynomial. To exploit parallelism, we split
295 // the polynomial into several portions.
296 //
297 // Let r = X
298 //
299 // If Flag is not 1 ...i.e. expl( argument )
300 //
301 // rsq := r * r;
302 // r4 := rsq*rsq
303 // poly_lo := P_3 + r*(P_4 + r*(P_5 + r*P_6))
304 // poly_hi := r + rsq*(P_1 + r*P_2)
305 // Y_lo := poly_hi + r4 * poly_lo
306 // set lsb(Y_lo) to 1
307 // Y_hi := 1.0
308 // Scale := 1.0
309 //
310 // Else ...i.e. expl( argument ) - 1
311 //
312 // rsq := r * r
313 // r4 := rsq * rsq
314 // r6 := rsq * r4
315 // poly_lo := r6*(Q_5 + r*(Q_6 + r*Q_7))
316 // poly_hi := Q_1 + r*(Q_2 + r*(Q_3 + r*Q_4))
317 // Y_lo := rsq*poly_hi + poly_lo
318 // set lsb(Y_lo) to 1
319 // Y_hi := X
320 // Scale := 1.0
321 //
322 // Endif
323 //
324 // Case exp_regular:
325 //
326 // The previous description contain enough information except the
327 // computation of poly and the final Y_hi and Y_lo in the case for
328 // expl(X)-1.
329 //
330 // The computation of poly for Step 2:
331 //
332 // rsq := r*r
333 // poly := r + rsq*(A_1 + r*(A_2 + r*A_3))
334 //
335 // For the case expl(X) - 1, we need to incorporate 2^(-K) into
336 // Y_hi and Y_lo at the end of Step 4.
337 //
338 // If K > 10 then
339 // Y_lo := Y_lo - 2^(-K)
340 // Else
341 // If K < -10 then
342 // Y_lo := Y_hi + Y_lo
343 // Y_hi := -2^(-K)
344 // Else
345 // Y_hi := Y_hi - 2^(-K)
346 // End If
347 // End If
348 //
350 #include "libm_support.h"
352 #ifdef _LIBC
353 .rodata
354 #else
355 .data
356 #endif
358 .align 64
359 Constants_exp_64_Arg:
360 ASM_TYPE_DIRECTIVE(Constants_exp_64_Arg,@object)
361 data4 0x5C17F0BC,0xB8AA3B29,0x0000400B,0x00000000
362 data4 0x00000000,0xB17217F4,0x00003FF2,0x00000000
363 data4 0xF278ECE6,0xF473DE6A,0x00003FD4,0x00000000
364 // /* Inv_L, L_hi, L_lo */
365 ASM_SIZE_DIRECTIVE(Constants_exp_64_Arg)
367 .align 64
368 Constants_exp_64_Exponents:
369 ASM_TYPE_DIRECTIVE(Constants_exp_64_Exponents,@object)
370 data4 0x0000007E,0x00000000,0xFFFFFF83,0xFFFFFFFF
371 data4 0x000003FE,0x00000000,0xFFFFFC03,0xFFFFFFFF
372 data4 0x00003FFE,0x00000000,0xFFFFC003,0xFFFFFFFF
373 data4 0x00003FFE,0x00000000,0xFFFFC003,0xFFFFFFFF
374 data4 0xFFFFFFE2,0xFFFFFFFF,0xFFFFFFC4,0xFFFFFFFF
375 data4 0xFFFFFFBA,0xFFFFFFFF,0xFFFFFFBA,0xFFFFFFFF
376 ASM_SIZE_DIRECTIVE(Constants_exp_64_Exponents)
378 .align 64
379 Constants_exp_64_A:
380 ASM_TYPE_DIRECTIVE(Constants_exp_64_A,@object)
381 data4 0xB1B736A0,0xAAAAAAAB,0x00003FFA,0x00000000
382 data4 0x90CD6327,0xAAAAAAAB,0x00003FFC,0x00000000
383 data4 0xFFFFFFFF,0xFFFFFFFF,0x00003FFD,0x00000000
384 // /* Reversed */
385 ASM_SIZE_DIRECTIVE(Constants_exp_64_A)
387 .align 64
388 Constants_exp_64_P:
389 ASM_TYPE_DIRECTIVE(Constants_exp_64_P,@object)
390 data4 0x43914A8A,0xD00D6C81,0x00003FF2,0x00000000
391 data4 0x30304B30,0xB60BC4AC,0x00003FF5,0x00000000
392 data4 0x7474C518,0x88888888,0x00003FF8,0x00000000
393 data4 0x8DAE729D,0xAAAAAAAA,0x00003FFA,0x00000000
394 data4 0xAAAAAF61,0xAAAAAAAA,0x00003FFC,0x00000000
395 data4 0x000004C7,0x80000000,0x00003FFE,0x00000000
396 // /* Reversed */
397 ASM_SIZE_DIRECTIVE(Constants_exp_64_P)
399 .align 64
400 Constants_exp_64_Q:
401 ASM_TYPE_DIRECTIVE(Constants_exp_64_Q,@object)
402 data4 0xA49EF6CA,0xD00D56F7,0x00003FEF,0x00000000
403 data4 0x1C63493D,0xD00D59AB,0x00003FF2,0x00000000
404 data4 0xFB50CDD2,0xB60B60B5,0x00003FF5,0x00000000
405 data4 0x7BA68DC8,0x88888888,0x00003FF8,0x00000000
406 data4 0xAAAAAC8D,0xAAAAAAAA,0x00003FFA,0x00000000
407 data4 0xAAAAACCA,0xAAAAAAAA,0x00003FFC,0x00000000
408 data4 0x00000000,0x80000000,0x00003FFE,0x00000000
409 // /* Reversed */
410 ASM_SIZE_DIRECTIVE(Constants_exp_64_Q)
412 .align 64
413 Constants_exp_64_T1:
414 ASM_TYPE_DIRECTIVE(Constants_exp_64_T1,@object)
415 data4 0x3F800000,0x3F8164D2,0x3F82CD87,0x3F843A29
416 data4 0x3F85AAC3,0x3F871F62,0x3F88980F,0x3F8A14D5
417 data4 0x3F8B95C2,0x3F8D1ADF,0x3F8EA43A,0x3F9031DC
418 data4 0x3F91C3D3,0x3F935A2B,0x3F94F4F0,0x3F96942D
419 data4 0x3F9837F0,0x3F99E046,0x3F9B8D3A,0x3F9D3EDA
420 data4 0x3F9EF532,0x3FA0B051,0x3FA27043,0x3FA43516
421 data4 0x3FA5FED7,0x3FA7CD94,0x3FA9A15B,0x3FAB7A3A
422 data4 0x3FAD583F,0x3FAF3B79,0x3FB123F6,0x3FB311C4
423 data4 0x3FB504F3,0x3FB6FD92,0x3FB8FBAF,0x3FBAFF5B
424 data4 0x3FBD08A4,0x3FBF179A,0x3FC12C4D,0x3FC346CD
425 data4 0x3FC5672A,0x3FC78D75,0x3FC9B9BE,0x3FCBEC15
426 data4 0x3FCE248C,0x3FD06334,0x3FD2A81E,0x3FD4F35B
427 data4 0x3FD744FD,0x3FD99D16,0x3FDBFBB8,0x3FDE60F5
428 data4 0x3FE0CCDF,0x3FE33F89,0x3FE5B907,0x3FE8396A
429 data4 0x3FEAC0C7,0x3FED4F30,0x3FEFE4BA,0x3FF28177
430 data4 0x3FF5257D,0x3FF7D0DF,0x3FFA83B3,0x3FFD3E0C
431 ASM_SIZE_DIRECTIVE(Constants_exp_64_T1)
433 .align 64
434 Constants_exp_64_T2:
435 ASM_TYPE_DIRECTIVE(Constants_exp_64_T2,@object)
436 data4 0x3F800000,0x3F80058C,0x3F800B18,0x3F8010A4
437 data4 0x3F801630,0x3F801BBD,0x3F80214A,0x3F8026D7
438 data4 0x3F802C64,0x3F8031F2,0x3F803780,0x3F803D0E
439 data4 0x3F80429C,0x3F80482B,0x3F804DB9,0x3F805349
440 data4 0x3F8058D8,0x3F805E67,0x3F8063F7,0x3F806987
441 data4 0x3F806F17,0x3F8074A8,0x3F807A39,0x3F807FCA
442 data4 0x3F80855B,0x3F808AEC,0x3F80907E,0x3F809610
443 data4 0x3F809BA2,0x3F80A135,0x3F80A6C7,0x3F80AC5A
444 data4 0x3F80B1ED,0x3F80B781,0x3F80BD14,0x3F80C2A8
445 data4 0x3F80C83C,0x3F80CDD1,0x3F80D365,0x3F80D8FA
446 data4 0x3F80DE8F,0x3F80E425,0x3F80E9BA,0x3F80EF50
447 data4 0x3F80F4E6,0x3F80FA7C,0x3F810013,0x3F8105AA
448 data4 0x3F810B41,0x3F8110D8,0x3F81166F,0x3F811C07
449 data4 0x3F81219F,0x3F812737,0x3F812CD0,0x3F813269
450 data4 0x3F813802,0x3F813D9B,0x3F814334,0x3F8148CE
451 data4 0x3F814E68,0x3F815402,0x3F81599C,0x3F815F37
452 ASM_SIZE_DIRECTIVE(Constants_exp_64_T2)
454 .align 64
455 Constants_exp_64_W1:
456 ASM_TYPE_DIRECTIVE(Constants_exp_64_W1,@object)
457 data4 0x00000000,0x00000000,0x171EC4B4,0xBE384454
458 data4 0x4AA72766,0xBE694741,0xD42518F8,0xBE5D32B6
459 data4 0x3A319149,0x3E68D96D,0x62415F36,0xBE68F4DA
460 data4 0xC9C86A3B,0xBE6DDA2F,0xF49228FE,0x3E6B2E50
461 data4 0x1188B886,0xBE49C0C2,0x1A4C2F1F,0x3E64BFC2
462 data4 0x2CB98B54,0xBE6A2FBB,0x9A55D329,0x3E5DC5DE
463 data4 0x39A7AACE,0x3E696490,0x5C66DBA5,0x3E54728B
464 data4 0xBA1C7D7D,0xBE62B0DB,0x09F1AF5F,0x3E576E04
465 data4 0x1A0DD6A1,0x3E612500,0x795FBDEF,0xBE66A419
466 data4 0xE1BD41FC,0xBE5CDE8C,0xEA54964F,0xBE621376
467 data4 0x476E76EE,0x3E6370BE,0x3427EB92,0x3E390D1A
468 data4 0x2BF82BF8,0x3E1336DE,0xD0F7BD9E,0xBE5FF1CB
469 data4 0x0CEB09DD,0xBE60A355,0x0980F30D,0xBE5CA37E
470 data4 0x4C082D25,0xBE5C541B,0x3B467D29,0xBE5BBECA
471 data4 0xB9D946C5,0xBE400D8A,0x07ED374A,0xBE5E2A08
472 data4 0x365C8B0A,0xBE66CB28,0xD3403BCA,0x3E3AAD5B
473 data4 0xC7EA21E0,0x3E526055,0xE72880D6,0xBE442C75
474 data4 0x85222A43,0x3E58B2BB,0x522C42BF,0xBE5AAB79
475 data4 0x469DC2BC,0xBE605CB4,0xA48C40DC,0xBE589FA7
476 data4 0x1AA42614,0xBE51C214,0xC37293F4,0xBE48D087
477 data4 0xA2D673E0,0x3E367A1C,0x114F7A38,0xBE51BEBB
478 data4 0x661A4B48,0xBE6348E5,0x1D3B9962,0xBDF52643
479 data4 0x35A78A53,0x3E3A3B5E,0x1CECD788,0xBE46C46C
480 data4 0x7857D689,0xBE60B7EC,0xD14F1AD7,0xBE594D3D
481 data4 0x4C9A8F60,0xBE4F9C30,0x02DFF9D2,0xBE521873
482 data4 0x55E6D68F,0xBE5E4C88,0x667F3DC4,0xBE62140F
483 data4 0x3BF88747,0xBE36961B,0xC96EC6AA,0x3E602861
484 data4 0xD57FD718,0xBE3B5151,0xFC4A627B,0x3E561CD0
485 data4 0xCA913FEA,0xBE3A5217,0x9A5D193A,0x3E40A3CC
486 data4 0x10A9C312,0xBE5AB713,0xC5F57719,0x3E4FDADB
487 data4 0xDBDF59D5,0x3E361428,0x61B4180D,0x3E5DB5DB
488 data4 0x7408D856,0xBE42AD5F,0x31B2B707,0x3E2A3148
489 ASM_SIZE_DIRECTIVE(Constants_exp_64_W1)
491 .align 64
492 Constants_exp_64_W2:
493 ASM_TYPE_DIRECTIVE(Constants_exp_64_W2,@object)
494 data4 0x00000000,0x00000000,0x37A3D7A2,0xBE641F25
495 data4 0xAD028C40,0xBE68DD57,0xF212B1B6,0xBE5C77D8
496 data4 0x1BA5B070,0x3E57878F,0x2ECAE6FE,0xBE55A36A
497 data4 0x569DFA3B,0xBE620608,0xA6D300A3,0xBE53B50E
498 data4 0x223F8F2C,0x3E5B5EF2,0xD6DE0DF4,0xBE56A0D9
499 data4 0xEAE28F51,0xBE64EEF3,0x367EA80B,0xBE5E5AE2
500 data4 0x5FCBC02D,0x3E47CB1A,0x9BDAFEB7,0xBE656BA0
501 data4 0x805AFEE7,0x3E6E70C6,0xA3415EBA,0xBE6E0509
502 data4 0x49BFF529,0xBE56856B,0x00508651,0x3E66DD33
503 data4 0xC114BC13,0x3E51165F,0xC453290F,0x3E53333D
504 data4 0x05539FDA,0x3E6A072B,0x7C0A7696,0xBE47CD87
505 data4 0xEB05C6D9,0xBE668BF4,0x6AE86C93,0xBE67C3E3
506 data4 0xD0B3E84B,0xBE533904,0x556B53CE,0x3E63E8D9
507 data4 0x63A98DC8,0x3E212C89,0x032A7A22,0xBE33138F
508 data4 0xBC584008,0x3E530FA9,0xCCB93C97,0xBE6ADF82
509 data4 0x8370EA39,0x3E5F9113,0xFB6A05D8,0x3E5443A4
510 data4 0x181FEE7A,0x3E63DACD,0xF0F67DEC,0xBE62B29D
511 data4 0x3DDE6307,0x3E65C483,0xD40A24C1,0x3E5BF030
512 data4 0x14E437BE,0x3E658B8F,0xED98B6C7,0xBE631C29
513 data4 0x04CF7C71,0x3E6335D2,0xE954A79D,0x3E529EED
514 data4 0xF64A2FB8,0x3E5D9257,0x854ED06C,0xBE6BED1B
515 data4 0xD71405CB,0x3E5096F6,0xACB9FDF5,0xBE3D4893
516 data4 0x01B68349,0xBDFEB158,0xC6A463B9,0x3E628D35
517 data4 0xADE45917,0xBE559725,0x042FC476,0xBE68C29C
518 data4 0x01E511FA,0xBE67593B,0x398801ED,0xBE4A4313
519 data4 0xDA7C3300,0x3E699571,0x08062A9E,0x3E5349BE
520 data4 0x755BB28E,0x3E5229C4,0x77A1F80D,0x3E67E426
521 data4 0x6B69C352,0xBE52B33F,0x084DA57F,0xBE6B3550
522 data4 0xD1D09A20,0xBE6DB03F,0x2161B2C1,0xBE60CBC4
523 data4 0x78A2B771,0x3E56ED9C,0x9D0FA795,0xBE508E31
524 data4 0xFD1A54E9,0xBE59482A,0xB07FD23E,0xBE2A17CE
525 data4 0x17365712,0x3E68BF5C,0xB3785569,0x3E3956F9
526 ASM_SIZE_DIRECTIVE(Constants_exp_64_W2)
528 GR_SAVE_PFS = r59
529 GR_SAVE_B0 = r60
530 GR_SAVE_GP = r61
531 GR_Parameter_X = r62
532 GR_Parameter_Y = r63
533 GR_Parameter_RESULT = r64
534 GR_Parameter_TAG = r65
536 FR_X = f9
537 FR_Y = f9
538 FR_RESULT = f99
540 .section .text
541 .proc expm1l#
542 .global expm1l#
543 .align 64
544 expm1l:
545 #ifdef _LIBC
546 .global __expm1l#
547 __expm1l:
548 #endif
549 { .mii
550 alloc r32 = ar.pfs,0,30,4,0
551 (p0) add r33 = 1, r0
552 (p0) cmp.eq.unc p7, p0 = r0, r0
553 }
554 { .mbb
555 nop.m 999
556 (p0) br.cond.sptk exp_continue
557 nop.b 999 ;;
558 }
560 //
561 // Set p7 true for expm1
562 // Set Flag = r33 = 1 for expm1
563 //
565 .endp expm1l
566 ASM_SIZE_DIRECTIVE(expm1l)
568 #ifdef _LIBC
569 libm_hidden_def (__expm1l)
570 #endif
572 .section .text
573 .proc expl#
574 .global expl#
575 .align 64
576 expl:
577 #ifdef _LIBC
578 .global __ieee754_expl#
579 __ieee754_expl:
580 #endif
581 { .mii
582 alloc r32 = ar.pfs,0,30,4,0
583 (p0) add r33 = r0, r0
584 (p0) cmp.eq.unc p0, p7 = r0, r0 ;;
585 }
586 exp_continue:
587 { .mfi
588 (p0) add r32 = 2,r0
589 (p0) fnorm.s1 f9 = f8
590 nop.i 0
591 }
592 { .mfi
593 (p0) nop.m 0
594 //
595 // Set p7 false for exp
596 // Set Flag = r33 = 0 for exp
597 //
598 (p0) fclass.m.unc p6, p8 = f8, 0x1E7
599 nop.i 0;;
600 }
601 { .mfi
602 nop.m 999
603 (p0) fclass.nm.unc p9, p0 = f8, 0x1FF
604 nop.i 0
605 }
606 { .mfi
607 nop.m 999
608 (p0) mov f36 = f1
609 nop.i 999 ;;
610 }
611 { .mfb
612 nop.m 999
613 //
614 // Identify NatVals, NaNs, Infs, and Zeros.
615 // Identify EM unsupporteds.
616 // Save special input registers
617 (p0) mov f32 = f0
618 //
619 // Create FR_X_cor = 0.0
620 // GR_Flag = 0
621 // GR_Expo_Range = 2 (r32) for double-extended precision
622 // FR_Scale = 1.0
623 //
624 (p6) br.cond.spnt EXPL_64_SPECIAL ;;
625 }
626 { .mib
627 nop.m 999
628 nop.i 999
629 (p9) br.cond.spnt EXPL_64_UNSUPPORTED ;;
630 }
631 { .mfi
632 (p0) cmp.ne.unc p12, p13 = 0x01, r33
633 //
634 // Branch out for special input values
635 //
636 (p0) fcmp.lt.unc.s0 p9,p0 = f8, f0
637 (p0) cmp.eq.unc p15, p0 = r0, r0
638 }
639 { .mmi
640 nop.m 999
641 //
642 // Raise possible denormal operand exception
643 // Normalize x
644 //
645 // This function computes expl( x + x_cor)
646 // Input FR 1: FR_X
647 // Input FR 2: FR_X_cor
648 // Input GR 1: GR_Flag
649 // Input GR 2: GR_Expo_Range
650 // Output FR 3: FR_Y_hi
651 // Output FR 4: FR_Y_lo
652 // Output FR 5: FR_Scale
653 // Output PR 1: PR_Safe
654 (p0) addl r34 = @ltoff(Constants_exp_64_Arg#),gp
655 (p0) addl r40 = @ltoff(Constants_exp_64_W1#),gp
656 };;
657 //
658 // Prepare to load constants
659 // Set Safe = True
660 //
662 { .mmi
663 ld8 r34 = [r34]
664 ld8 r40 = [r40]
665 (p0) addl r41 = @ltoff(Constants_exp_64_W2#),gp
666 };;
668 { .mmi
669 (p0) ldfe f37 = [r34],16
670 (p0) ld8 r41 = [r41] ;;
671 }
673 //
674 // N = fcvt.fx(float_N)
675 // Set p14 if -6 > expo_X
676 //
677 //
678 // Bias = 0x0FFFF
679 // expo_X = expo_X and Mask
680 //
682 { .mmi
683 (p0) ldfe f40 = [r34],16
684 nop.m 999
685 //
686 // Load L_lo
687 // Set p10 if 14 < expo_X
688 //
689 (p0) addl r50 = @ltoff(Constants_exp_64_T1#),gp
690 }
691 { .mmi
692 nop.m 999
693 nop.m 999
694 (p0) addl r51 = @ltoff(Constants_exp_64_T2#),gp ;;
695 }
696 //
697 // Load W2_ptr
698 // Branch to SMALL is expo_X < -6
699 //
701 {.mmi
702 (p0) ld8 r50 = [r50]
703 (p0) ld8 r51 = [r51]
704 };;
706 { .mlx
707 (p0) ldfe f41 = [r34],16
708 //
709 // float_N = X * L_Inv
710 // expo_X = exponent of X
711 // Mask = 0x1FFFF
712 //
713 (p0) movl r58 = 0x0FFFF
714 }
715 { .mlx
716 nop.m 999
717 (p0) movl r39 = 0x1FFFF ;;
718 }
719 { .mmi
720 (p0) getf.exp r37 = f9
721 nop.m 999
722 (p0) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp ;;
723 }
724 { .mii
725 (p0) ld8 r34 = [r34]
726 nop.i 999
727 (p0) and r37 = r37, r39 ;;
728 }
729 { .mmi
730 (p0) sub r37 = r37, r58 ;;
731 (p0) cmp.gt.unc p14, p0 = -6, r37
732 (p0) cmp.lt.unc p10, p0 = 14, r37 ;;
733 }
734 { .mfi
735 (p0) nop.m 0
736 //
737 // Load L_inv
738 // Set p12 true for Flag = 0 (exp)
739 // Set p13 true for Flag = 1 (expm1)
740 //
741 (p0) fmpy.s1 f38 = f9, f37
742 nop.i 999 ;;
743 }
744 { .mfb
745 nop.m 999
746 //
747 // Load L_hi
748 // expo_X = expo_X - Bias
749 // get W1_ptr
750 //
751 (p0) fcvt.fx.s1 f39 = f38
752 (p14) br.cond.spnt EXPL_SMALL ;;
753 }
754 { .mib
755 nop.m 999
756 nop.i 999
757 (p10) br.cond.spnt EXPL_HUGE ;;
758 }
759 { .mmi
760 (p0) shladd r34 = r32,4,r34
761 nop.m 999
762 (p0) addl r35 = @ltoff(Constants_exp_64_A#),gp ;;
763 }
764 //
765 // Load T_1,T_2
766 //
767 { .mmi
768 nop.m 999
769 ld8 r35 =[r35]
770 nop.i 99
771 };;
772 { .mmb
773 (p0) ldfe f51 = [r35],16
774 (p0) ld8 r45 = [r34],8
775 nop.b 999 ;;
776 }
777 //
778 // Set Safe = True if k >= big_expo_neg
779 // Set Safe = False if k < big_expo_neg
780 //
781 { .mmb
782 (p0) ldfe f49 = [r35],16
783 (p0) ld8 r48 = [r34],0
784 nop.b 999 ;;
785 }
786 { .mfi
787 nop.m 999
788 //
789 // Branch to HUGE is expo_X > 14
790 //
791 (p0) fcvt.xf f38 = f39
792 nop.i 999 ;;
793 }
794 { .mfi
795 (p0) getf.sig r52 = f39
796 nop.f 999
797 nop.i 999 ;;
798 }
799 { .mii
800 nop.m 999
801 (p0) extr.u r43 = r52, 6, 6 ;;
802 //
803 // r = r - float_N * L_lo
804 // K = extr(N_fix,12,52)
805 //
806 (p0) shladd r40 = r43,3,r40 ;;
807 }
808 { .mfi
809 (p0) shladd r50 = r43,2,r50
810 (p0) fnma.s1 f42 = f40, f38, f9
811 //
812 // float_N = float(N)
813 // N_fix = signficand N
814 //
815 (p0) extr.u r42 = r52, 0, 6
816 }
817 { .mmi
818 (p0) ldfd f43 = [r40],0 ;;
819 (p0) shladd r41 = r42,3,r41
820 (p0) shladd r51 = r42,2,r51
821 }
822 //
823 // W_1_p1 = 1 + W_1
824 //
825 { .mmi
826 (p0) ldfs f44 = [r50],0 ;;
827 (p0) ldfd f45 = [r41],0
828 //
829 // M_2 = extr(N_fix,0,6)
830 // M_1 = extr(N_fix,6,6)
831 // r = X - float_N * L_hi
832 //
833 (p0) extr r44 = r52, 12, 52
834 }
835 { .mmi
836 (p0) ldfs f46 = [r51],0 ;;
837 (p0) sub r46 = r58, r44
838 (p0) cmp.gt.unc p8, p15 = r44, r45
839 }
840 //
841 // W = W_1 + W_1_p1*W_2
842 // Load A_2
843 // Bias_m_K = Bias - K
844 //
845 { .mii
846 (p0) ldfe f40 = [r35],16
847 //
848 // load A_1
849 // poly = A_2 + r*A_3
850 // rsq = r * r
851 // neg_2_mK = exponent of Bias_m_k
852 //
853 (p0) add r47 = r58, r44 ;;
854 //
855 // Set Safe = True if k <= big_expo_pos
856 // Set Safe = False if k > big_expo_pos
857 // Load A_3
858 //
859 (p15) cmp.lt p8,p15 = r44,r48 ;;
860 }
861 { .mmf
862 (p0) setf.exp f61 = r46
863 //
864 // Bias_p + K = Bias + K
865 // T = T_1 * T_2
866 //
867 (p0) setf.exp f36 = r47
868 (p0) fnma.s1 f42 = f41, f38, f42 ;;
869 }
870 { .mfi
871 nop.m 999
872 //
873 // Load W_1,W_2
874 // Load big_exp_pos, load big_exp_neg
875 //
876 (p0) fadd.s1 f47 = f43, f1
877 nop.i 999 ;;
878 }
879 { .mfi
880 nop.m 999
881 (p0) fma.s1 f52 = f42, f51, f49
882 nop.i 999
883 }
884 { .mfi
885 nop.m 999
886 (p0) fmpy.s1 f48 = f42, f42
887 nop.i 999 ;;
888 }
889 { .mfi
890 nop.m 999
891 (p0) fmpy.s1 f53 = f44, f46
892 nop.i 999 ;;
893 }
894 { .mfi
895 nop.m 999
896 (p0) fma.s1 f54 = f45, f47, f43
897 nop.i 999
898 }
899 { .mfi
900 nop.m 999
901 (p0) fneg f61 = f61
902 nop.i 999 ;;
903 }
904 { .mfi
905 nop.m 999
906 (p0) fma.s1 f52 = f42, f52, f40
907 nop.i 999 ;;
908 }
909 { .mfi
910 nop.m 999
911 (p0) fadd.s1 f55 = f54, f1
912 nop.i 999
913 }
914 { .mfi
915 nop.m 999
916 //
917 // W + Wp1 * poly
918 //
919 (p0) mov f34 = f53
920 nop.i 999 ;;
921 }
922 { .mfi
923 nop.m 999
924 //
925 // A_1 + r * poly
926 // Scale = setf_expl(Bias_p_k)
927 //
928 (p0) fma.s1 f52 = f48, f52, f42
929 nop.i 999 ;;
930 }
931 { .mfi
932 nop.m 999
933 //
934 // poly = r + rsq(A_1 + r*poly)
935 // Wp1 = 1 + W
936 // neg_2_mK = -neg_2_mK
937 //
938 (p0) fma.s1 f35 = f55, f52, f54
939 nop.i 999 ;;
940 }
941 { .mfb
942 nop.m 999
943 (p0) fmpy.s1 f35 = f35, f53
944 //
945 // Y_hi = T
946 // Y_lo = T * (W + Wp1*poly)
947 //
948 (p12) br.cond.sptk EXPL_MAIN ;;
949 }
950 //
951 // Branch if expl(x)
952 // Continue for expl(x-1)
953 //
954 { .mii
955 (p0) cmp.lt.unc p12, p13 = 10, r44
956 nop.i 999 ;;
957 //
958 // Set p12 if 10 < K, Else p13
959 //
960 (p13) cmp.gt.unc p13, p14 = -10, r44 ;;
961 }
962 //
963 // K > 10: Y_lo = Y_lo + neg_2_mK
964 // K <=10: Set p13 if -10 > K, Else set p14
965 //
966 { .mfi
967 (p13) cmp.eq p15, p0 = r0, r0
968 (p14) fadd.s1 f34 = f61, f34
969 nop.i 999 ;;
970 }
971 { .mfi
972 nop.m 999
973 (p12) fadd.s1 f35 = f35, f61
974 nop.i 999 ;;
975 }
976 { .mfi
977 nop.m 999
978 (p13) fadd.s1 f35 = f35, f34
979 nop.i 999
980 }
981 { .mfb
982 nop.m 999
983 //
984 // K <= 10 and K < -10, Set Safe = True
985 // K <= 10 and K < 10, Y_lo = Y_hi + Y_lo
986 // K <= 10 and K > =-10, Y_hi = Y_hi + neg_2_mk
987 //
988 (p13) mov f34 = f61
989 (p0) br.cond.sptk EXPL_MAIN ;;
990 }
991 EXPL_SMALL:
992 { .mmi
993 nop.m 999
994 (p0) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp
995 (p12) addl r35 = @ltoff(Constants_exp_64_P#),gp ;;
996 }
997 .pred.rel "mutex",p12,p13
998 { .mmi
999 (p12) ld8 r35=[r35]
1000 nop.m 999
1001 (p13) addl r35 = @ltoff(Constants_exp_64_Q#),gp
1002 };;
1003 { .mmi
1004 (p13) ld8 r35=[r35]
1005 (p0) ld8 r34=[r34]
1006 nop.i 999
1007 };;
1008 { .mfi
1009 (p0) add r34 = 0x48,r34
1010 //
1011 // Return
1012 // K <= 10 and K < 10, Y_hi = neg_2_mk
1013 //
1014 // /*******************************************************/
1015 // /*********** Branch EXPL_SMALL ************************/
1016 // /*******************************************************/
1017 (p0) mov f42 = f9
1018 nop.i 999 ;;
1019 }
1020 //
1021 // Flag = 0
1022 // r4 = rsq * rsq
1023 //
1024 { .mfi
1025 (p0) ld8 r49 =[r34],0
1026 nop.f 999
1027 nop.i 999 ;;
1028 }
1029 { .mii
1030 nop.m 999
1031 nop.i 999 ;;
1032 //
1033 // Flag = 1
1034 //
1035 (p0) cmp.lt.unc p14, p0 = r37, r49 ;;
1036 }
1037 { .mfi
1038 nop.m 999
1039 //
1040 // r = X
1041 //
1042 (p0) fmpy.s1 f48 = f42, f42
1043 nop.i 999 ;;
1044 }
1045 { .mfb
1046 nop.m 999
1047 //
1048 // rsq = r * r
1049 //
1050 (p0) fmpy.s1 f50 = f48, f48
1051 //
1052 // Is input very small?
1053 //
1054 (p14) br.cond.spnt EXPL_VERY_SMALL ;;
1055 }
1056 //
1057 // Flag_not1: Y_hi = 1.0
1058 // Flag is 1: r6 = rsq * r4
1059 //
1060 { .mfi
1061 (p12) ldfe f52 = [r35],16
1062 (p12) mov f34 = f1
1063 (p0) add r53 = 0x1,r0 ;;
1064 }
1065 { .mfi
1066 (p13) ldfe f51 = [r35],16
1067 //
1068 // Flag_not_1: Y_lo = poly_hi + r4 * poly_lo
1069 //
1070 (p13) mov f34 = f9
1071 nop.i 999 ;;
1072 }
1073 { .mmf
1074 (p12) ldfe f53 = [r35],16
1075 //
1076 // For Flag_not_1, Y_hi = X
1077 // Scale = 1
1078 // Create 0x000...01
1079 //
1080 (p0) setf.sig f37 = r53
1081 (p0) mov f36 = f1 ;;
1082 }
1083 { .mmi
1084 (p13) ldfe f52 = [r35],16 ;;
1085 (p12) ldfe f54 = [r35],16
1086 nop.i 999 ;;
1087 }
1088 { .mfi
1089 (p13) ldfe f53 = [r35],16
1090 (p13) fmpy.s1 f58 = f48, f50
1091 nop.i 999 ;;
1092 }
1093 //
1094 // Flag_not1: poly_lo = P_5 + r*P_6
1095 // Flag_1: poly_lo = Q_6 + r*Q_7
1096 //
1097 { .mmi
1098 (p13) ldfe f54 = [r35],16 ;;
1099 (p12) ldfe f55 = [r35],16
1100 nop.i 999 ;;
1101 }
1102 { .mmi
1103 (p12) ldfe f56 = [r35],16 ;;
1104 (p13) ldfe f55 = [r35],16
1105 nop.i 999 ;;
1106 }
1107 { .mmi
1108 (p12) ldfe f57 = [r35],0 ;;
1109 (p13) ldfe f56 = [r35],16
1110 nop.i 999 ;;
1111 }
1112 { .mfi
1113 (p13) ldfe f57 = [r35],0
1114 nop.f 999
1115 nop.i 999 ;;
1116 }
1117 { .mfi
1118 nop.m 999
1119 //
1120 // For Flag_not_1, load p5,p6,p1,p2
1121 // Else load p5,p6,p1,p2
1122 //
1123 (p12) fma.s1 f60 = f52, f42, f53
1124 nop.i 999 ;;
1125 }
1126 { .mfi
1127 nop.m 999
1128 (p13) fma.s1 f60 = f51, f42, f52
1129 nop.i 999 ;;
1130 }
1131 { .mfi
1132 nop.m 999
1133 (p12) fma.s1 f60 = f60, f42, f54
1134 nop.i 999 ;;
1135 }
1136 { .mfi
1137 nop.m 999
1138 (p12) fma.s1 f59 = f56, f42, f57
1139 nop.i 999 ;;
1140 }
1141 { .mfi
1142 nop.m 999
1143 (p13) fma.s1 f60 = f42, f60, f53
1144 nop.i 999 ;;
1145 }
1146 { .mfi
1147 nop.m 999
1148 (p12) fma.s1 f59 = f59, f48, f42
1149 nop.i 999 ;;
1150 }
1151 { .mfi
1152 nop.m 999
1153 //
1154 // Flag_1: poly_lo = Q_5 + r*(Q_6 + r*Q_7)
1155 // Flag_not1: poly_lo = P_4 + r*(P_5 + r*P_6)
1156 // Flag_not1: poly_hi = (P_1 + r*P_2)
1157 //
1158 (p13) fmpy.s1 f60 = f60, f58
1159 nop.i 999 ;;
1160 }
1161 { .mfi
1162 nop.m 999
1163 (p12) fma.s1 f60 = f60, f42, f55
1164 nop.i 999 ;;
1165 }
1166 { .mfi
1167 nop.m 999
1168 //
1169 // Flag_1: poly_lo = r6 *(Q_5 + ....)
1170 // Flag_not1: poly_hi = r + rsq *(P_1 + r*P_2)
1171 //
1172 (p12) fma.s1 f35 = f60, f50, f59
1173 nop.i 999
1174 }
1175 { .mfi
1176 nop.m 999
1177 (p13) fma.s1 f59 = f54, f42, f55
1178 nop.i 999 ;;
1179 }
1180 { .mfi
1181 nop.m 999
1182 //
1183 // Flag_not1: Y_lo = rsq* poly_hi + poly_lo
1184 // Flag_1: poly_lo = rsq* poly_hi + poly_lo
1185 //
1186 (p13) fma.s1 f59 = f59, f42, f56
1187 nop.i 999 ;;
1188 }
1189 { .mfi
1190 nop.m 999
1191 //
1192 // Flag_not_1: (P_1 + r*P_2)
1193 //
1194 (p13) fma.s1 f59 = f59, f42, f57
1195 nop.i 999 ;;
1196 }
1197 { .mfi
1198 nop.m 999
1199 //
1200 // Flag_not_1: poly_hi = r + rsq * (P_1 + r*P_2)
1201 //
1202 (p13) fma.s1 f35 = f59, f48, f60
1203 nop.i 999 ;;
1204 }
1205 { .mfi
1206 nop.m 999
1207 //
1208 // Create 0.000...01
1209 //
1210 (p0) for f37 = f35, f37
1211 nop.i 999 ;;
1212 }
1213 { .mfb
1214 nop.m 999
1215 //
1216 // Set lsb of Y_lo to 1
1217 //
1218 (p0) fmerge.se f35 = f35,f37
1219 (p0) br.cond.sptk EXPL_MAIN ;;
1220 }
1221 EXPL_VERY_SMALL:
1222 { .mmi
1223 nop.m 999
1224 nop.m 999
1225 (p13) addl r34 = @ltoff(Constants_exp_64_Exponents#),gp
1226 }
1227 { .mfi
1228 nop.m 999
1229 (p12) mov f35 = f9
1230 nop.i 999 ;;
1231 }
1232 { .mfb
1233 (p13) ld8 r34 = [r34]
1234 (p12) mov f34 = f1
1235 (p12) br.cond.sptk EXPL_MAIN ;;
1236 }
1237 { .mlx
1238 (p13) add r34 = 8,r34
1239 (p13) movl r39 = 0x0FFFE ;;
1240 }
1241 //
1242 // Load big_exp_neg
1243 // Create 1/2's exponent
1244 //
1245 { .mii
1246 (p13) setf.exp f56 = r39
1247 (p13) shladd r34 = r32,4,r34 ;;
1248 nop.i 999
1249 }
1250 //
1251 // Negative exponents are stored after positive
1252 //
1253 { .mfi
1254 (p13) ld8 r45 = [r34],0
1255 //
1256 // Y_hi = x
1257 // Scale = 1
1258 //
1259 (p13) fmpy.s1 f35 = f9, f9
1260 nop.i 999 ;;
1261 }
1262 { .mfi
1263 nop.m 999
1264 //
1265 // Reset Safe if necessary
1266 // Create 1/2
1267 //
1268 (p13) mov f34 = f9
1269 nop.i 999 ;;
1270 }
1271 { .mfi
1272 (p13) cmp.lt.unc p0, p15 = r37, r45
1273 (p13) mov f36 = f1
1274 nop.i 999 ;;
1275 }
1276 { .mfb
1277 nop.m 999
1278 //
1279 // Y_lo = x * x
1280 //
1281 (p13) fmpy.s1 f35 = f35, f56
1282 //
1283 // Y_lo = x*x/2
1284 //
1285 (p13) br.cond.sptk EXPL_MAIN ;;
1286 }
1287 EXPL_HUGE:
1288 { .mfi
1289 nop.m 999
1290 (p0) fcmp.gt.unc.s1 p14, p0 = f9, f0
1291 nop.i 999
1292 }
1293 { .mlx
1294 nop.m 999
1295 (p0) movl r39 = 0x15DC0 ;;
1296 }
1297 { .mfi
1298 (p14) setf.exp f34 = r39
1299 (p14) mov f35 = f1
1300 (p14) cmp.eq p0, p15 = r0, r0 ;;
1301 }
1302 { .mfb
1303 nop.m 999
1304 (p14) mov f36 = f34
1305 //
1306 // If x > 0, Set Safe = False
1307 // If x > 0, Y_hi = 2**(24,000)
1308 // If x > 0, Y_lo = 1.0
1309 // If x > 0, Scale = 2**(24,000)
1310 //
1311 (p14) br.cond.sptk EXPL_MAIN ;;
1312 }
1313 { .mlx
1314 nop.m 999
1315 (p12) movl r39 = 0xA240
1316 }
1317 { .mlx
1318 nop.m 999
1319 (p12) movl r38 = 0xA1DC ;;
1320 }
1321 { .mmb
1322 (p13) cmp.eq p15, p14 = r0, r0
1323 (p12) setf.exp f34 = r39
1324 nop.b 999 ;;
1325 }
1326 { .mlx
1327 (p12) setf.exp f35 = r38
1328 (p13) movl r39 = 0xFF9C
1329 }
1330 { .mfi
1331 nop.m 999
1332 (p13) fsub.s1 f34 = f0, f1
1333 nop.i 999 ;;
1334 }
1335 { .mfi
1336 nop.m 999
1337 (p12) mov f36 = f34
1338 (p12) cmp.eq p0, p15 = r0, r0 ;;
1339 }
1340 { .mfi
1341 (p13) setf.exp f35 = r39
1342 (p13) mov f36 = f1
1343 nop.i 999 ;;
1344 }
1345 EXPL_MAIN:
1346 { .mfi
1347 (p0) cmp.ne.unc p12, p0 = 0x01, r33
1348 (p0) fmpy.s1 f101 = f36, f35
1349 nop.i 999 ;;
1350 }
1351 { .mfb
1352 nop.m 999
1353 (p0) fma.s0 f99 = f34, f36, f101
1354 (p15) br.cond.sptk EXPL_64_RETURN ;;
1355 }
1356 { .mfi
1357 nop.m 999
1358 (p0) fsetc.s3 0x7F,0x01
1359 nop.i 999
1360 }
1361 { .mlx
1362 nop.m 999
1363 (p0) movl r50 = 0x00000000013FFF ;;
1364 }
1365 //
1366 // S0 user supplied status
1367 // S2 user supplied status + WRE + TD (Overflows)
1368 // S3 user supplied status + RZ + TD (Underflows)
1369 //
1370 //
1371 // If (Safe) is true, then
1372 // Compute result using user supplied status field.
1373 // No overflow or underflow here, but perhaps inexact.
1374 // Return
1375 // Else
1376 // Determine if overflow or underflow was raised.
1377 // Fetch +/- overflow threshold for IEEE single, double,
1378 // double extended
1379 //
1380 { .mfi
1381 (p0) setf.exp f60 = r50
1382 (p0) fma.s3 f102 = f34, f36, f101
1383 nop.i 999
1384 }
1385 { .mfi
1386 nop.m 999
1387 (p0) fsetc.s3 0x7F,0x40
1388 nop.i 999 ;;
1389 }
1390 { .mfi
1391 nop.m 999
1392 //
1393 // For Safe, no need to check for over/under.
1394 // For expm1, handle errors like exp.
1395 //
1396 (p0) fsetc.s2 0x7F,0x42
1397 nop.i 999;;
1398 }
1399 { .mfi
1400 nop.m 999
1401 (p0) fma.s2 f100 = f34, f36, f101
1402 nop.i 999 ;;
1403 }
1404 { .mfi
1405 nop.m 999
1406 (p0) fsetc.s2 0x7F,0x40
1407 nop.i 999 ;;
1408 }
1409 { .mfi
1410 nop.m 999
1411 (p7) fclass.m.unc p12, p0 = f102, 0x00F
1412 nop.i 999
1413 }
1414 { .mfi
1415 nop.m 999
1416 (p0) fclass.m.unc p11, p0 = f102, 0x00F
1417 nop.i 999 ;;
1418 }
1419 { .mfi
1420 nop.m 999
1421 (p7) fcmp.ge.unc.s1 p10, p0 = f100, f60
1422 nop.i 999
1423 }
1424 { .mfi
1425 nop.m 999
1426 //
1427 // Create largest double exponent + 1.
1428 // Create smallest double exponent - 1.
1429 //
1430 (p0) fcmp.ge.unc.s1 p8, p0 = f100, f60
1431 nop.i 999 ;;
1432 }
1433 //
1434 // fcmp: resultS2 >= + overflow threshold -> set (a) if true
1435 // fcmp: resultS2 <= - overflow threshold -> set (b) if true
1436 // fclass: resultS3 is denorm/unorm/0 -> set (d) if true
1437 //
1438 { .mib
1439 (p10) mov GR_Parameter_TAG = 39
1440 nop.i 999
1441 (p10) br.cond.sptk __libm_error_region ;;
1442 }
1443 { .mib
1444 (p8) mov GR_Parameter_TAG = 12
1445 nop.i 999
1446 (p8) br.cond.sptk __libm_error_region ;;
1447 }
1448 //
1449 // Report that exp overflowed
1450 //
1451 { .mib
1452 (p12) mov GR_Parameter_TAG = 40
1453 nop.i 999
1454 (p12) br.cond.sptk __libm_error_region ;;
1455 }
1456 { .mib
1457 (p11) mov GR_Parameter_TAG = 13
1458 nop.i 999
1459 (p11) br.cond.sptk __libm_error_region ;;
1460 }
1461 { .mib
1462 nop.m 999
1463 nop.i 999
1464 //
1465 // Report that exp underflowed
1466 //
1467 (p0) br.cond.sptk EXPL_64_RETURN ;;
1468 }
1469 EXPL_64_SPECIAL:
1470 { .mfi
1471 nop.m 999
1472 (p0) fclass.m.unc p6, p0 = f8, 0x0c3
1473 nop.i 999
1474 }
1475 { .mfi
1476 nop.m 999
1477 (p0) fclass.m.unc p13, p8 = f8, 0x007
1478 nop.i 999 ;;
1479 }
1480 { .mfi
1481 nop.m 999
1482 (p7) fclass.m.unc p14, p0 = f8, 0x007
1483 nop.i 999
1484 }
1485 { .mfi
1486 nop.m 999
1487 (p0) fclass.m.unc p12, p9 = f8, 0x021
1488 nop.i 999 ;;
1489 }
1490 { .mfi
1491 nop.m 999
1492 (p0) fclass.m.unc p11, p0 = f8, 0x022
1493 nop.i 999
1494 }
1495 { .mfi
1496 nop.m 999
1497 (p7) fclass.m.unc p10, p0 = f8, 0x022
1498 nop.i 999 ;;
1499 }
1500 { .mfi
1501 nop.m 999
1502 //
1503 // Identify +/- 0, Inf, or -Inf
1504 // Generate the right kind of NaN.
1505 //
1506 (p13) fadd.s0 f99 = f0, f1
1507 nop.i 999 ;;
1508 }
1509 { .mfi
1510 nop.m 999
1511 (p14) mov f99 = f8
1512 nop.i 999 ;;
1513 }
1514 { .mfb
1515 nop.m 999
1516 (p6) fadd.s0 f99 = f8, f1
1517 //
1518 // expl(+/-0) = 1
1519 // expm1l(+/-0) = +/-0
1520 // No exceptions raised
1521 //
1522 (p6) br.cond.sptk EXPL_64_RETURN ;;
1523 }
1524 { .mib
1525 nop.m 999
1526 nop.i 999
1527 (p14) br.cond.sptk EXPL_64_RETURN ;;
1528 }
1529 { .mfi
1530 nop.m 999
1531 (p11) mov f99 = f0
1532 nop.i 999 ;;
1533 }
1534 { .mfb
1535 nop.m 999
1536 (p10) fsub.s1 f99 = f0, f1
1537 //
1538 // expl(-Inf) = 0
1539 // expm1l(-Inf) = -1
1540 // No exceptions raised.
1541 //
1542 (p10) br.cond.sptk EXPL_64_RETURN ;;
1543 }
1544 { .mfb
1545 nop.m 999
1546 (p12) fmpy.s1 f99 = f8, f1
1547 //
1548 // expl(+Inf) = Inf
1549 // No exceptions raised.
1550 //
1551 (p0) br.cond.sptk EXPL_64_RETURN ;;
1552 }
1553 EXPL_64_UNSUPPORTED:
1554 { .mfb
1555 nop.m 999
1556 (p0) fmpy.s0 f99 = f8, f0
1557 (p0) br.cond.sptk EXPL_64_RETURN ;;
1558 }
1559 EXPL_64_RETURN:
1560 { .mfb
1561 nop.m 999
1562 (p0) mov f8 = f99
1563 (p0) br.ret.sptk b0
1564 }
1565 .endp
1566 ASM_SIZE_DIRECTIVE(expl)
1568 .proc __libm_error_region
1569 __libm_error_region:
1570 .prologue
1571 { .mfi
1572 add GR_Parameter_Y=-32,sp // Parameter 2 value
1573 nop.f 0
1574 .save ar.pfs,GR_SAVE_PFS
1575 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
1576 }
1577 { .mfi
1578 .fframe 64
1579 add sp=-64,sp // Create new stack
1580 nop.f 0
1581 mov GR_SAVE_GP=gp // Save gp
1582 };;
1583 { .mmi
1584 stfe [GR_Parameter_Y] = FR_Y,16 // Save Parameter 2 on stack
1585 add GR_Parameter_X = 16,sp // Parameter 1 address
1586 .save b0, GR_SAVE_B0
1587 mov GR_SAVE_B0=b0 // Save b0
1588 };;
1589 .body
1590 { .mib
1591 stfe [GR_Parameter_X] = FR_X // Store Parameter 1 on stack
1592 add GR_Parameter_RESULT = 0,GR_Parameter_Y
1593 nop.b 0 // Parameter 3 address
1594 }
1595 { .mib
1596 stfe [GR_Parameter_Y] = FR_RESULT // Store Parameter 3 on stack
1597 add GR_Parameter_Y = -16,GR_Parameter_Y
1598 br.call.sptk b0=__libm_error_support# // Call error handling function
1599 };;
1600 { .mmi
1601 nop.m 0
1602 nop.m 0
1603 add GR_Parameter_RESULT = 48,sp
1604 };;
1605 { .mmi
1606 ldfe f8 = [GR_Parameter_RESULT] // Get return result off stack
1607 .restore sp
1608 add sp = 64,sp // Restore stack pointer
1609 mov b0 = GR_SAVE_B0 // Restore return address
1610 };;
1611 { .mib
1612 mov gp = GR_SAVE_GP // Restore gp
1613 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
1614 br.ret.sptk b0 // Return
1615 };;
1616 .endp __libm_error_region
1617 ASM_SIZE_DIRECTIVE(__libm_error_region)
1619 .type __libm_error_support#,@function
1620 .global __libm_error_support#