From 96b0b1bd078020834648b5ab83b23248e4217ef6 Mon Sep 17 00:00:00 2001
From: Jessica Hamilton <jessica.l.hamilton@gmail.com>
Date: Mon, 2 Jun 2014 12:57:35 +1200
Subject: [PATCH] libroot.so: update glibc's e_pow.S on x86. Fixes #9962
---
src/system/libroot/posix/glibc/arch/x86/e_pow.S | 213 +++++++++++++++++-------
1 file changed, 150 insertions(+), 63 deletions(-)
diff --git a/src/system/libroot/posix/glibc/arch/x86/e_pow.S b/src/system/libroot/posix/glibc/arch/x86/e_pow.S
index 997dd30..76be3d1 100644
a
|
b
|
|
1 | 1 | /* ix87 specific implementation of pow function. |
2 | | Copyright (C) 1996, 1997, 1998, 1999, 2001 Free Software Foundation, Inc. |
| 2 | Copyright (C) 1996-2014 Free Software Foundation, Inc. |
| 3 | This file is part of the GNU C Library. |
| 4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. |
| 5 | |
| 6 | The GNU C Library is free software; you can redistribute it and/or |
| 7 | modify it under the terms of the GNU Lesser General Public |
| 8 | License as published by the Free Software Foundation; either |
| 9 | /* ix87 specific implementation of pow function. |
| 10 | Copyright (C) 1996-2014 Free Software Foundation, Inc. |
3 | 11 | This file is part of the GNU C Library. |
4 | 12 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. |
5 | 13 | |
… |
… |
|
14 | 22 | Lesser General Public License for more details. |
15 | 23 | |
16 | 24 | You should have received a copy of the GNU Lesser General Public |
17 | | License along with the GNU C Library; if not, write to the Free |
18 | | Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA |
19 | | 02111-1307 USA. */ |
| 25 | License along with the GNU C Library; if not, see |
| 26 | <http://www.gnu.org/licenses/>. */ |
20 | 27 | |
21 | 28 | #include <machine/asm.h> |
22 | 29 | |
23 | | #ifdef __ELF__ |
24 | | .section .rodata |
25 | | #else |
26 | | .text |
27 | | #endif |
| 30 | .section .rodata.cst8,"aM",@progbits,8 |
28 | 31 | |
29 | | .align ALIGNARG(4) |
30 | | ASM_TYPE_DIRECTIVE(infinity,@object) |
| 32 | .p2align 3 |
| 33 | .type one,@object |
| 34 | one: .double 1.0 |
| 35 | ASM_SIZE_DIRECTIVE(one) |
| 36 | .type limit,@object |
| 37 | limit: .double 0.29 |
| 38 | ASM_SIZE_DIRECTIVE(limit) |
| 39 | .type p63,@object |
| 40 | p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43 |
| 41 | ASM_SIZE_DIRECTIVE(p63) |
| 42 | .type p10,@object |
| 43 | p10: .byte 0, 0, 0, 0, 0, 0, 0x90, 0x40 |
| 44 | ASM_SIZE_DIRECTIVE(p10) |
| 45 | |
| 46 | .section .rodata.cst16,"aM",@progbits,16 |
| 47 | |
| 48 | .p2align 3 |
| 49 | .type infinity,@object |
31 | 50 | inf_zero: |
32 | 51 | infinity: |
33 | 52 | .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f |
34 | 53 | ASM_SIZE_DIRECTIVE(infinity) |
35 | | ASM_TYPE_DIRECTIVE(zero,@object) |
| 54 | .type zero,@object |
36 | 55 | zero: .double 0.0 |
37 | 56 | ASM_SIZE_DIRECTIVE(zero) |
38 | | ASM_TYPE_DIRECTIVE(minf_mzero,@object) |
| 57 | .type minf_mzero,@object |
39 | 58 | minf_mzero: |
40 | 59 | minfinity: |
41 | 60 | .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff |
42 | 61 | mzero: |
43 | 62 | .byte 0, 0, 0, 0, 0, 0, 0, 0x80 |
44 | 63 | ASM_SIZE_DIRECTIVE(minf_mzero) |
45 | | ASM_TYPE_DIRECTIVE(one,@object) |
46 | | one: .double 1.0 |
47 | | ASM_SIZE_DIRECTIVE(one) |
48 | | ASM_TYPE_DIRECTIVE(limit,@object) |
49 | | limit: .double 0.29 |
50 | | ASM_SIZE_DIRECTIVE(limit) |
51 | 64 | |
52 | 65 | #ifdef PIC |
53 | | #define MO(op) op##@GOTOFF(%ecx) |
54 | | #define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) |
| 66 | # define MO(op) op##@GOTOFF(%ecx) |
| 67 | # define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) |
55 | 68 | #else |
56 | | #define MO(op) op |
57 | | #define MOX(op,x,f) op(,x,f) |
| 69 | # define MO(op) op |
| 70 | # define MOX(op,x,f) op(,x,f) |
58 | 71 | #endif |
59 | 72 | |
60 | 73 | .text |
… |
… |
ENTRY(__ieee754_pow)
|
63 | 76 | fxam |
64 | 77 | |
65 | 78 | #ifdef PIC |
66 | | call 1f |
67 | | 1: popl %ecx |
68 | | addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx |
| 79 | LOAD_PIC_REG (cx) |
69 | 80 | #endif |
70 | 81 | |
71 | 82 | fnstsw |
… |
… |
ENTRY(__ieee754_pow)
|
74 | 85 | cmpb $0x40, %ah // is y == 0 ? |
75 | 86 | je 11f |
76 | 87 | |
77 | | cmpb $0x05, %ah // is y == ±inf ? |
| 88 | cmpb $0x05, %ah // is y == ±inf ? |
78 | 89 | je 12f |
79 | 90 | |
80 | 91 | cmpb $0x01, %ah // is y == NaN ? |
… |
… |
ENTRY(__ieee754_pow)
|
89 | 100 | movb %ah, %dh |
90 | 101 | andb $0x45, %ah |
91 | 102 | cmpb $0x40, %ah |
92 | | je 20f // x is ±0 |
| 103 | je 20f // x is ±0 |
93 | 104 | |
94 | 105 | cmpb $0x05, %ah |
95 | | je 15f // x is ±inf |
| 106 | je 15f // x is ±inf |
| 107 | |
| 108 | cmpb $0x01, %ah |
| 109 | je 32f // x is NaN |
96 | 110 | |
97 | 111 | fxch // y : x |
98 | 112 | |
| 113 | /* fistpll raises invalid exception for |y| >= 1L<<63. */ |
| 114 | fld %st // y : y : x |
| 115 | fabs // |y| : y : x |
| 116 | fcompl MO(p63) // y : x |
| 117 | fnstsw |
| 118 | sahf |
| 119 | jnc 2f |
| 120 | |
99 | 121 | /* First see whether `y' is a natural number. In this case we |
100 | 122 | can use a more precise algorithm. */ |
101 | 123 | fld %st // y : y : x |
… |
… |
ENTRY(__ieee754_pow)
|
104 | 126 | fucomp %st(1) // y : x |
105 | 127 | fnstsw |
106 | 128 | sahf |
107 | | jne 2f |
| 129 | jne 3f |
108 | 130 | |
109 | | /* OK, we have an integer value for y. */ |
| 131 | /* OK, we have an integer value for y. If large enough that |
| 132 | errors may propagate out of the 11 bits excess precision, use |
| 133 | the algorithm for real exponent instead. */ |
| 134 | fld %st // y : y : x |
| 135 | fabs // |y| : y : x |
| 136 | fcompl MO(p10) // y : x |
| 137 | fnstsw |
| 138 | sahf |
| 139 | jnc 2f |
110 | 140 | popl %eax |
111 | 141 | popl %edx |
112 | 142 | orl $0, %edx |
… |
… |
ENTRY(__ieee754_pow)
|
132 | 162 | fstp %st(0) // ST*x |
133 | 163 | ret |
134 | 164 | |
135 | | /* y is ±NAN */ |
| 165 | /* y is ±NAN */ |
136 | 166 | 30: fldl 4(%esp) // x : y |
137 | 167 | fldl MO(one) // 1.0 : x : y |
138 | 168 | fucomp %st(1) // x : y |
… |
… |
ENTRY(__ieee754_pow)
|
143 | 173 | 31: fstp %st(1) |
144 | 174 | ret |
145 | 175 | |
| 176 | 32: addl $8, %esp |
| 177 | fstp %st(1) |
| 178 | ret |
| 179 | |
| 180 | .align ALIGNARG(4) |
| 181 | 2: // y is a large integer (absolute value at least 1L<<10), but |
| 182 | // may be odd unless at least 1L<<64. So it may be necessary |
| 183 | // to adjust the sign of a negative result afterwards. |
| 184 | fxch // x : y |
| 185 | fabs // |x| : y |
| 186 | fxch // y : x |
146 | 187 | .align ALIGNARG(4) |
147 | | 2: /* y is a real number. */ |
| 188 | 3: /* y is a real number. */ |
148 | 189 | fxch // x : y |
149 | 190 | fldl MO(one) // 1.0 : x : y |
150 | | fld %st(1) // x : 1.0 : x : y |
151 | | fsub %st(1) // x-1 : 1.0 : x : y |
152 | | fabs // |x-1| : 1.0 : x : y |
153 | | fcompl MO(limit) // 1.0 : x : y |
| 191 | fldl MO(limit) // 0.29 : 1.0 : x : y |
| 192 | fld %st(2) // x : 0.29 : 1.0 : x : y |
| 193 | fsub %st(2) // x-1 : 0.29 : 1.0 : x : y |
| 194 | fabs // |x-1| : 0.29 : 1.0 : x : y |
| 195 | fucompp // 1.0 : x : y |
154 | 196 | fnstsw |
155 | 197 | fxch // x : 1.0 : y |
156 | 198 | sahf |
… |
… |
ENTRY(__ieee754_pow)
|
168 | 210 | f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) |
169 | 211 | faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) |
170 | 212 | fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) |
171 | | addl $8, %esp |
172 | 213 | fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) |
| 214 | testb $2, %dh |
| 215 | jz 292f |
| 216 | // x is negative. If y is an odd integer, negate the result. |
| 217 | fldl 20(%esp) // y : abs(result) |
| 218 | fld %st // y : y : abs(result) |
| 219 | fabs // |y| : y : abs(result) |
| 220 | fcompl MO(p63) // y : abs(result) |
| 221 | fnstsw |
| 222 | sahf |
| 223 | jnc 291f |
| 224 | |
| 225 | // We must find out whether y is an odd integer. |
| 226 | fld %st // y : y : abs(result) |
| 227 | fistpll (%esp) // y : abs(result) |
| 228 | fildll (%esp) // int(y) : y : abs(result) |
| 229 | fucompp // abs(result) |
| 230 | fnstsw |
| 231 | sahf |
| 232 | jne 292f |
| 233 | |
| 234 | // OK, the value is an integer, but is it odd? |
| 235 | popl %eax |
| 236 | popl %edx |
| 237 | andb $1, %al |
| 238 | jz 290f // jump if not odd |
| 239 | // It's an odd integer. |
| 240 | fchs |
| 241 | 290: ret |
| 242 | 291: fstp %st(0) // abs(result) |
| 243 | 292: addl $8, %esp |
173 | 244 | ret |
174 | 245 | |
175 | 246 | |
176 | | // pow(x,±0) = 1 |
| 247 | // pow(x,±0) = 1 |
177 | 248 | .align ALIGNARG(4) |
178 | 249 | 11: fstp %st(0) // pop y |
179 | 250 | fldl MO(one) |
180 | 251 | ret |
181 | 252 | |
182 | | // y == ±inf |
| 253 | // y == ±inf |
183 | 254 | .align ALIGNARG(4) |
184 | 255 | 12: fstp %st(0) // pop y |
185 | | fldl 4(%esp) // x |
186 | | fabs |
187 | | fcompl MO(one) // < 1, == 1, or > 1 |
| 256 | fldl MO(one) // 1 |
| 257 | fldl 4(%esp) // x : 1 |
| 258 | fabs // abs(x) : 1 |
| 259 | fucompp // < 1, == 1, or > 1 |
188 | 260 | fnstsw |
189 | 261 | andb $0x45, %ah |
190 | 262 | cmpb $0x45, %ah |
… |
… |
ENTRY(__ieee754_pow)
|
208 | 280 | ret |
209 | 281 | |
210 | 282 | .align ALIGNARG(4) |
211 | | // x is ±inf |
| 283 | // x is ±inf |
212 | 284 | 15: fstp %st(0) // y |
213 | 285 | testb $2, %dh |
214 | 286 | jz 16f // jump if x == +inf |
215 | 287 | |
| 288 | // fistpll raises invalid exception for |y| >= 1L<<63, so test |
| 289 | // that (in which case y is certainly even) before testing |
| 290 | // whether y is odd. |
| 291 | fld %st // y : y |
| 292 | fabs // |y| : y |
| 293 | fcompl MO(p63) // y |
| 294 | fnstsw |
| 295 | sahf |
| 296 | jnc 16f |
| 297 | |
216 | 298 | // We must find out whether y is an odd integer. |
217 | 299 | fld %st // y : y |
218 | 300 | fistpll (%esp) // y |
… |
… |
ENTRY(__ieee754_pow)
|
222 | 304 | sahf |
223 | 305 | jne 17f |
224 | 306 | |
225 | | // OK, the value is an integer, but is the number of bits small |
226 | | // enough so that all are coming from the mantissa? |
| 307 | // OK, the value is an integer. |
227 | 308 | popl %eax |
228 | 309 | popl %edx |
229 | 310 | andb $1, %al |
230 | 311 | jz 18f // jump if not odd |
231 | | movl %edx, %eax |
232 | | orl %edx, %edx |
233 | | jns 155f |
234 | | negl %eax |
235 | | 155: cmpl $0x00200000, %eax |
236 | | ja 18f // does not fit in mantissa bits |
237 | 312 | // It's an odd integer. |
238 | 313 | shrl $31, %edx |
239 | 314 | fldl MOX(minf_mzero, %edx, 8) |
… |
… |
ENTRY(__ieee754_pow)
|
256 | 331 | ret |
257 | 332 | |
258 | 333 | .align ALIGNARG(4) |
259 | | // x is ±0 |
| 334 | // x is ±0 |
260 | 335 | 20: fstp %st(0) // y |
261 | 336 | testb $2, %dl |
262 | 337 | jz 21f // y > 0 |
263 | 338 | |
264 | | // x is ±0 and y is < 0. We must find out whether y is an odd integer. |
| 339 | // x is ±0 and y is < 0. We must find out whether y is an odd integer. |
265 | 340 | testb $2, %dh |
266 | 341 | jz 25f |
267 | 342 | |
| 343 | // fistpll raises invalid exception for |y| >= 1L<<63, so test |
| 344 | // that (in which case y is certainly even) before testing |
| 345 | // whether y is odd. |
| 346 | fld %st // y : y |
| 347 | fabs // |y| : y |
| 348 | fcompl MO(p63) // y |
| 349 | fnstsw |
| 350 | sahf |
| 351 | jnc 25f |
| 352 | |
268 | 353 | fld %st // y : y |
269 | 354 | fistpll (%esp) // y |
270 | 355 | fildll (%esp) // int(y) : y |
… |
… |
ENTRY(__ieee754_pow)
|
273 | 358 | sahf |
274 | 359 | jne 26f |
275 | 360 | |
276 | | // OK, the value is an integer, but is the number of bits small |
277 | | // enough so that all are coming from the mantissa? |
| 361 | // OK, the value is an integer. |
278 | 362 | popl %eax |
279 | 363 | popl %edx |
280 | 364 | andb $1, %al |
281 | 365 | jz 27f // jump if not odd |
282 | | cmpl $0xffe00000, %edx |
283 | | jbe 27f // does not fit in mantissa bits |
284 | 366 | // It's an odd integer. |
285 | 367 | // Raise divide-by-zero exception and get minus infinity value. |
286 | 368 | fldl MO(one) |
… |
… |
ENTRY(__ieee754_pow)
|
296 | 378 | ret |
297 | 379 | |
298 | 380 | .align ALIGNARG(4) |
299 | | // x is ±0 and y is > 0. We must find out whether y is an odd integer. |
| 381 | // x is ±0 and y is > 0. We must find out whether y is an odd integer. |
300 | 382 | 21: testb $2, %dh |
301 | 383 | jz 22f |
302 | 384 | |
| 385 | // fistpll raises invalid exception for |y| >= 1L<<63, so test |
| 386 | // that (in which case y is certainly even) before testing |
| 387 | // whether y is odd. |
| 388 | fcoml MO(p63) // y |
| 389 | fnstsw |
| 390 | sahf |
| 391 | jnc 22f |
| 392 | |
303 | 393 | fld %st // y : y |
304 | 394 | fistpll (%esp) // y |
305 | 395 | fildll (%esp) // int(y) : y |
… |
… |
ENTRY(__ieee754_pow)
|
308 | 398 | sahf |
309 | 399 | jne 23f |
310 | 400 | |
311 | | // OK, the value is an integer, but is the number of bits small |
312 | | // enough so that all are coming from the mantissa? |
| 401 | // OK, the value is an integer. |
313 | 402 | popl %eax |
314 | 403 | popl %edx |
315 | 404 | andb $1, %al |
316 | 405 | jz 24f // jump if not odd |
317 | | cmpl $0xffe00000, %edx |
318 | | jae 24f // does not fit in mantissa bits |
319 | 406 | // It's an odd integer. |
320 | 407 | fldl MO(mzero) |
321 | 408 | ret |