1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
/* ix87 specific implementation of pow function.
Copyright (C) 1996 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with the GNU C Library; see the file COPYING.LIB. If not,
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA. */
#include <machine/asm.h>
#ifdef __ELF__
.section .rodata
#else
.text
#endif
.align ALIGNARG(4)
ASM_TYPE_DIRECTIVE(one,@object)
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
ASM_TYPE_DIRECTIVE(limit,@object)
limit: .double 0.29
ASM_SIZE_DIRECTIVE(limit)
#ifdef PIC
#define MO(op) op##@GOTOFF(%ecx)
#else
#define MO(op) op
#endif
.text
ENTRY(__ieee754_pow)
fldl 4(%esp) // x
fldl 12(%esp) // y : x
#ifdef PIC
call 1f
1: popl %ecx
addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
#endif
subl $8,%esp
/* First see whether `y' is a natural number. In this case we
can use a more precise algorithm. */
fld %st // y : y : x
fistpll (%esp) // y : x
fildll (%esp) // int(y) : y : x
fucomp %st(1) // y : x
fnstsw
sahf
jne 2f
/* OK, we have an integer value for y. */
ftst // y : x
fstp %st(0) // x
fnstsw
sahf
popl %eax
popl %edx
jnc 4f // y >= 0, jump
fdivrl MO(one) // 1/x (now referred to as x)
negl %eax
adcl $0, %edx
negl %edx
4: fldl MO(one) // 1 : x
fxch
6: shrdl $1, %edx, %eax
jnc 5f
fxch
fmul %st(1) // x : ST*x
fxch
5: fmul %st(0), %st // x*x : ST*x
movl %eax, %ecx
orl %edx, %ecx
jnz 6b
fstp %st(0) // ST*x
ret
.align ALIGNARG(4)
2: /* y is a real number. */
fxch // x : y
fldl MO(one) // 1.0 : x : y
fld %st(1) // x : 1.0 : x : y
fsub %st(1) // x-1 : 1.0 : x : y
fabs // |x-1| : 1.0 : x : y
fcompl MO(limit) // 1.0 : x : y
fnstsw
fxch // x : 1.0 : y
sahf
ja 7f
fsub %st(1) // x-1 : 1.0 : y
fyl2xp1 // log2(x) : y
jmp 8f
7: fyl2x // log2(x) : y
8: fmul %st(1) // y*log2(x) : y
fst %st(1) // y*log2(x) : y*log2(x)
frndint // int(y*log2(x)) : y*log2(x)
fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
fxch // fract(y*log2(x)) : int(y*log2(x))
f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
addl $8, %esp
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
ret
END(__ieee754_pow)
|
