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) 19962014 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) 19962014 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   021111307 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) // x1 : 1.0 : x : y 
152   fabs // x1 : 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) // x1 : 0.29 : 1.0 : x : y 
 194  fabs // x1 : 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 dividebyzero 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 