1 /* SSE2 version of __ieee754_expf and __expf_finite
2 Copyright (C) 2012-2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <http://www.gnu.org/licenses/>. */
22 /* Short algorithm description:
24 * Let K = 64 (table size).
25 * e^x = 2^(x/log(2)) = 2^n * T[j] * (1 + P(y))
27 * x = m*log(2)/K + y, y in [0.0..log(2)/K]
28 * m = n*K + j, m,n,j - signed integer, j in [0..K-1]
29 * values of 2^(j/K) are tabulated as T[j].
31 * P(y) is a minimax polynomial approximation of expf(x)-1
32 * on small interval [0.0..log(2)/K].
34 * P(y) = P3*y*y*y*y + P2*y*y*y + P1*y*y + P0*y, calculated as
35 * z = y*y; P(y) = (P3*z + P1)*z + (P2*z + P0)*y
38 * __ieee754_expf_sse2(NaN) = NaN
39 * __ieee754_expf_sse2(+INF) = +INF
40 * __ieee754_expf_sse2(-INF) = 0
41 * __ieee754_expf_sse2(x) = 1 for subnormals
42 * for finite argument, only __ieee754_expf_sse2(0)=1 is exact
43 * __ieee754_expf_sse2(x) overflows if x>700
44 * __ieee754_expf_sse2(x) underflows if x<-700
47 * For |x|<700, __ieee754_expf_sse2 computes result in double precision,
48 * with accuracy a bit more than needed for expf, and does not round it
49 * to single precision.
54 # define MO1(symbol) L(symbol)##@GOTOFF(%edx)
55 # define MO2(symbol,reg2,_scale) L(symbol)##@GOTOFF(%edx,reg2,_scale)
57 # define MO1(symbol) L(symbol)
58 # define MO2(symbol,reg2,_scale) L(symbol)(,reg2,_scale)
62 ENTRY(__ieee754_expf_sse2)
63 /* Input: single precision x on stack at address 4(%esp) */
69 cvtss2sd 4(%esp), %xmm1 /* Convert x to double precision */
70 mov 4(%esp), %ecx /* Copy x */
71 movsd MO1(DP_KLN2), %xmm2 /* DP K/log(2) */
72 movsd MO1(DP_P2), %xmm3 /* DP P2 */
73 movl %ecx, %eax /* x */
74 mulsd %xmm1, %xmm2 /* DP x*K/log(2) */
75 andl $0x7fffffff, %ecx /* |x| */
76 cmpl $0x442f0000, %ecx /* |x|<700 ? */
77 movsd MO1(DP_P3), %xmm4 /* DP P3 */
78 addsd MO1(DP_RS), %xmm2 /* DP x*K/log(2)+RS */
82 cmpl $0x31800000, %ecx /* |x|<2^(-28) ? */
85 /* Main path: here if 2^(-28)<=|x|<700 */
86 cvtsd2ss %xmm2, %xmm2 /* SP x*K/log(2)+RS */
87 movd %xmm2, %eax /* bits of n*K+j with trash */
88 subss MO1(SP_RS), %xmm2 /* SP t=round(x*K/log(2)) */
89 movl %eax, %ecx /* n*K+j with trash */
90 cvtss2sd %xmm2, %xmm2 /* DP t */
91 andl $0x3f, %eax /* bits of j */
92 mulsd MO1(DP_NLN2K), %xmm2 /* DP -t*log(2)/K */
93 andl $0xffffffc0, %ecx /* bits of n */
95 vaddsd %xmm1, %xmm2, %xmm0 /* DP y=x-t*log(2)/K */
96 vmulsd %xmm0, %xmm0, %xmm2 /* DP z=y*y */
98 addsd %xmm1, %xmm2 /* DP y=x-t*log(2)/K */
99 movaps %xmm2, %xmm0 /* DP y */
100 mulsd %xmm2, %xmm2 /* DP z=y*y */
102 mulsd %xmm2, %xmm4 /* DP P3*z */
103 addl $0xffc0, %ecx /* bits of n + DP exponent bias */
104 mulsd %xmm2, %xmm3 /* DP P2*z */
105 shrl $2, %ecx /* High 2 bytes of DP 2^n */
106 pxor %xmm1, %xmm1 /* clear %xmm1 */
107 addsd MO1(DP_P1), %xmm4 /* DP P3*z+P1 */
108 addsd MO1(DP_P0), %xmm3 /* DP P2*z+P0 */
109 pinsrw $3, %ecx, %xmm1 /* DP 2^n */
110 mulsd %xmm2, %xmm4 /* DP (P3*z+P1)*z */
111 mulsd %xmm3, %xmm0 /* DP (P2*z+P0)*y */
112 addsd %xmm4, %xmm0 /* DP P(y) */
113 mulsd MO2(DP_T,%eax,8), %xmm0 /* DP P(y)*T[j] */
114 addsd MO2(DP_T,%eax,8), %xmm0 /* DP T[j]*(P(y)+1) */
115 mulsd %xmm1, %xmm0 /* DP result=2^n*(T[j]*(P(y)+1)) */
116 cvtsd2ss %xmm0, %xmm1
118 lea -4(%esp), %esp /* Borrow 4 bytes of stack frame */
119 movss %xmm1, 0(%esp) /* Move result from sse... */
120 flds 0(%esp) /* ...to FPU. */
121 lea 4(%esp), %esp /* Return back 4 bytes of stack frame */
126 /* Here if 0<=|x|<2^(-28) */
127 movss 4(%esp), %xmm0 /* load x */
128 addss MO1(SP_ONE), %xmm0 /* 1.0 + x */
129 /* Return 1.0 with inexact raised, except for x==0 */
134 /* Here if x is NaN, or Inf, or finite |x|>=700 */
135 movss 4(%esp), %xmm0 /* load x */
137 cmpl $0x7f800000, %ecx /* |x| is finite ? */
138 jae L(arg_inf_or_nan)
140 /* Here if finite |x|>=700 */
141 testl $0x80000000, %eax /* sign of x nonzero ? */
144 /* Here if finite x<=-700 */
145 movss MO1(SP_SMALL), %xmm0 /* load small value 2^(-100) */
146 mulss %xmm0, %xmm0 /* Return underflowed result (zero or subnormal) */
151 /* Here if finite x>=700 */
152 movss MO1(SP_LARGE), %xmm0 /* load large value 2^100 */
153 mulss %xmm0, %xmm0 /* Return overflowed result (Inf or max normal) */
158 /* Here if |x| is Inf or NAN */
159 jne L(arg_nan) /* |x| is Inf ? */
161 /* Here if |x| is Inf */
162 shrl $31, %eax /* Get sign bit of x */
163 movss MO2(SP_INF_0,%eax,4), %xmm0/* return zero or Inf, depending on sign of x */
168 /* Here if |x| is NaN */
169 addss %xmm0, %xmm0 /* Return x+x (raise invalid) */
173 lea -4(%esp), %esp /* Borrow 4 bytes of stack frame */
174 movss %xmm0, 0(%esp) /* Move result from sse... */
175 flds 0(%esp) /* ...to FPU. */
176 lea 4(%esp), %esp /* Return back 4 bytes of stack frame */
178 END(__ieee754_expf_sse2)
180 .section .rodata, "a"
182 L(DP_T): /* table of double precision values 2^(j/K) for j=[0..K-1] */
183 .long 0x00000000, 0x3ff00000
184 .long 0x3e778061, 0x3ff02c9a
185 .long 0xd3158574, 0x3ff059b0
186 .long 0x18759bc8, 0x3ff08745
187 .long 0x6cf9890f, 0x3ff0b558
188 .long 0x32d3d1a2, 0x3ff0e3ec
189 .long 0xd0125b51, 0x3ff11301
190 .long 0xaea92de0, 0x3ff1429a
191 .long 0x3c7d517b, 0x3ff172b8
192 .long 0xeb6fcb75, 0x3ff1a35b
193 .long 0x3168b9aa, 0x3ff1d487
194 .long 0x88628cd6, 0x3ff2063b
195 .long 0x6e756238, 0x3ff2387a
196 .long 0x65e27cdd, 0x3ff26b45
197 .long 0xf51fdee1, 0x3ff29e9d
198 .long 0xa6e4030b, 0x3ff2d285
199 .long 0x0a31b715, 0x3ff306fe
200 .long 0xb26416ff, 0x3ff33c08
201 .long 0x373aa9cb, 0x3ff371a7
202 .long 0x34e59ff7, 0x3ff3a7db
203 .long 0x4c123422, 0x3ff3dea6
204 .long 0x21f72e2a, 0x3ff4160a
205 .long 0x6061892d, 0x3ff44e08
206 .long 0xb5c13cd0, 0x3ff486a2
207 .long 0xd5362a27, 0x3ff4bfda
208 .long 0x769d2ca7, 0x3ff4f9b2
209 .long 0x569d4f82, 0x3ff5342b
210 .long 0x36b527da, 0x3ff56f47
211 .long 0xdd485429, 0x3ff5ab07
212 .long 0x15ad2148, 0x3ff5e76f
213 .long 0xb03a5585, 0x3ff6247e
214 .long 0x82552225, 0x3ff66238
215 .long 0x667f3bcd, 0x3ff6a09e
216 .long 0x3c651a2f, 0x3ff6dfb2
217 .long 0xe8ec5f74, 0x3ff71f75
218 .long 0x564267c9, 0x3ff75feb
219 .long 0x73eb0187, 0x3ff7a114
220 .long 0x36cf4e62, 0x3ff7e2f3
221 .long 0x994cce13, 0x3ff82589
222 .long 0x9b4492ed, 0x3ff868d9
223 .long 0x422aa0db, 0x3ff8ace5
224 .long 0x99157736, 0x3ff8f1ae
225 .long 0xb0cdc5e5, 0x3ff93737
226 .long 0x9fde4e50, 0x3ff97d82
227 .long 0x82a3f090, 0x3ff9c491
228 .long 0x7b5de565, 0x3ffa0c66
229 .long 0xb23e255d, 0x3ffa5503
230 .long 0x5579fdbf, 0x3ffa9e6b
231 .long 0x995ad3ad, 0x3ffae89f
232 .long 0xb84f15fb, 0x3ffb33a2
233 .long 0xf2fb5e47, 0x3ffb7f76
234 .long 0x904bc1d2, 0x3ffbcc1e
235 .long 0xdd85529c, 0x3ffc199b
236 .long 0x2e57d14b, 0x3ffc67f1
237 .long 0xdcef9069, 0x3ffcb720
238 .long 0x4a07897c, 0x3ffd072d
239 .long 0xdcfba487, 0x3ffd5818
240 .long 0x03db3285, 0x3ffda9e6
241 .long 0x337b9b5f, 0x3ffdfc97
242 .long 0xe78b3ff6, 0x3ffe502e
243 .long 0xa2a490da, 0x3ffea4af
244 .long 0xee615a27, 0x3ffefa1b
245 .long 0x5b6e4540, 0x3fff5076
246 .long 0x819e90d8, 0x3fffa7c1
247 .type L(DP_T), @object
248 ASM_SIZE_DIRECTIVE(L(DP_T))
250 .section .rodata.cst8,"aM",@progbits,8
252 L(DP_KLN2): /* double precision K/log(2) */
253 .long 0x652b82fe, 0x40571547
254 .type L(DP_KLN2), @object
255 ASM_SIZE_DIRECTIVE(L(DP_KLN2))
258 L(DP_NLN2K): /* double precision -log(2)/K */
259 .long 0xfefa39ef, 0xbf862e42
260 .type L(DP_NLN2K), @object
261 ASM_SIZE_DIRECTIVE(L(DP_NLN2K))
264 L(DP_RS): /* double precision 2^23+2^22 */
265 .long 0x00000000, 0x41680000
266 .type L(DP_RS), @object
267 ASM_SIZE_DIRECTIVE(L(DP_RS))
270 L(DP_P3): /* double precision polynomial coefficient P3 */
271 .long 0xeb78fa85, 0x3fa56420
272 .type L(DP_P3), @object
273 ASM_SIZE_DIRECTIVE(L(DP_P3))
276 L(DP_P1): /* double precision polynomial coefficient P1 */
277 .long 0x008d6118, 0x3fe00000
278 .type L(DP_P1), @object
279 ASM_SIZE_DIRECTIVE(L(DP_P1))
282 L(DP_P2): /* double precision polynomial coefficient P2 */
283 .long 0xda752d4f, 0x3fc55550
284 .type L(DP_P2), @object
285 ASM_SIZE_DIRECTIVE(L(DP_P2))
288 L(DP_P0): /* double precision polynomial coefficient P0 */
289 .long 0xffffe7c6, 0x3fefffff
290 .type L(DP_P0), @object
291 ASM_SIZE_DIRECTIVE(L(DP_P0))
295 .long 0x7f800000 /* single precision Inf */
296 .long 0 /* single precision zero */
297 .type L(SP_INF_0), @object
298 ASM_SIZE_DIRECTIVE(L(SP_INF_0))
300 .section .rodata.cst4,"aM",@progbits,4
302 L(SP_RS): /* single precision 2^23+2^22 */
304 .type L(SP_RS), @object
305 ASM_SIZE_DIRECTIVE(L(SP_RS))
308 L(SP_SMALL): /* single precision small value 2^(-100) */
310 .type L(SP_SMALL), @object
311 ASM_SIZE_DIRECTIVE(L(SP_SMALL))
314 L(SP_LARGE): /* single precision large value 2^100 */
316 .type L(SP_LARGE), @object
317 ASM_SIZE_DIRECTIVE(L(SP_LARGE))
320 L(SP_ONE): /* single precision 1.0 */
322 .type L(SP_ONE), @object
323 ASM_SIZE_DIRECTIVE(L(SP_ONE))
325 strong_alias (__ieee754_expf_sse2, __expf_finite_sse2)