| blob | 41b9954ee88981794b4f787101e3e2daf4095a8a |
1 .file "exp_m1.s"
4 // Copyright (c) 2000 - 2002, Intel Corporation
5 // All rights reserved.
6 //
7 // Contributed 2000 by the Intel Numerics Group, 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.
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://www.intel.com/software/products/opensource/libraries/num.htm.
39 //
40 // History
41 //==============================================================
42 // 02/02/00 Initial Version
43 // 04/04/00 Unwind support added
44 // 08/15/00 Bundle added after call to __libm_error_support to properly
45 // set [the previously overwritten] GR_Parameter_RESULT.
46 // 07/07/01 Improved speed of all paths
47 // 05/20/02 Cleaned up namespace and sf0 syntax
48 // 11/20/02 Improved speed, algorithm based on exp
50 // API
51 //==============================================================
52 // double expm1(double)
54 // Overview of operation
55 //==============================================================
56 // 1. Inputs of Nan, Inf, Zero, NatVal handled with special paths
57 //
58 // 2. |x| < 2^-60
59 // Result = x, computed by x + x*x to handle appropriate flags and rounding
60 //
61 // 3. 2^-60 <= |x| < 2^-2
62 // Result determined by 13th order Taylor series polynomial
63 // expm1f(x) = x + Q2*x^2 + ... + Q13*x^13
64 //
65 // 4. x < -48.0
66 // Here we know result is essentially -1 + eps, where eps only affects
67 // rounded result. Set I.
68 //
69 // 5. x >= 709.7827
70 // Result overflows. Set I, O, and call error support
71 //
72 // 6. 2^-2 <= x < 709.7827 or -48.0 <= x < -2^-2
73 // This is the main path. The algorithm is described below:
75 // Take the input x. w is "how many log2/128 in x?"
76 // w = x * 128/log2
77 // n = int(w)
78 // x = n log2/128 + r + delta
80 // n = 128M + index_1 + 2^4 index_2
81 // x = M log2 + (log2/128) index_1 + (log2/8) index_2 + r + delta
83 // exp(x) = 2^M 2^(index_1/128) 2^(index_2/8) exp(r) exp(delta)
84 // Construct 2^M
85 // Get 2^(index_1/128) from table_1;
86 // Get 2^(index_2/8) from table_2;
87 // Calculate exp(r) by series by 5th order polynomial
88 // r = x - n (log2/128)_high
89 // delta = - n (log2/128)_low
90 // Calculate exp(delta) as 1 + delta
93 // Special values
94 //==============================================================
95 // expm1(+0) = +0.0
96 // expm1(-0) = -0.0
98 // expm1(+qnan) = +qnan
99 // expm1(-qnan) = -qnan
100 // expm1(+snan) = +qnan
101 // expm1(-snan) = -qnan
103 // expm1(-inf) = -1.0
104 // expm1(+inf) = +inf
106 // Overflow and Underflow
107 //=======================
108 // expm1(x) = largest double normal when
109 // x = 709.7827 = 40862e42fefa39ef
110 //
111 // Underflow is handled as described in case 2 above.
114 // Registers used
115 //==============================================================
116 // Floating Point registers used:
117 // f8, input
118 // f9 -> f15, f32 -> f75
120 // General registers used:
121 // r14 -> r40
123 // Predicate registers used:
124 // p6 -> p15
126 // Assembly macros
127 //==============================================================
129 rRshf = r14
130 rAD_TB1 = r15
131 rAD_T1 = r15
132 rAD_TB2 = r16
133 rAD_T2 = r16
134 rAD_Ln2_lo = r17
135 rAD_P = r17
137 rN = r18
138 rIndex_1 = r19
139 rIndex_2_16 = r20
141 rM = r21
142 rBiased_M = r21
143 rIndex_1_16 = r22
144 rSignexp_x = r23
145 rExp_x = r24
146 rSig_inv_ln2 = r25
148 rAD_Q1 = r26
149 rAD_Q2 = r27
150 rTmp = r27
151 rExp_bias = r28
152 rExp_mask = r29
153 rRshf_2to56 = r30
155 rGt_ln = r31
156 rExp_2tom56 = r31
159 GR_SAVE_B0 = r33
160 GR_SAVE_PFS = r34
161 GR_SAVE_GP = r35
162 GR_SAVE_SP = r36
164 GR_Parameter_X = r37
165 GR_Parameter_Y = r38
166 GR_Parameter_RESULT = r39
167 GR_Parameter_TAG = r40
170 FR_X = f10
171 FR_Y = f1
172 FR_RESULT = f8
174 fRSHF_2TO56 = f6
175 fINV_LN2_2TO63 = f7
176 fW_2TO56_RSH = f9
177 f2TOM56 = f11
178 fP5 = f12
179 fP54 = f50
180 fP5432 = f50
181 fP4 = f13
182 fP3 = f14
183 fP32 = f14
184 fP2 = f15
186 fLn2_by_128_hi = f33
187 fLn2_by_128_lo = f34
189 fRSHF = f35
190 fNfloat = f36
191 fW = f37
192 fR = f38
193 fF = f39
195 fRsq = f40
196 fRcube = f41
198 f2M = f42
199 fS1 = f43
200 fT1 = f44
202 fMIN_DBL_OFLOW_ARG = f45
203 fMAX_DBL_MINUS_1_ARG = f46
204 fMAX_DBL_NORM_ARG = f47
205 fP_lo = f51
206 fP_hi = f52
207 fP = f53
208 fS = f54
210 fNormX = f56
212 fWre_urm_f8 = f57
214 fGt_pln = f58
215 fTmp = f58
217 fS2 = f59
218 fT2 = f60
219 fSm1 = f61
221 fXsq = f62
222 fX6 = f63
223 fX4 = f63
224 fQ7 = f64
225 fQ76 = f64
226 fQ7654 = f64
227 fQ765432 = f64
228 fQ6 = f65
229 fQ5 = f66
230 fQ54 = f66
231 fQ4 = f67
232 fQ3 = f68
233 fQ32 = f68
234 fQ2 = f69
235 fQD = f70
236 fQDC = f70
237 fQDCBA = f70
238 fQDCBA98 = f70
239 fQDCBA98765432 = f70
240 fQC = f71
241 fQB = f72
242 fQBA = f72
243 fQA = f73
244 fQ9 = f74
245 fQ98 = f74
246 fQ8 = f75
248 // Data tables
249 //==============================================================
251 RODATA
252 .align 16
254 // ************* DO NOT CHANGE ORDER OF THESE TABLES ********************
256 // double-extended 1/ln(2)
257 // 3fff b8aa 3b29 5c17 f0bb be87fed0691d3e88
258 // 3fff b8aa 3b29 5c17 f0bc
259 // For speed the significand will be loaded directly with a movl and setf.sig
260 // and the exponent will be bias+63 instead of bias+0. Thus subsequent
261 // computations need to scale appropriately.
262 // The constant 128/ln(2) is needed for the computation of w. This is also
263 // obtained by scaling the computations.
264 //
265 // Two shifting constants are loaded directly with movl and setf.d.
266 // 1. fRSHF_2TO56 = 1.1000..00 * 2^(63-7)
267 // This constant is added to x*1/ln2 to shift the integer part of
268 // x*128/ln2 into the rightmost bits of the significand.
269 // The result of this fma is fW_2TO56_RSH.
270 // 2. fRSHF = 1.1000..00 * 2^(63)
271 // This constant is subtracted from fW_2TO56_RSH * 2^(-56) to give
272 // the integer part of w, n, as a floating-point number.
273 // The result of this fms is fNfloat.
276 LOCAL_OBJECT_START(exp_Table_1)
277 data8 0x40862e42fefa39f0 // smallest dbl overflow arg
278 data8 0xc048000000000000 // approx largest arg for minus one result
279 data8 0x40862e42fefa39ef // largest dbl arg to give normal dbl result
280 data8 0x0 // pad
281 data8 0xb17217f7d1cf79ab , 0x00003ff7 // ln2/128 hi
282 data8 0xc9e3b39803f2f6af , 0x00003fb7 // ln2/128 lo
283 //
284 // Table 1 is 2^(index_1/128) where
285 // index_1 goes from 0 to 15
286 //
287 data8 0x8000000000000000 , 0x00003FFF
288 data8 0x80B1ED4FD999AB6C , 0x00003FFF
289 data8 0x8164D1F3BC030773 , 0x00003FFF
290 data8 0x8218AF4373FC25EC , 0x00003FFF
291 data8 0x82CD8698AC2BA1D7 , 0x00003FFF
292 data8 0x8383594EEFB6EE37 , 0x00003FFF
293 data8 0x843A28C3ACDE4046 , 0x00003FFF
294 data8 0x84F1F656379C1A29 , 0x00003FFF
295 data8 0x85AAC367CC487B15 , 0x00003FFF
296 data8 0x8664915B923FBA04 , 0x00003FFF
297 data8 0x871F61969E8D1010 , 0x00003FFF
298 data8 0x87DB357FF698D792 , 0x00003FFF
299 data8 0x88980E8092DA8527 , 0x00003FFF
300 data8 0x8955EE03618E5FDD , 0x00003FFF
301 data8 0x8A14D575496EFD9A , 0x00003FFF
302 data8 0x8AD4C6452C728924 , 0x00003FFF
303 LOCAL_OBJECT_END(exp_Table_1)
305 // Table 2 is 2^(index_1/8) where
306 // index_2 goes from 0 to 7
307 LOCAL_OBJECT_START(exp_Table_2)
308 data8 0x8000000000000000 , 0x00003FFF
309 data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF
310 data8 0x9837F0518DB8A96F , 0x00003FFF
311 data8 0xA5FED6A9B15138EA , 0x00003FFF
312 data8 0xB504F333F9DE6484 , 0x00003FFF
313 data8 0xC5672A115506DADD , 0x00003FFF
314 data8 0xD744FCCAD69D6AF4 , 0x00003FFF
315 data8 0xEAC0C6E7DD24392F , 0x00003FFF
316 LOCAL_OBJECT_END(exp_Table_2)
319 LOCAL_OBJECT_START(exp_p_table)
320 data8 0x3f8111116da21757 //P5
321 data8 0x3fa55555d787761c //P4
322 data8 0x3fc5555555555414 //P3
323 data8 0x3fdffffffffffd6a //P2
324 LOCAL_OBJECT_END(exp_p_table)
326 LOCAL_OBJECT_START(exp_Q1_table)
327 data8 0x3de6124613a86d09 // QD = 1/13!
328 data8 0x3e21eed8eff8d898 // QC = 1/12!
329 data8 0x3ec71de3a556c734 // Q9 = 1/9!
330 data8 0x3efa01a01a01a01a // Q8 = 1/8!
331 data8 0x8888888888888889,0x3ff8 // Q5 = 1/5!
332 data8 0xaaaaaaaaaaaaaaab,0x3ffc // Q3 = 1/3!
333 data8 0x0,0x0 // Pad to avoid bank conflicts
334 LOCAL_OBJECT_END(exp_Q1_table)
336 LOCAL_OBJECT_START(exp_Q2_table)
337 data8 0x3e5ae64567f544e4 // QB = 1/11!
338 data8 0x3e927e4fb7789f5c // QA = 1/10!
339 data8 0x3f2a01a01a01a01a // Q7 = 1/7!
340 data8 0x3f56c16c16c16c17 // Q6 = 1/6!
341 data8 0xaaaaaaaaaaaaaaab,0x3ffa // Q4 = 1/4!
342 data8 0x8000000000000000,0x3ffe // Q2 = 1/2!
343 LOCAL_OBJECT_END(exp_Q2_table)
346 .section .text
347 GLOBAL_IEEE754_ENTRY(expm1)
349 { .mlx
350 getf.exp rSignexp_x = f8 // Must recompute if x unorm
351 movl rSig_inv_ln2 = 0xb8aa3b295c17f0bc // signif of 1/ln2
352 }
353 { .mlx
354 addl rAD_TB1 = @ltoff(exp_Table_1), gp
355 movl rRshf_2to56 = 0x4768000000000000 // 1.10000 2^(63+56)
356 }
357 ;;
359 // We do this fnorm right at the beginning to normalize
360 // any input unnormals so that SWA is not taken.
361 { .mfi
362 ld8 rAD_TB1 = [rAD_TB1]
363 fclass.m p6,p0 = f8,0x0b // Test for x=unorm
364 mov rExp_mask = 0x1ffff
365 }
366 { .mfi
367 mov rExp_bias = 0xffff
368 fnorm.s1 fNormX = f8
369 mov rExp_2tom56 = 0xffff-56
370 }
371 ;;
373 // Form two constants we need
374 // 1/ln2 * 2^63 to compute w = x * 1/ln2 * 128
375 // 1.1000..000 * 2^(63+63-7) to right shift int(w) into the significand
377 { .mfi
378 setf.sig fINV_LN2_2TO63 = rSig_inv_ln2 // form 1/ln2 * 2^63
379 fclass.m p8,p0 = f8,0x07 // Test for x=0
380 nop.i 0
381 }
382 { .mlx
383 setf.d fRSHF_2TO56 = rRshf_2to56 // Form 1.100 * 2^(63+56)
384 movl rRshf = 0x43e8000000000000 // 1.10000 2^63 for rshift
385 }
386 ;;
388 { .mfi
389 setf.exp f2TOM56 = rExp_2tom56 // form 2^-56 for scaling Nfloat
390 fclass.m p9,p0 = f8,0x22 // Test for x=-inf
391 add rAD_TB2 = 0x140, rAD_TB1 // Point to Table 2
392 }
393 { .mib
394 add rAD_Q1 = 0x1e0, rAD_TB1 // Point to Q table for small path
395 add rAD_Ln2_lo = 0x30, rAD_TB1 // Point to ln2_by_128_lo
396 (p6) br.cond.spnt EXPM1_UNORM // Branch if x unorm
397 }
398 ;;
400 EXPM1_COMMON:
401 { .mfi
402 ldfpd fMIN_DBL_OFLOW_ARG, fMAX_DBL_MINUS_1_ARG = [rAD_TB1],16
403 fclass.m p10,p0 = f8,0x1e1 // Test for x=+inf, NaN, NaT
404 add rAD_Q2 = 0x50, rAD_Q1 // Point to Q table for small path
405 }
406 { .mfb
407 nop.m 0
408 nop.f 0
409 (p8) br.ret.spnt b0 // Exit for x=0, return x
410 }
411 ;;
413 { .mfi
414 ldfd fMAX_DBL_NORM_ARG = [rAD_TB1],16
415 nop.f 0
416 and rExp_x = rExp_mask, rSignexp_x // Biased exponent of x
417 }
418 { .mfb
419 setf.d fRSHF = rRshf // Form right shift const 1.100 * 2^63
420 (p9) fms.d.s0 f8 = f0,f0,f1 // quick exit for x=-inf
421 (p9) br.ret.spnt b0
422 }
423 ;;
425 { .mfi
426 ldfpd fQD, fQC = [rAD_Q1], 16 // Load coeff for small path
427 nop.f 0
428 sub rExp_x = rExp_x, rExp_bias // True exponent of x
429 }
430 { .mfb
431 ldfpd fQB, fQA = [rAD_Q2], 16 // Load coeff for small path
432 (p10) fma.d.s0 f8 = f8, f1, f0 // For x=+inf, NaN, NaT
433 (p10) br.ret.spnt b0 // Exit for x=+inf, NaN, NaT
434 }
435 ;;
437 { .mfi
438 ldfpd fQ9, fQ8 = [rAD_Q1], 16 // Load coeff for small path
439 fma.s1 fXsq = fNormX, fNormX, f0 // x*x for small path
440 cmp.gt p7, p8 = -2, rExp_x // Test |x| < 2^(-2)
441 }
442 { .mfi
443 ldfpd fQ7, fQ6 = [rAD_Q2], 16 // Load coeff for small path
444 nop.f 0
445 nop.i 0
446 }
447 ;;
449 { .mfi
450 ldfe fQ5 = [rAD_Q1], 16 // Load coeff for small path
451 nop.f 0
452 nop.i 0
453 }
454 { .mib
455 ldfe fQ4 = [rAD_Q2], 16 // Load coeff for small path
456 (p7) cmp.gt.unc p6, p7 = -60, rExp_x // Test |x| < 2^(-60)
457 (p7) br.cond.spnt EXPM1_SMALL // Branch if 2^-60 <= |x| < 2^-2
458 }
459 ;;
461 // W = X * Inv_log2_by_128
462 // By adding 1.10...0*2^63 we shift and get round_int(W) in significand.
463 // We actually add 1.10...0*2^56 to X * Inv_log2 to do the same thing.
465 { .mfi
466 ldfe fLn2_by_128_hi = [rAD_TB1],32
467 fma.s1 fW_2TO56_RSH = fNormX, fINV_LN2_2TO63, fRSHF_2TO56
468 nop.i 0
469 }
470 { .mfb
471 ldfe fLn2_by_128_lo = [rAD_Ln2_lo]
472 (p6) fma.d.s0 f8 = f8, f8, f8 // If x < 2^-60, result=x+x*x
473 (p6) br.ret.spnt b0 // Exit if x < 2^-60
474 }
475 ;;
477 // Divide arguments into the following categories:
478 // Certain minus one p11 - -inf < x <= MAX_DBL_MINUS_1_ARG
479 // Possible Overflow p14 - MAX_DBL_NORM_ARG < x < MIN_DBL_OFLOW_ARG
480 // Certain Overflow p15 - MIN_DBL_OFLOW_ARG <= x < +inf
481 //
482 // If the input is really a double arg, then there will never be "Possible
483 // Overflow" arguments.
484 //
486 // After that last load, rAD_TB1 points to the beginning of table 1
488 { .mfi
489 nop.m 0
490 fcmp.ge.s1 p15,p14 = fNormX,fMIN_DBL_OFLOW_ARG
491 nop.i 0
492 }
493 ;;
495 { .mfi
496 add rAD_P = 0x80, rAD_TB2
497 fcmp.le.s1 p11,p0 = fNormX,fMAX_DBL_MINUS_1_ARG
498 nop.i 0
499 }
500 ;;
502 { .mfb
503 ldfpd fP5, fP4 = [rAD_P] ,16
504 (p14) fcmp.gt.unc.s1 p14,p0 = fNormX,fMAX_DBL_NORM_ARG
505 (p15) br.cond.spnt EXPM1_CERTAIN_OVERFLOW
506 }
507 ;;
509 // Nfloat = round_int(W)
510 // The signficand of fW_2TO56_RSH contains the rounded integer part of W,
511 // as a twos complement number in the lower bits (that is, it may be negative).
512 // That twos complement number (called N) is put into rN.
514 // Since fW_2TO56_RSH is scaled by 2^56, it must be multiplied by 2^-56
515 // before the shift constant 1.10000 * 2^63 is subtracted to yield fNfloat.
516 // Thus, fNfloat contains the floating point version of N
518 { .mfb
519 ldfpd fP3, fP2 = [rAD_P]
520 fms.s1 fNfloat = fW_2TO56_RSH, f2TOM56, fRSHF
521 (p11) br.cond.spnt EXPM1_CERTAIN_MINUS_ONE
522 }
523 ;;
525 { .mfi
526 getf.sig rN = fW_2TO56_RSH
527 nop.f 0
528 nop.i 0
529 }
530 ;;
532 // rIndex_1 has index_1
533 // rIndex_2_16 has index_2 * 16
534 // rBiased_M has M
535 // rIndex_1_16 has index_1 * 16
537 // r = x - Nfloat * ln2_by_128_hi
538 // f = 1 - Nfloat * ln2_by_128_lo
539 { .mfi
540 and rIndex_1 = 0x0f, rN
541 fnma.s1 fR = fNfloat, fLn2_by_128_hi, fNormX
542 shr rM = rN, 0x7
543 }
544 { .mfi
545 and rIndex_2_16 = 0x70, rN
546 fnma.s1 fF = fNfloat, fLn2_by_128_lo, f1
547 nop.i 0
548 }
549 ;;
551 // rAD_T1 has address of T1
552 // rAD_T2 has address if T2
554 { .mmi
555 add rBiased_M = rExp_bias, rM
556 add rAD_T2 = rAD_TB2, rIndex_2_16
557 shladd rAD_T1 = rIndex_1, 4, rAD_TB1
558 }
559 ;;
561 // Create Scale = 2^M
562 // Load T1 and T2
563 { .mmi
564 setf.exp f2M = rBiased_M
565 ldfe fT2 = [rAD_T2]
566 nop.i 0
567 }
568 ;;
570 { .mfi
571 ldfe fT1 = [rAD_T1]
572 fmpy.s0 fTmp = fLn2_by_128_lo, fLn2_by_128_lo // Force inexact
573 nop.i 0
574 }
575 ;;
577 { .mfi
578 nop.m 0
579 fma.s1 fP54 = fR, fP5, fP4
580 nop.i 0
581 }
582 { .mfi
583 nop.m 0
584 fma.s1 fP32 = fR, fP3, fP2
585 nop.i 0
586 }
587 ;;
589 { .mfi
590 nop.m 0
591 fma.s1 fRsq = fR, fR, f0
592 nop.i 0
593 }
594 ;;
596 { .mfi
597 nop.m 0
598 fma.s1 fP5432 = fRsq, fP54, fP32
599 nop.i 0
600 }
601 ;;
603 { .mfi
604 nop.m 0
605 fma.s1 fS2 = fF,fT2,f0
606 nop.i 0
607 }
608 { .mfi
609 nop.m 0
610 fma.s1 fS1 = f2M,fT1,f0
611 nop.i 0
612 }
613 ;;
615 { .mfi
616 nop.m 0
617 fma.s1 fP = fRsq, fP5432, fR
618 nop.i 0
619 }
620 ;;
622 { .mfi
623 nop.m 0
624 fms.s1 fSm1 = fS1,fS2,f1 // S - 1.0
625 nop.i 0
626 }
627 { .mfb
628 nop.m 0
629 fma.s1 fS = fS1,fS2,f0
630 (p14) br.cond.spnt EXPM1_POSSIBLE_OVERFLOW
631 }
632 ;;
634 { .mfb
635 nop.m 0
636 fma.d.s0 f8 = fS, fP, fSm1
637 br.ret.sptk b0 // Normal path exit
638 }
639 ;;
641 // Here if 2^-60 <= |x| <2^-2
642 // Compute 13th order polynomial
643 EXPM1_SMALL:
644 { .mmf
645 ldfe fQ3 = [rAD_Q1], 16
646 ldfe fQ2 = [rAD_Q2], 16
647 fma.s1 fX4 = fXsq, fXsq, f0
648 }
649 ;;
651 { .mfi
652 nop.m 0
653 fma.s1 fQDC = fQD, fNormX, fQC
654 nop.i 0
655 }
656 { .mfi
657 nop.m 0
658 fma.s1 fQBA = fQB, fNormX, fQA
659 nop.i 0
660 }
661 ;;
663 { .mfi
664 nop.m 0
665 fma.s1 fQ98 = fQ9, fNormX, fQ8
666 nop.i 0
667 }
668 { .mfi
669 nop.m 0
670 fma.s1 fQ76= fQ7, fNormX, fQ6
671 nop.i 0
672 }
673 ;;
675 { .mfi
676 nop.m 0
677 fma.s1 fQ54 = fQ5, fNormX, fQ4
678 nop.i 0
679 }
680 ;;
682 { .mfi
683 nop.m 0
684 fma.s1 fX6 = fX4, fXsq, f0
685 nop.i 0
686 }
687 { .mfi
688 nop.m 0
689 fma.s1 fQ32= fQ3, fNormX, fQ2
690 nop.i 0
691 }
692 ;;
694 { .mfi
695 nop.m 0
696 fma.s1 fQDCBA = fQDC, fXsq, fQBA
697 nop.i 0
698 }
699 { .mfi
700 nop.m 0
701 fma.s1 fQ7654 = fQ76, fXsq, fQ54
702 nop.i 0
703 }
704 ;;
706 { .mfi
707 nop.m 0
708 fma.s1 fQDCBA98 = fQDCBA, fXsq, fQ98
709 nop.i 0
710 }
711 { .mfi
712 nop.m 0
713 fma.s1 fQ765432 = fQ7654, fXsq, fQ32
714 nop.i 0
715 }
716 ;;
718 { .mfi
719 nop.m 0
720 fma.s1 fQDCBA98765432 = fQDCBA98, fX6, fQ765432
721 nop.i 0
722 }
723 ;;
725 { .mfb
726 nop.m 0
727 fma.d.s0 f8 = fQDCBA98765432, fXsq, fNormX
728 br.ret.sptk b0 // Exit small branch
729 }
730 ;;
733 EXPM1_POSSIBLE_OVERFLOW:
735 // Here if fMAX_DBL_NORM_ARG < x < fMIN_DBL_OFLOW_ARG
736 // This cannot happen if input is a double, only if input higher precision.
737 // Overflow is a possibility, not a certainty.
739 // Recompute result using status field 2 with user's rounding mode,
740 // and wre set. If result is larger than largest double, then we have
741 // overflow
743 { .mfi
744 mov rGt_ln = 0x103ff // Exponent for largest dbl + 1 ulp
745 fsetc.s2 0x7F,0x42 // Get user's round mode, set wre
746 nop.i 0
747 }
748 ;;
750 { .mfi
751 setf.exp fGt_pln = rGt_ln // Create largest double + 1 ulp
752 fma.d.s2 fWre_urm_f8 = fS, fP, fSm1 // Result with wre set
753 nop.i 0
754 }
755 ;;
757 { .mfi
758 nop.m 0
759 fsetc.s2 0x7F,0x40 // Turn off wre in sf2
760 nop.i 0
761 }
762 ;;
764 { .mfi
765 nop.m 0
766 fcmp.ge.s1 p6, p0 = fWre_urm_f8, fGt_pln // Test for overflow
767 nop.i 0
768 }
769 ;;
771 { .mfb
772 nop.m 0
773 nop.f 0
774 (p6) br.cond.spnt EXPM1_CERTAIN_OVERFLOW // Branch if overflow
775 }
776 ;;
778 { .mfb
779 nop.m 0
780 fma.d.s0 f8 = fS, fP, fSm1
781 br.ret.sptk b0 // Exit if really no overflow
782 }
783 ;;
785 EXPM1_CERTAIN_OVERFLOW:
786 { .mmi
787 sub rTmp = rExp_mask, r0, 1
788 ;;
789 setf.exp fTmp = rTmp
790 nop.i 0
791 }
792 ;;
794 { .mfi
795 alloc r32=ar.pfs,1,4,4,0
796 fmerge.s FR_X = f8,f8
797 nop.i 0
798 }
799 { .mfb
800 mov GR_Parameter_TAG = 41
801 fma.d.s0 FR_RESULT = fTmp, fTmp, f0 // Set I,O and +INF result
802 br.cond.sptk __libm_error_region
803 }
804 ;;
806 // Here if x unorm
807 EXPM1_UNORM:
808 { .mfb
809 getf.exp rSignexp_x = fNormX // Must recompute if x unorm
810 fcmp.eq.s0 p6, p0 = f8, f0 // Set D flag
811 br.cond.sptk EXPM1_COMMON
812 }
813 ;;
815 // here if result will be -1 and inexact, x <= -48.0
816 EXPM1_CERTAIN_MINUS_ONE:
817 { .mmi
818 mov rTmp = 1
819 ;;
820 setf.exp fTmp = rTmp
821 nop.i 0
822 }
823 ;;
825 { .mfb
826 nop.m 0
827 fms.d.s0 FR_RESULT = fTmp, fTmp, f1 // Set I, rounded -1+eps result
828 br.ret.sptk b0
829 }
830 ;;
832 GLOBAL_IEEE754_END(expm1)
834 LOCAL_LIBM_ENTRY(__libm_error_region)
835 .prologue
836 { .mfi
837 add GR_Parameter_Y=-32,sp // Parameter 2 value
838 nop.f 0
839 .save ar.pfs,GR_SAVE_PFS
840 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
841 }
842 { .mfi
843 .fframe 64
844 add sp=-64,sp // Create new stack
845 nop.f 0
846 mov GR_SAVE_GP=gp // Save gp
847 };;
848 { .mmi
849 stfd [GR_Parameter_Y] = FR_Y,16 // STORE Parameter 2 on stack
850 add GR_Parameter_X = 16,sp // Parameter 1 address
851 .save b0, GR_SAVE_B0
852 mov GR_SAVE_B0=b0 // Save b0
853 };;
854 .body
855 { .mib
856 stfd [GR_Parameter_X] = FR_X // STORE Parameter 1 on stack
857 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address
858 nop.b 0
859 }
860 { .mib
861 stfd [GR_Parameter_Y] = FR_RESULT // STORE Parameter 3 on stack
862 add GR_Parameter_Y = -16,GR_Parameter_Y
863 br.call.sptk b0=__libm_error_support# // Call error handling function
864 };;
865 { .mmi
866 add GR_Parameter_RESULT = 48,sp
867 nop.m 0
868 nop.i 0
869 };;
870 { .mmi
871 ldfd f8 = [GR_Parameter_RESULT] // Get return result off stack
872 .restore sp
873 add sp = 64,sp // Restore stack pointer
874 mov b0 = GR_SAVE_B0 // Restore return address
875 };;
876 { .mib
877 mov gp = GR_SAVE_GP // Restore gp
878 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
879 br.ret.sptk b0 // Return
880 };;
882 LOCAL_LIBM_END(__libm_error_region)
883 .type __libm_error_support#,@function
884 .global __libm_error_support#