]> git.ipfire.org Git - thirdparty/glibc.git/blame - sysdeps/ia64/fpu/s_expm1.S
projects / thirdparty / glibc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Remove "Contributed by" lines
[thirdparty/glibc.git] / sysdeps / ia64 / fpu / s_expm1.S
CommitLineData
d5efd131
MF		 1.file "exp_m1.s"
2
3
4// Copyright (c) 2000 - 2005, Intel Corporation
5// All rights reserved.
6//
d5efd131
MF		 7//
8// Redistribution and use in source and binary forms, with or without
9// modification, are permitted provided that the following conditions are
10// met:
11//
12// * Redistributions of source code must retain the above copyright
13// notice, this list of conditions and the following disclaimer.
14//
15// * Redistributions in binary form must reproduce the above copyright
16// notice, this list of conditions and the following disclaimer in the
17// documentation and/or other materials provided with the distribution.
18//
19// * The name of Intel Corporation may not be used to endorse or promote
20// products derived from this software without specific prior written
21// permission.
22
23// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
26// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
27// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
28// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
29// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
30// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
31// OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
32// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34//
35// Intel Corporation is the author of this code, and requests that all
36// problem reports or change requests be submitted to it directly at
37// http://www.intel.com/software/products/opensource/libraries/num.htm.
38//
39// History
40//==============================================================
41// 02/02/00 Initial Version
42// 04/04/00 Unwind support added
43// 08/15/00 Bundle added after call to __libm_error_support to properly
44// set [the previously overwritten] GR_Parameter_RESULT.
45// 07/07/01 Improved speed of all paths
46// 05/20/02 Cleaned up namespace and sf0 syntax
47// 11/20/02 Improved speed, algorithm based on exp
48// 03/31/05 Reformatted delimiters between data tables
49
50// API
51//==============================================================
52// double expm1(double)
53
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//
0347518d 72// 6. 2^-2 <= x < 709.7827 or -48.0 <= x < -2^-2
d5efd131
MF		 73// This is the main path. The algorithm is described below:
74
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
79
80// n = 128M + index_1 + 2^4 index_2
81// x = M log2 + (log2/128) index_1 + (log2/8) index_2 + r + delta
82
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
91
92
93// Special values
94//==============================================================
95// expm1(+0) = +0.0
96// expm1(-0) = -0.0
97
98// expm1(+qnan) = +qnan
99// expm1(-qnan) = -qnan
100// expm1(+snan) = +qnan
101// expm1(-snan) = -qnan
102
103// expm1(-inf) = -1.0
104// expm1(+inf) = +inf
105
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.
112
113
114// Registers used
115//==============================================================
116// Floating Point registers used:
117// f8, input
118// f9 -> f15, f32 -> f75
119
120// General registers used:
121// r14 -> r40
122
123// Predicate registers used:
124// p6 -> p15
125
126// Assembly macros
127//==============================================================
128
129rRshf = r14
130rAD_TB1 = r15
131rAD_T1 = r15
132rAD_TB2 = r16
133rAD_T2 = r16
134rAD_Ln2_lo = r17
135rAD_P = r17
136
137rN = r18
138rIndex_1 = r19
139rIndex_2_16 = r20
140
141rM = r21
142rBiased_M = r21
143rIndex_1_16 = r22
144rSignexp_x = r23
145rExp_x = r24
146rSig_inv_ln2 = r25
147
148rAD_Q1 = r26
149rAD_Q2 = r27
150rTmp = r27
151rExp_bias = r28
152rExp_mask = r29
153rRshf_2to56 = r30
154
155rGt_ln = r31
156rExp_2tom56 = r31
157
158
159GR_SAVE_B0 = r33
160GR_SAVE_PFS = r34
161GR_SAVE_GP = r35
162GR_SAVE_SP = r36
163
164GR_Parameter_X = r37
165GR_Parameter_Y = r38
166GR_Parameter_RESULT = r39
167GR_Parameter_TAG = r40
168
169
170FR_X = f10
171FR_Y = f1
172FR_RESULT = f8
173
174fRSHF_2TO56 = f6
175fINV_LN2_2TO63 = f7
176fW_2TO56_RSH = f9
177f2TOM56 = f11
178fP5 = f12
179fP54 = f50
180fP5432 = f50
181fP4 = f13
182fP3 = f14
183fP32 = f14
184fP2 = f15
185
186fLn2_by_128_hi = f33
187fLn2_by_128_lo = f34
188
189fRSHF = f35
190fNfloat = f36
191fW = f37
192fR = f38
193fF = f39
194
195fRsq = f40
196fRcube = f41
197
198f2M = f42
199fS1 = f43
200fT1 = f44
201
202fMIN_DBL_OFLOW_ARG = f45
203fMAX_DBL_MINUS_1_ARG = f46
204fMAX_DBL_NORM_ARG = f47
205fP_lo = f51
206fP_hi = f52
207fP = f53
208fS = f54
209
210fNormX = f56
211
212fWre_urm_f8 = f57
213
214fGt_pln = f58
215fTmp = f58
216
217fS2 = f59
218fT2 = f60
219fSm1 = f61
220
221fXsq = f62
222fX6 = f63
223fX4 = f63
224fQ7 = f64
225fQ76 = f64
226fQ7654 = f64
227fQ765432 = f64
228fQ6 = f65
229fQ5 = f66
230fQ54 = f66
231fQ4 = f67
232fQ3 = f68
233fQ32 = f68
234fQ2 = f69
235fQD = f70
236fQDC = f70
237fQDCBA = f70
238fQDCBA98 = f70
239fQDCBA98765432 = f70
240fQC = f71
241fQB = f72
242fQBA = f72
243fQA = f73
244fQ9 = f74
245fQ98 = f74
246fQ8 = f75
247
248// Data tables
249//==============================================================
250
251RODATA
252.align 16
253
254// ************* DO NOT CHANGE ORDER OF THESE TABLES ********************
255
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.
274
275
276LOCAL_OBJECT_START(exp_Table_1)
277data8 0x40862e42fefa39f0 // smallest dbl overflow arg
278data8 0xc048000000000000 // approx largest arg for minus one result
279data8 0x40862e42fefa39ef // largest dbl arg to give normal dbl result
280data8 0x0 // pad
281data8 0xb17217f7d1cf79ab , 0x00003ff7 // ln2/128 hi
282data8 0xc9e3b39803f2f6af , 0x00003fb7 // ln2/128 lo
283//
284// Table 1 is 2^(index_1/128) where
285// index_1 goes from 0 to 15
286//
287data8 0x8000000000000000 , 0x00003FFF
288data8 0x80B1ED4FD999AB6C , 0x00003FFF
289data8 0x8164D1F3BC030773 , 0x00003FFF
290data8 0x8218AF4373FC25EC , 0x00003FFF
291data8 0x82CD8698AC2BA1D7 , 0x00003FFF
292data8 0x8383594EEFB6EE37 , 0x00003FFF
293data8 0x843A28C3ACDE4046 , 0x00003FFF
294data8 0x84F1F656379C1A29 , 0x00003FFF
295data8 0x85AAC367CC487B15 , 0x00003FFF
296data8 0x8664915B923FBA04 , 0x00003FFF
297data8 0x871F61969E8D1010 , 0x00003FFF
298data8 0x87DB357FF698D792 , 0x00003FFF
299data8 0x88980E8092DA8527 , 0x00003FFF
300data8 0x8955EE03618E5FDD , 0x00003FFF
301data8 0x8A14D575496EFD9A , 0x00003FFF
302data8 0x8AD4C6452C728924 , 0x00003FFF
303LOCAL_OBJECT_END(exp_Table_1)
304
305// Table 2 is 2^(index_1/8) where
306// index_2 goes from 0 to 7
307LOCAL_OBJECT_START(exp_Table_2)
308data8 0x8000000000000000 , 0x00003FFF
309data8 0x8B95C1E3EA8BD6E7 , 0x00003FFF
310data8 0x9837F0518DB8A96F , 0x00003FFF
311data8 0xA5FED6A9B15138EA , 0x00003FFF
312data8 0xB504F333F9DE6484 , 0x00003FFF
313data8 0xC5672A115506DADD , 0x00003FFF
314data8 0xD744FCCAD69D6AF4 , 0x00003FFF
315data8 0xEAC0C6E7DD24392F , 0x00003FFF
316LOCAL_OBJECT_END(exp_Table_2)
317
318
319LOCAL_OBJECT_START(exp_p_table)
320data8 0x3f8111116da21757 //P5
321data8 0x3fa55555d787761c //P4
322data8 0x3fc5555555555414 //P3
323data8 0x3fdffffffffffd6a //P2
324LOCAL_OBJECT_END(exp_p_table)
325
326LOCAL_OBJECT_START(exp_Q1_table)
327data8 0x3de6124613a86d09 // QD = 1/13!
328data8 0x3e21eed8eff8d898 // QC = 1/12!
329data8 0x3ec71de3a556c734 // Q9 = 1/9!
330data8 0x3efa01a01a01a01a // Q8 = 1/8!
331data8 0x8888888888888889,0x3ff8 // Q5 = 1/5!
332data8 0xaaaaaaaaaaaaaaab,0x3ffc // Q3 = 1/3!
333data8 0x0,0x0 // Pad to avoid bank conflicts
334LOCAL_OBJECT_END(exp_Q1_table)
335
336LOCAL_OBJECT_START(exp_Q2_table)
337data8 0x3e5ae64567f544e4 // QB = 1/11!
338data8 0x3e927e4fb7789f5c // QA = 1/10!
339data8 0x3f2a01a01a01a01a // Q7 = 1/7!
340data8 0x3f56c16c16c16c17 // Q6 = 1/6!
341data8 0xaaaaaaaaaaaaaaab,0x3ffa // Q4 = 1/4!
342data8 0x8000000000000000,0x3ffe // Q2 = 1/2!
343LOCAL_OBJECT_END(exp_Q2_table)
344
345
346.section .text
347GLOBAL_IEEE754_ENTRY(expm1)
348
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;;
358
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;;
372
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
376
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;;
387
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;;
399
400EXPM1_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;;
412
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;;
424
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;;
436
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;;
448
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;;
460
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.
464
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;;
476
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//
485
486// After that last load, rAD_TB1 points to the beginning of table 1
487
488{ .mfi
489 nop.m 0
490 fcmp.ge.s1 p15,p14 = fNormX,fMIN_DBL_OFLOW_ARG
491 nop.i 0
492}
493;;
494
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;;
501
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;;
508
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.
513
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
517
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;;
524
525{ .mfi
526 getf.sig rN = fW_2TO56_RSH
527 nop.f 0
528 nop.i 0
529}
530;;
531
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
536
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;;
550
551// rAD_T1 has address of T1
552// rAD_T2 has address if T2
553
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;;
560
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;;
569
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;;
576
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;;
588
589{ .mfi
590 nop.m 0
591 fma.s1 fRsq = fR, fR, f0
592 nop.i 0
593}
594;;
595
596{ .mfi
597 nop.m 0
598 fma.s1 fP5432 = fRsq, fP54, fP32
599 nop.i 0
600}
601;;
602
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;;
614
615{ .mfi
616 nop.m 0
617 fma.s1 fP = fRsq, fP5432, fR
618 nop.i 0
619}
620;;
621
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;;
633
634{ .mfb
635 nop.m 0
636 fma.d.s0 f8 = fS, fP, fSm1
637 br.ret.sptk b0 // Normal path exit
638}
639;;
640
641// Here if 2^-60 <= |x| <2^-2
642// Compute 13th order polynomial
643EXPM1_SMALL:
644{ .mmf
645 ldfe fQ3 = [rAD_Q1], 16
646 ldfe fQ2 = [rAD_Q2], 16
647 fma.s1 fX4 = fXsq, fXsq, f0
648}
649;;
650
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;;
662
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;;
674
675{ .mfi
676 nop.m 0
677 fma.s1 fQ54 = fQ5, fNormX, fQ4
678 nop.i 0
679}
680;;
681
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;;
693
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;;
705
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;;
717
718{ .mfi
719 nop.m 0
720 fma.s1 fQDCBA98765432 = fQDCBA98, fX6, fQ765432
721 nop.i 0
722}
723;;
724
725{ .mfb
726 nop.m 0
727 fma.d.s0 f8 = fQDCBA98765432, fXsq, fNormX
728 br.ret.sptk b0 // Exit small branch
729}
730;;
731
732
733EXPM1_POSSIBLE_OVERFLOW:
734
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.
738
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
742
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;;
749
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;;
756
757{ .mfi
758 nop.m 0
759 fsetc.s2 0x7F,0x40 // Turn off wre in sf2
760 nop.i 0
761}
762;;
763
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;;
770
771{ .mfb
772 nop.m 0
773 nop.f 0
774(p6) br.cond.spnt EXPM1_CERTAIN_OVERFLOW // Branch if overflow
775}
776;;
777
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;;
784
785EXPM1_CERTAIN_OVERFLOW:
786{ .mmi
787 sub rTmp = rExp_mask, r0, 1
788;;
789 setf.exp fTmp = rTmp
790 nop.i 0
791}
792;;
793
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;;
805
806// Here if x unorm
807EXPM1_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;;
814
815// here if result will be -1 and inexact, x <= -48.0
816EXPM1_CERTAIN_MINUS_ONE:
817{ .mmi
818 mov rTmp = 1
819;;
820 setf.exp fTmp = rTmp
821 nop.i 0
822}
823;;
824
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;;
831
832GLOBAL_IEEE754_END(expm1)
0609ec0a 833libm_alias_double_other (__expm1, expm1)
d5efd131
MF		 834
835
836LOCAL_LIBM_ENTRY(__libm_error_region)
837.prologue
838{ .mfi
839 add GR_Parameter_Y=-32,sp // Parameter 2 value
840 nop.f 0
841.save ar.pfs,GR_SAVE_PFS
842 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
843}
844{ .mfi
845.fframe 64
846 add sp=-64,sp // Create new stack
847 nop.f 0
848 mov GR_SAVE_GP=gp // Save gp
849};;
850{ .mmi
851 stfd [GR_Parameter_Y] = FR_Y,16 // STORE Parameter 2 on stack
852 add GR_Parameter_X = 16,sp // Parameter 1 address
853.save b0, GR_SAVE_B0
854 mov GR_SAVE_B0=b0 // Save b0
855};;
856.body
857{ .mib
858 stfd [GR_Parameter_X] = FR_X // STORE Parameter 1 on stack
859 add GR_Parameter_RESULT = 0,GR_Parameter_Y // Parameter 3 address
860 nop.b 0
861}
862{ .mib
863 stfd [GR_Parameter_Y] = FR_RESULT // STORE Parameter 3 on stack
864 add GR_Parameter_Y = -16,GR_Parameter_Y
865 br.call.sptk b0=__libm_error_support# // Call error handling function
866};;
867{ .mmi
868 add GR_Parameter_RESULT = 48,sp
869 nop.m 0
870 nop.i 0
871};;
872{ .mmi
873 ldfd f8 = [GR_Parameter_RESULT] // Get return result off stack
874.restore sp
875 add sp = 64,sp // Restore stack pointer
876 mov b0 = GR_SAVE_B0 // Restore return address
877};;
878{ .mib
879 mov gp = GR_SAVE_GP // Restore gp
880 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
881 br.ret.sptk b0 // Return
882};;
883
884LOCAL_LIBM_END(__libm_error_region)
885.type __libm_error_support#,@function
886.global __libm_error_support#
A mirror of the official glibc repository