5261 libm should stop using synonyms.h
5298 fabs is 0-sized, confuses dis(1) and others
Reviewed by: Josef 'Jeff' Sipek <jeffpc@josefsipek.net>
Approved by: Gordon Ross <gwr@nexenta.com>
1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21 /*
22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 .file "expm1.s"
30
31 #include "libm.h"
32 LIBM_ANSI_PRAGMA_WEAK(expm1,function)
33 #include "libm_synonyms.h"
34
35 .data
36 .align 4
37 .mhundred: .float -100.0
38
39 ENTRY(expm1)
40 movl 8(%esp),%ecx / ecx <-- hi_32(x)
41 andl $0x7fffffff,%ecx / ecx <-- hi_32(|x|)
42 cmpl $0x3fe62e42,%ecx / Is |x| < ln(2)?
43 jb .shortcut / If so, take a shortcut.
44 je .check_tail / |x| may be only slightly < ln(2)
45 cmpl $0x7ff00000,%ecx / hi_32(|x|) >= hi_32(INF)?
46 jae .not_finite / if so, x is not finite
47 .finite_non_special: / Here, ln(2) < |x| < INF
48 fldl 4(%esp) / push x
49
50 subl $8,%esp / save RP and set round-to-64-bits
51 fstcw (%esp)
52 movw (%esp),%ax
53 movw %ax,4(%esp)
54 orw $0x0300,%ax
55 movw %ax,(%esp)
56 fldcw (%esp)
57
58 fldl2e / push log2e }not for xtndd_dbl
59 fmulp %st,%st(1) / z = x*log2e }not for xtndd_dbl
60 fld %st(0) / duplicate stack top
61 frndint / [z],z
62 / [z] != 0, compute exp(x) and then subtract one to get expm1(x)
63 fxch / z,[z]
64 fsub %st(1),%st / z-[z],[z]
65 f2xm1 / 2**(z-[z])-1,[z]
66 / avoid spurious underflow when scaling to compute exp(x)
67 PIC_SETUP(1)
68 flds PIC_L(.mhundred)
69 PIC_WRAPUP
70 fucom %st(2) / if -100 !< [z], then use -100
71 fstsw %ax
72 sahf
73 jb .got_int_part
74 fxch %st(2)
75 .got_int_part:
76 fstp %st(0) / 2**(z-[z])-1,max([z],-100)
77 fld1 / 1,2**(z-[z])-1,max([z],-100)
78 faddp %st,%st(1) / 2**(z-[z]) ,max([z],-100)
79 fscale / exp(x) ,max([z],-100)
80 fld1 / 1,exp(x) ,max([z],-100)
81 fxch / exp(x),1 ,max([z],-100)
82 fsubp %st,%st(1) / exp(x)-1 ,max([z],-100)
83 fstp %st(1)
84
85 fstcw (%esp) / restore old RP
86 movw (%esp),%dx
87 andw $0xfcff,%dx
88 movw 4(%esp),%cx
89 andw $0x0300,%cx
90 orw %dx,%cx
91 movw %cx,(%esp)
92 fldcw (%esp)
93 add $8,%esp
94
95 ret
96
97 .check_tail:
98 movl 4(%esp),%edx / edx <-- lo_32(x)
99 cmpl $0xfefa39ef,%edx / Is |x| slightly < ln(2)?
100 ja .finite_non_special / branch if |x| slightly > ln(2)
101 .shortcut:
102 / Here, |x| < ln(2), so |z| = |x*log2(e)| < 1,
103 / whence z is in f2xm1's domain.
104 fldl 4(%esp) / push x
105 fldl2e / push log2e }not for xtndd_dbl
106 fmulp %st,%st(1) / z = x*log2e }not for xtndd_dbl
107 f2xm1 / 2**(x*log2(e))-1 = e**x - 1
108 ret
109
110 .not_finite:
111 / Here, flags still have settings from execution of
112 / cmpl $0x7ff00000,%ecx / hi_32(|x|) > hi_32(INF)?
113 ja .NaN_or_pinf / if not, x may be +/- INF
114 movl 4(%esp),%edx / edx <-- lo_32(x)
115 cmpl $0,%edx / lo_32(x) = 0?
116 jne .NaN_or_pinf / if not, x is NaN
117 movl 8(%esp),%eax / eax <-- hi_32(x)
118 andl $0x80000000,%eax / here, x is infinite, but +/-?
119 jz .NaN_or_pinf / branch if x = +INF
120 fld1 / Here, x = -inf, so return -1
121 fchs
122 ret
123
124 .NaN_or_pinf:
125 / Here, x = NaN or +inf, so load x and return immediately.
126 fldl 4(%esp)
127 fwait
128 ret
129 .align 4
130 SET_SIZE(expm1)
--- EOF ---